Articles in this section

See more

In Detail: General RNA-Seq analysis pipeline (STAR)

  • Updated

Description

This workflow processes Illumina RNA-Seq sequencing reads by cleaning up reads, aligning them to a reference genome and quantifying gene expression. Analysis steps thus include trimming (BBDuk), mapping (STAR), expression quantification (STAR) and normalization (rnanorm). Reference genomes for mouse, rat, or human (Ensembl, version 109, 100 and 92) are optionally 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 interfere with the analysis. 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

The main pipeline steps, tool versions and parameters used are listed below. Pipeline details for both single-end (SE) and paired-end (PE) library preparation types are shared to account for diverse input data options.

 

Reads trimming using BBDuk (BBMap 37.90, Bushnell, 2018)

bbduk.sh (SE) in=input_reads.fastq.gz ref=input_references.fasta.gz out=output_reads.fastq stats=output_statistics.txt statscolumns=5 k=23 rcomp=t maskmiddle=t minkmerhits=1 minkmerfraction=0.0 mincovfraction=0.0 hammingdistance=1 qhdist=0 editdistance=0 hammingdistance2=0 qhdist2=0 editdistance2=0 forbidn=f findbestmatch=t ktrim=r maskfullycovered=f mink=11 qtrim=r trimq=28 trimpolya=0 minlength=30 minlengthfraction=0.0 minavgquality=0 maqb=0 minbasequality=0 minconsecutivebases=0 trimpad=0 minoverlap=14 mininsert=40 forcetrimleft=0 forcetrimright=0 forcetrimright2=0 forcetrimmod=0 restrictleft=0 restrictright=0 mingc=0.0 maxgc=1.0 maxns=-1 tossjunk=f chastityfilter=f barcodefilter=f xmin=-1 ymin=-1 xmax=-1 ymax=-1 entropy=-1.0 entropywindow=50 entropyk=5 minbasefrequency=0 entropymask=f threads=4
bbduk.sh (PE) in=input_fw_reads.fastq.gz in2=input_rv_reads.fastq.gz ref=input_references.fasta.gz out=output_fw_reads.fastq out2=output_rv_reads.fastq stats=output_statistics.txt statscolumns=5 k=23 rcomp=t maskmiddle=t minkmerhits=1 minkmerfraction=0.0 mincovfraction=0.0 hammingdistance=1 qhdist=0 editdistance=0 hammingdistance2=0 qhdist2=0 editdistance2=0 forbidn=f removeifeitherbad=t findbestmatch=t ecco=f ktrim=r maskfullycovered=f mink=11 qtrim=r trimq=28 trimpolya=0 minlength=30 minlengthfraction=0.0 minavgquality=0 maqb=0 minbasequality=0 minconsecutivebases=0 trimpad=0 trimbyoverlap=t minoverlap=14 mininsert=40 trimpairsevenly=t forcetrimleft=0 forcetrimright=0 forcetrimright2=0 forcetrimmod=0 restrictleft=0 restrictright=0 mingc=0.0 maxgc=1.0 maxns=-1 tossjunk=f chastityfilter=f barcodefilter=f xmin=-1 ymin=-1 xmax=-1 ymax=-1 entropy=-1.0 entropywindow=50 entropyk=5 minbasefrequency=0 entropymask=f threads=4

Read alignment using STAR (v2.7.10b, Dobin et al. 2013)

Reference genome: Ensembl, version 109)

STAR --runThreadN 1 --runMode genomeGenerate --genomeSAsparseD 5 --genomeDir ./star_index/ --genomeFastaFiles genome.fasta --sjdbGTFfile annotation.gtf --sjdbOverhang 100 --sjdbGTFfeatureExon exon
STAR (SE) --runThreadN 4 --genomeDir /star_index --readFilesIn mate1.fastq.gz --readFilesCommand zcat --limitIObufferSize 30000000 50000000 --limitOutSAMoneReadBytes 100000 --limitOutSJcollapsed 1000000 --limitOutSJoneRead 1000 --limitSjdbInsertNsj 1000000 --outReadsUnmapped Fastx --outSAMtype BAM Unsorted --outSAMattributes Standard --outSAMunmapped Within --outFilterType Normal --alignMatesGapMax 0 --quantMode GeneCounts --twopassMode Basic
STAR (PE) --runThreadN 4 --genomeDir /star_index --readFilesIn mate1.fastq mate2.fastq --readFilesCommand zcat --limitIObufferSize 30000000 50000000      --limitOutSAMoneReadBytes 100000 --limitOutSJcollapsed 1000000 --limitOutSJoneRead 1000 --limitSjdbInsertNsj 1000000 --outReadsUnmapped Fastx --outSAMtype BAM Unsorted      --outSAMattributes Standard --outSAMunmapped Within --outFilterType Normal --alignMatesGapMax 0 --quantMode GeneCounts --twopassMode Basic

Sorting of BAM files using Samtools (v1.14, Danecek et. al., 2021)

Samtools sort Aligned.out.bam -o sorted.bam -@ 4

Normalization of expression values using rnanorm (1.3.1)

rnanorm raw_counts.txt --tpm-output tpm_expression_values.txt --cpm-output cpm_expression_values.txt --annotation annotation.gtf

Seqtk (1.2-r94)

Subsampling to 1M reads, using a fixed random seed

seqtk (SE)sample -s 11 input_reads.fastq.gz 1000000 > output_reads.fastq
Seqtk (PE)sample -s 11 input_mate1_reads.fastq.gz 1000000 > mate1_out.fastq
seqtk (PE)sample -s 11 input_mate2_reads.fastq.gz 1000000 > mate2_out.fastq

 

References

Related articles