4.4.1 RNA-seq¶

This section runs the alignment and quantification pipeline on the cleaned RNA-seq reads. The commands are the same as those in one_step/run_step2.sh.

Step 1: Classify reads with Bowtie¶

Note: This step is not strictly required for RNA-seq, because RNA-seq libraries are usually constructed with oligo(dT) priming, which selectively enriches polyadenylated mRNA. It becomes necessary only when the library has been prepared using rRNA depletion, in which case the remaining rRNA and other contaminating reads must be removed before genome alignment.

Reads are aligned sequentially against the rRNA, tRNA, ncRNA, mRNA, and genome references. Reads mapped to a database are retained (--al) for later inspection, while the unmapped reads (--un) are passed on to the next database. In this way, contaminating reads are removed step by step and only the reads of interest are kept for genome alignment.

1.1 Create directories¶

Create the 2.bowtie directory that will hold all classification results, and enter it.

mkdir -p ./sce/3.rna-seq/2.bowtie/
cd ./sce/3.rna-seq/2.bowtie/

1.2 Classify reads against sequential references¶

For each sample, the reads are aligned against the five references in turn, allowing at most $mismatch mismatches (-v $mismatch).

# set database
rrna='../../1.reference/rrna/rrna'
trna='../../1.reference/trna/trna'
ncrna='../../1.reference/ncrna/ncrna'
mrna='../../1.reference/mrna/mrna'
chrom='../../1.reference/genome/genome'
threads=12
mismatch=1

# align reads to the reference databases
for fq in ../1.cleandata/*fastq.gz
do
  fqname=$(basename $fq .clean.fastq.gz)

  # rRNA
  bowtie -p $threads -v $mismatch --un="$fqname".norrna.fq --al="$fqname".rrna.fq \
    -x $rrna $fq -S "$fqname".rrna.sam 2>> "$fqname".log

  # tRNA
  bowtie -p $threads -v $mismatch --un="$fqname".notrna.fq --al="$fqname".trna.fq \
    -x $trna "$fqname".norrna.fq -S "$fqname".trna.sam 2>> "$fqname".log

  # ncRNA
  bowtie -p $threads -v $mismatch --un="$fqname".noncrna.fq --al="$fqname".ncrna.fq \
    -x $ncrna "$fqname".notrna.fq -S "$fqname".ncrna.sam 2>> "$fqname".log

  # mRNA
  bowtie -p $threads -v $mismatch --un="$fqname".nomrna.fq --al="$fqname".mrna.fq \
    -x $mrna "$fqname".noncrna.fq -S "$fqname".mrna.sam 2>> "$fqname".log

  # genome
  bowtie -p $threads -v $mismatch --un="$fqname".nogenome.fq --al="$fqname".genome.fq \
    -x $chrom "$fqname".nomrna.fq -S "$fqname".genome.sam 2>> "$fqname".log

  # compress FASTQ files
  pigz *fq

  # compress SAM files
  for sam in *.sam
  do
    samtools view -h -F 4 $sam | samtools sort -@ $threads -o $(basename $sam sam)bam
    rm $sam
  done

done

Step 2: Merge Bowtie mapping statistics¶

merge_bwt_log parses the Bowtie log files generated in Step 1 and merges the mapping statistics of all samples into a single table, which can then be inspected to check the read composition of each sample.

2.1 Parameters¶

Parameter	Meaning
`-l`, `--list`	Bowtie log files, such as `*log`
`-o`	output prefix
`-n`, `--name`	comma-separated database labels, such as `rRNA,tRNA,ncRNA,mRNA,Genome`

2.2 Example¶

merge_bwt_log \
  -n rRNA,tRNA,ncRNA,mRNA,Genome \
  -l *log \
  -o RNA_seq \
  &>> merge_bowtie.log

2.3 Output¶

File	Description
`RNA_seq_mapping.txt`	merged mapping statistics
`RNA_seq_mapping_barplot.png` / `RNA_seq_mapping_barplot.pdf`	stacked barplot of read classification across databases
`merge_bowtie.log`	running log

The stacked barplot summarizes the proportion of reads assigned to each database for every sample: refer to 4.4.2 Ribo-seq for details.

Step 3: Align mRNA reads with STAR¶

Reads that fail to map to any reference (*.noncrna.fq.gz) are aligned to the genome with STAR.

mkdir -p ./sce/3.rna-seq/3.star/

cd ./sce/3.rna-seq/3.star/

genome='../../1.reference/star-index/'
threads=12

for fastq in ../2.bowtie/*.noncrna.fq.gz
do
  output=$(basename $fastq .noncrna.fq.gz)

  STAR --runThreadN $threads \
    --readFilesCommand zcat \
    --genomeDir $genome \
    --readFilesIn $fastq \
    --outFileNamePrefix $output \
    --outSAMtype BAM Unsorted \
    --outFilterType BySJout \
    --quantMode TranscriptomeSAM GeneCounts \
    --outReadsUnmapped Fastx \
    --outSAMattributes All \
    --alignEndsType Local \
    --outFilterMultimapNmax 3 \
    --outFilterMismatchNmax 1 \
    --alignIntronMax 10000 \
    --outFilterMatchNmin 20

  # compress unmapped reads
  pigz *mate1

  # sort and index the BAM file
  samtools sort -@ $threads \
    $output"Aligned.out.bam" \
    -o $output"Aligned.sortedByCoord.out.bam"
  samtools index -@ $threads $output"Aligned.sortedByCoord.out.bam"
  rm $output"Aligned.out.bam"

done

Step 4: Quantify expression with RSEM¶

RSEM estimates gene and isoform expression from the transcriptome-aligned BAM files produced by STAR. For each sample it outputs a .genes.results and an .isoforms.results file, containing expected_count, TPM, and FPKM values for every gene or isoform.

mkdir -p ./sce/3.rna-seq/4.quantification/

cd ./sce/3.rna-seq/4.quantification/

for bam in ../3.star/*Aligned.toTranscriptome.out.bam
do
  rsem-calculate-expression -p 12 \
    --no-bam-output \
    --alignments \
    -q $bam \
    ../../1.reference/rsem-index/rsem \
    $(basename $bam Aligned.toTranscriptome.out.bam)
done

Step 5: Merge RSEM results¶

merge_rsem merges a selected column (for example expected_count, TPM, or FPKM) from multiple RSEM result files into a single table, in which each row is a gene or isoform and each column is a sample.

5.1 Parameters¶

Parameter	Meaning
`-l`, `--list`	RSEM result files, such as `.genes.results` or `.isoforms.results`
`-o`	output file
`-c`, `--column`	one of `expected_count`, `TPM`, or `FPKM`

5.2 Examples¶

# merge gene-level expression
merge_rsem -c expected_count -l *.genes.results -o gene.expected_count.txt
merge_rsem -c TPM -l *.genes.results -o gene.TPM.txt
merge_rsem -c FPKM -l *.genes.results -o gene.FPKM.txt

# merge isoform-level expression
merge_rsem -c expected_count -l *.isoforms.results -o isoforms.expected_count.txt
merge_rsem -c TPM -l *.isoforms.results -o isoforms.TPM.txt
merge_rsem -c FPKM -l *.isoforms.results -o isoforms.FPKM.txt

5.3 Output¶

File	Description
`gene.expected_count.txt` / `gene.TPM.txt` / `gene.FPKM.txt`	merged gene-level quantification tables
`isoforms.expected_count.txt` / `isoforms.TPM.txt` / `isoforms.FPKM.txt`	merged isoform-level quantification tables