In Detail: General 3’ mRNA-Seq pipeline

Description

Analysis steps include trimming (BBDuk or cutadapt), mapping (STAR), and expression quantification (featureCounts). Reference genomes for mouse, rat, or human (Ensembl, version 100 and 92) are extended with references for common spike-in standards (ERCC and SIRV). You can examine “expressions” of internal spike-in controls in the same way as interrogating expression of any other gene, but you can rest assured that the presence of spike-ins in the reference does not compromise the analysis even if spike-in controls have not been used. As an additional quality control step, a sample of a million reads (Seqtk tool) is mapped to rRNA and globin sequences of the selected species to determine the overall proportion of these kinds of reads in the sample. Results are reported in the summary table of the MultiQC report.

Pipeline Details - Tools and Parameters

Unless stated otherwise, all parameter values are set to tool defaults. Example command-line calls to the major analysis tools utilized in the pipeline are listed below.

Reads trimming using Cutadapt (v2.4, M. Martin, 2011)

cutadapt \
-m 20 \
-O 20 \
-n 2 \
-a “polyA=A{20}” \
-a “QUALITY=G{20}” \
<input.fastq.gz> \
| cutadapt \
-m 20 \
--nextseq-trim=10 \
-a “truseq=A{18}AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC” \
- \
| cutadapt \
-m 20 \
-O 20 \
-g “truseq=A{18}AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC” \
--discard-trimmed \
-o <input_trimmed.fastq.gz> -

By default, the Cutadapt trimming step applies the option --nextseq-trim, which is to be used with the data from the recent Illumina machines that utilize 2-color chemistry to encode the four bases. If the data is obtained with older illumina machines, the analysis protocol can be configured so that Cutadapt uses --quality-cutoff option instead. To change the Cutadapt read trimming mode from the set default, please toggle the Advanced options section of the pipeline input parameters, and set the Quality cutoff option to the desired value.

Read alignment using STAR (2.7.0f, Dobin et al. 2013)

STAR \
--readFilesCommand zcat \
--genomeDir <STAR index> \
--limitBAMsort RAM 10000000000 \
--readFilesIn <input_trimmed.fastq.gz> \
--outFilterType BySJout \
--outFilterMultimapNmax 20 \
--alignSJoverhangMin 8 \
--alignSJDBoverhangMin 1 \
--outFilterMismatchNmax 999 \
--outFilterMismatchNoverLmax 0.6 \
--alignIntronMin 20 \
--alignIntronMax 1000000 \
--alignMatesGapMax 1000000 \
--outSAMattributes NH HI NM MD \
--outSAMstrandField intronMotif \
--outSAMtype BAM SortedByCoordinate \

Quantification using featureCounts (1.6.3, Liao et al. 2014)

featureCounts \
-s 1 \
-t exon \
-g gene_id \
-a <annotation.gtf> \
-o <featureCounts_results_file>
-M <aligned_reads.bam>

QC Steps

Normalization of Expression Values

A Custom script is used for calculating final expression values (TPM) from the featureCounts-derived raw read count results.

Downsampling of input reads using Seqtk (1.2-r94)

References