RNA-seq expression#

Miniprotocol Timing#

Timing <3.5 hours

Overview#

This miniprotocol shows the use of various modules to prepare reference data, perform RNA-seq calling, quantify expression, conduct quality control and normalize data. The modules are as follows:

  1. RNA_calling.ipynb (steps i-v): Quantifying expression from RNA-seq data

  2. bulk_expression_QC.ipynb (step vi): Sample level RNA-seq quality control

  3. bulk_expression_normalization.ipynb (step vii): Bulk RNA-seq counts normalization

Steps#

i. Perform data quality summary via fastqc#

sos run pipeline/RNA_calling.ipynb fastqc \
    --cwd output/rnaseq/fastqc \
    --sample-list data/fastq.list.txt \
    --data-dir data/fastq

ii. Cut adaptor (Optional)#

sos run pipeline/RNA_calling.ipynb fastp_trim_adaptor \
    --cwd output/rnaseq --sample-list data/fastq.list.txt \
    --data-dir data/fastq --STAR-index reference_data/STAR_Index/ \
    --gtf reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf \
    --reference-fasta reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
    --ref-flat reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.ref.flat

iii. Read alignment via STAR and QC via Picard#

sos run pipeline/RNA_calling.ipynb STAR_align \
    --cwd output/rnaseq/bam --sample-list data/fastq.list.txt \
    --data-dir data/fastq --STAR-index reference_data/STAR_Index/ \
    --gtf reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf \
    --reference-fasta reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
    --ref-flat reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.ref.flat \
    --chimSegmentMin 0 \
    -J 50 --mem 200G --numThreads 8    

iv. Call gene-level RNA expression via rnaseqc#

sos run pipeline/RNA_calling.ipynb rnaseqc_call \
    --cwd data/bam \
    --sample-list data/fastq.list.txt \
    --data-dir data/fastq \
    --gtf reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.gtf \
    --reference-fasta reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
    --varVCFfile reference_data/ZOD14598_AD_GRM_WGS_2021-04-29_all.recalibrated_variants.leftnorm.filtered.AF.WASP.vcf \
    --bam_list data/bam/sample_bam_list.txt

v. Call transcript level RNA expression via RSEM#

sos run pipeline/RNA_calling.ipynb rsem_call \
    --cwd data/bam \
    --sample-list data/fastq.list.txt \
    --data-dir data/fastq \
    --STAR-index reference_data/STAR_Index/ \
    --gtf reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.gtf \
    --reference-fasta reference_data/GRCh38_full_analysis_set_plus_decoy_hla.noALT_noHLA_noDecoy_ERCC.fasta \
    --ref-flat reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.ERCC.ref.flat \
    --varVCFfile reference_data/ZOD14598_AD_GRM_WGS_2021-04-29_all.recalibrated_variants.leftnorm.filtered.AF.WASP.vcf \
    --bam_list data/bam/sample_bam_list.txt \
    --RSEM-index reference_data/RSEM_Index

vi. Multi-sample RNA-seq QC#

sos run pipeline/bulk_expression_QC.ipynb qc \
    --cwd output/rnaseq \
    --tpm-gct data/rnaseq/bulk_rnaseq_tmp_matrix.bed \
    --counts-gct data/rnaseq/bulk_rnaseq_count_matrix.bed

vii. Multi-sample read count normalization#

sos run pipeline/bulk_expression_normalization.ipynb normalize \
    --cwd output/rnaseq \
    --tpm-gct output/rnaseq/bulk_rnaseq_tmp_matrix.low_expression_filtered.outlier_removed.tpm.gct.gz \
    --counts-gct output/rnaseq/bulk_rnaseq_tmp_matrix.low_expression_filtered.outlier_removed.geneCount.gct.gz \
    --annotation-gtf reference_data/Homo_sapiens.GRCh38.103.chr.reformatted.collapse_only.gene.ERCC.gtf  \
    --count-threshold 1 --sample_participant_lookup data/rnaseq/sample_participant_lookup.txt

Anticipated Results#

The final output contained QCed and normalized expression data in a bed.gz file. This file is ready for use in TensorQTL.

Figure 1A. Bulk RNA-Seq Quality Control D-Statistic Distribution.

Figure 1B. Bulk RNA-Seq Quality Control Relative Log Expression Residuals.

Figure 1C. Bulk RNA-Seq Quality Control Mahalanobis Distance P-Value Clustering.