A comprehensive quality-control and quantification RNA-seq pipeline

public 1yr ago Version: v1.9.0 0 bookmarks

View Workflow

Help improve this workflow!

This workflow has been published but could be further improved with some additional meta data:

Keyword(s) in categories input, output, operation, topic

You can help improve this workflow by suggesting the addition or removal of keywords, suggest changes and report issues, or request to become a maintainer of the Workflow .

An open-source, reproducible, and scalable solution for analyzing RNA-seq data.

Introduction
Overview
2.1 RNA-seek Pipeline
2.2 Reference Genomes
2.3 Dependencies
2.4 Installation
Run RNA-seek pipeline
3.1 Using Singularity
3.2 Using Docker
3.3 Biowulf
Contribute
References

1. Introduction

RNA-sequencing ( RNA-seq ) has a wide variety of applications. This popular transcriptome profiling technique can be used to quantify gene and isoform expression, detect alternative splicing events, predict gene-fusions, call variants and much more.

RNA-seek is a comprehensive, open-source RNA-seq pipeline that relies on technologies like Docker20 and Singularity21 to maintain the highest-level of reproducibility. The pipeline consists of a series of data processing and quality-control steps orchestrated by Snakemake19 , a flexible and scalable workflow management system, to submit jobs to a cluster or cloud provider.

RNA-seek_overview_diagram
Fig 1. Run locally on a compute instance, on-premise using a cluster, or on the cloud using AWS. A user can define the method or mode of execution. The pipeline can submit jobs to a cluster using a job scheduler like SLURM, or run on AWS using Tibanna (feature coming soon!). A hybrid approach ensures the pipeline is accessible to all users. As an optional step, relevelant output files and metadata can be stored in object storage using HPC DME (NIH users) or Amazon S3 for archival purposes (coming soon!).

2. Overview

2.1 RNA-seek Pipeline

A bioinformatics pipeline is more than the sum of its data processing steps. A pipeline without quality-control steps provides a myopic view of the potential sources of variation within your data (i.e., biological verses technical sources of variation). RNA-seek pipeline is composed of a series of quality-control and data processing steps.

The accuracy of the downstream interpretations made from transcriptomic data are highly dependent on initial sample library. Unwanted sources of technical variation, which if not accounted for properly, can influence the results. RNA-seek's comprehensive quality-control helps ensure your results are reliable and reproducible across experiments . In the data processing steps, RNA-seek quantifies gene and isoform expression and predicts gene fusions. Please note that the detection of alternative splicing events and variant calling will be incorporated in a later release.

RNA-seq quantification pipeline Fig 2. An Overview of RNA-seek Pipeline. Gene and isoform counts are quantified and a series of QC-checks are performed to assess the quality of the data. This pipeline stops at the generation of a raw counts matrix and gene-fusion calling. To run the pipeline, a user must select their raw data, a reference genome, and output directory (i.e., the location where the pipeline performs the analysis). Quality-control information is summarized across all samples in a MultiQC report.

Quality Control
FastQC 2 is used to assess the sequencing quality. FastQC is run twice, before and after adapter trimming. It generates a set of basic statistics to identify problems that can arise during sequencing or library preparation. FastQC will summarize per base and per read QC metrics such as quality scores and GC content. It will also summarize the distribution of sequence lengths and will report the presence of adapter sequences.

Kraken2 14 and FastQ Screen 17 are used to screen for various sources of contamination. During the process of sample collection to library preparation, there is a risk for introducing wanted sources of DNA. FastQ Screen compares your sequencing data to a set of different reference genomes to determine if there is contamination. It allows a user to see if the composition of your library matches what you expect. Also, if there are high levels of microbial contamination, Kraken can provide an estimation of the taxonomic composition. Kraken can be used in conjunction with Krona 15 to produce interactive reports.

Preseq 1 is used to estimate the complexity of a library for each samples. If the duplication rate is very high, the overall library complexity will be low. Low library complexity could signal an issue with library preparation where very little input RNA was over-amplified or the sample may be degraded.

Picard 10 can be used to estimate the duplication rate, and it has another particularly useful sub-command called CollectRNAseqMetrics which reports the number and percentage of reads that align to various regions: such as coding, intronic, UTR, intergenic and ribosomal regions. This is particularly useful as you would expect a library constructed with ploy(A)-selection to have a high percentage of reads that map to coding regions. Picard CollectRNAseqMetrics will also report the uniformity of coverage across all genes, which is useful for determining whether a sample has a 3' bias (observed in ploy(A)-selection libraries containing degraded RNA).

RSeQC 9 is another particularity useful package that is tailored for RNA-seq data. It is used to calculate the inner distance between paired-end reads and calculate TIN values for a set of canonical protein-coding transcripts. A median TIN value is calucated for each sample, which analogous to a computationally derived RIN.

MultiQC11 is used to aggreate the results of each tool into a single interactive report.

Quantification
Cutadapt 3 is used to remove adapter sequences, perform quality trimming, and remove very short sequences that would otherwise multi-map all over the genome prior to alignment.

STAR 4 is used to align reads to the reference genome. The RNA-seek pipeline runs STAR in a two-passes where splice-junctions are collected and aggregated across all samples and provided to the second-pass of STAR. In the second pass of STAR, the splice-junctions detected in the first pass are inserted into the genome indices prior to alignment.

RSEM 5 is used to quantify gene and isoform expression. The expected counts from RSEM are merged across samples to create a two counts matrices for gene counts and isoform counts.

Arriba 22 is used to predict gene-fusion events. The pre-built human and mouse reference genomes use Arriba blacklists to reduce the false-positive rate.

2.2 Reference Genomes

Reference files are pulled from an S3 bucket to the compute instance or local filesystem prior to execution.
RNA-seek comes bundled with pre-built reference files for the following genomes:

Name	Species	Genome	Annotation
hg38_30	Homo sapiens (human)	GRCh38	Gencode6 Release 30
mm10_M21	Mus musculus (mouse)	GRCm38	Gencode Release M21

Warning: This section contains FTP links for downloading each reference file. Open the link in a new tab to start a download.

2.3 Dependencies

Requires: singularity>=3.5 snakemake>=6.0

Snakemake and singularity must be installed on the target system. Snakemake orchestrates the execution of each step in the pipeline. To guarantee reproducibility, each step relies on pre-built images from DockerHub . Snakemake uses singaularity to pull these images onto the local filesystem prior to job execution, and as so, snakemake and singularity are the only two dependencies.

2.4 Installation

Please clone this repository to your local filesystem using the following command:

# Clone Repository from Github
git clone https://github.com/skchronicles/RNA-seek.git
# Change your working directory to the RNA-seek repo
cd RNA-seek/

3. Run RNA-seek pipeline

3.1 Using Singularity

# Coming Soon!

3.2 Using Docker

# Coming Soon!

3.3 Biowulf

# rna-seek is configured to use different execution backends: local or slurm
# view the help page for more information
./rna-seek run --help
# @local: uses local singularity execution method
# The local MODE will run serially on compute
# instance. This is useful for testing, debugging,
# or when a users does not have access to a high
# performance computing environment.
# Please note that you can dry-run the command below
# by providing the --dry-run flag
# Do not run this on the head node!
# Grab an interactive node
sinteractive --mem=110g --cpus-per-task=12 --gres=lscratch:200
module purge
module load singularity snakemake
./rna-seek run --input .tests/*.R?.fastq.gz --output /data/$USER/RNA_hg38 --genome hg38_30 --mode local
# @slurm: uses slurm and singularity execution method
# The slurm MODE will submit jobs to the cluster.
# It is recommended running rna-seek in this mode.
module purge
module load singularity snakemake
./rna-seek run --input .tests/*.R?.fastq.gz --output /data/$USER/RNA_hg38 --genome hg38_30 --mode slurm

4. Contribute

This section is for new developers working with the RNA-seek pipeline. If you have added new features or adding new changes, please consider contributing them back to the original repository:

Fork the original repo to a personal or org account.
Clone the fork to your local filesystem.
Copy the modified files to the cloned fork.
Commit and push your changes to your fork.
Create a pull request to this repository.

5. References

1. Daley, T. and A.D. Smith, Predicting the molecular complexity of sequencing libraries. Nat Methods, 2013. 10(4): p. 325-7.
2. Andrews, S. (2010). FastQC: a quality control tool for high throughput sequence data.
3. Martin, M. (2011). "Cutadapt removes adapter sequences from high-throughput sequencing reads." EMBnet 17(1): 10-12.
4. Dobin, A., et al., STAR: ultrafast universal RNA-seq aligner. Bioinformatics, 2013. 29(1): p. 15-21.
5. Li, B. and C.N. Dewey, RSEM: accurate transcript quantification from RNA-Seq data with or without a reference genome. BMC Bioinformatics, 2011. 12: p. 323.
6. Harrow, J., et al., GENCODE: the reference human genome annotation for The ENCODE Project. Genome Res, 2012. 22(9): p. 1760-74.
7. Law, C.W., et al., voom: Precision weights unlock linear model analysis tools for RNA-seq read counts. Genome Biol, 2014. 15(2): p. R29.
8. Smyth, G.K., Linear models and empirical bayes methods for assessing differential expression in microarray experiments. Stat Appl Genet Mol Biol, 2004. 3: p. Article3.
9. Wang, L., et al. (2012). "RSeQC: quality control of RNA-seq experiments." Bioinformatics 28(16): 2184-2185.
10. The Picard toolkit. https://broadinstitute.github.io/picard/.
11. Ewels, P., et al. (2016). "MultiQC: summarize analysis results for multiple tools and samples in a single report." Bioinformatics 32(19): 3047-3048.
12. R Core Team (2018). R: A Language and Environment for Statistical Computing. Vienna, Austria, R Foundation for Statistical Computing.
13. Li, H., et al. (2009). "The Sequence Alignment/Map format and SAMtools." Bioinformatics 25(16): 2078-2079.
14. Wood, D. E. and S. L. Salzberg (2014). "Kraken: ultrafast metagenomic sequence classification using exact alignments." Genome Biol 15(3): R46.
15. Ondov, B. D., et al. (2011). "Interactive metagenomic visualization in a Web browser." BMC Bioinformatics 12(1): 385.
16. Okonechnikov, K., et al. (2015). "Qualimap 2: advanced multi-sample quality control for high-throughput sequencing data." Bioinformatics 32(2): 292-294.
17. Wingett, S. and S. Andrews (2018). "FastQ Screen: A tool for multi-genome mapping and quality control." F1000Research 7(2): 1338.
18. Robinson, M. D., et al. (2009). "edgeR: a Bioconductor package for differential expression analysis of digital gene expression data." Bioinformatics 26(1): 139-140.
19. Koster, J. and S. Rahmann (2018). "Snakemake-a scalable bioinformatics workflow engine." Bioinformatics 34(20): 3600.
20. Merkel, D. (2014). Docker: lightweight linux containers for consistent development and deployment. Linux Journal, 2014(239), 2.
21. Kurtzer GM, Sochat V, Bauer MW (2017). Singularity: Scientific containers for mobility of compute. PLoS ONE 12(5): e0177459.
22. Haas, B. J., et al. (2019). "Accuracy assessment of fusion transcript detection via read-mapping and de novo fusion transcript assembly-based methods." Genome Biology 20(1): 213.

Code Snippets

shell: """
python {params.get_flowcell_lanes} {input.R1} {wildcards.name} > {output.fqinfo}
"""

SnakeMake From line 36 of rules/common.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

java -Xmx{params.memory}g  -XX:ParallelGCThreads={threads} -jar ${{PICARDJARPATH}}/picard.jar AddOrReplaceReadGroups \
    I={input.file1} O=${{tmp}}/{params.sampleName}.star_rg_added.sorted.bam \
    TMP_DIR=${{tmp}} RGID=id RGLB=library RGPL=illumina RGPU=machine RGSM=sample VALIDATION_STRINGENCY=SILENT;
java -Xmx{params.memory}g -XX:ParallelGCThreads={threads} -jar ${{PICARDJARPATH}}/picard.jar MarkDuplicates \
    I=${{tmp}}/{params.sampleName}.star_rg_added.sorted.bam \
    O=${{tmp}}/{params.sampleName}.star_rg_added.sorted.dmark.bam \
    TMP_DIR=${{tmp}} CREATE_INDEX=true VALIDATION_STRINGENCY=SILENT METRICS_FILE={output.metrics};

mv ${{tmp}}/{params.sampleName}.star_rg_added.sorted.dmark.bam {output.bam};
mv ${{tmp}}/{params.sampleName}.star_rg_added.sorted.dmark.bai {output.bai};
sed -i 's/MarkDuplicates/picard.sam.MarkDuplicates/g' {output.metrics};
"""

SnakeMake Picard From line 63 of rules/common.smk

shell:"""
preseq c_curve -B -o {output.ccurve} {input.bam}
"""

SnakeMake preseq From line 103 of rules/common.smk

shell: """
unset DISPLAY;
qualimap bamqc -bam {input.bamfile} --feature-file {params.gtfFile} \
    -outdir {params.outdir} -nt {threads} --java-mem-size={params.memory}G
"""

SnakeMake QualiMap From line 130 of rules/common.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

java -Xmx{params.memory}g -jar ${{PICARDJARPATH}}/picard.jar CollectRnaSeqMetrics \
    REF_FLAT={params.refflat} I={input.file1} O={output.outstar1} \
    RIBOSOMAL_INTERVALS={params.rrnalist} \
    STRAND_SPECIFICITY=SECOND_READ_TRANSCRIPTION_STRAND \
    TMP_DIR=${{tmp}}  VALIDATION_STRINGENCY=SILENT;
sed -i 's/CollectRnaSeqMetrics/picard.analysis.CollectRnaSeqMetrics/g' {output.outstar1}
samtools flagstat {input.file1} > {output.outstar2};
python3 {params.statscript} {input.file1} >> {output.outstar2}
"""

SnakeMake SAMtools Picard From line 168 of rules/common.smk

shell: """
python {params.pythonscript} {params.annotate} {params.inputdir} {params.inputdir}
sed 's/\\t/|/1' {output.gene_counts_matrix} | \
    sed '1 s/^gene_id|GeneName/symbol/' > {output.reformatted}
"""

SnakeMake From line 210 of rules/common.smk

shell: """
infer_experiment.py -r {params.bedref} -i {input.file1} -s 1000000 > {output.out1}
read_distribution.py -i {input.file1} -r {params.bedref} > {output.out2}
"""

SnakeMake From line 236 of rules/common.smk

shell: """
# tin.py writes to current working directory
cd {params.outdir}
tin.py -i {input.bam} -r {params.bedref}
"""

SnakeMake From line 266 of rules/common.smk

shell: """
python {params.create_matrix} {input.tins} > {output.matrix}
"""

SnakeMake From line 291 of rules/common.smk

shell: """
mkdir -p {params.outdir}
fastQValidator --noeof --file {input.R1} > {output.out1}
fastQValidator --noeof --file {input.R2} > {output.out2}
"""

SnakeMake From line 30 of rules/paired-end.smk

shell: """
fastqc {input.R1} {input.R2} -t {threads} -o {params.outdir};
"""

SnakeMake FastQC From line 60 of rules/paired-end.smk

shell: """
cutadapt --pair-filter=any --nextseq-trim=2 --trim-n \
    -n 5 -O 5 -q {params.leadingquality},{params.trailingquality} \
    -m {params.minlen}:{params.minlen} \
    -b file:{params.fastawithadaptersetd} -B file:{params.fastawithadaptersetd} \
    -j {threads} -o {output.out1} -p {output.out2} {input.file1} {input.file2}
"""

SnakeMake Cutadapt From line 93 of rules/paired-end.smk

shell: """
fastqc {input.R1} {input.R2} -t {threads} -o {params.outdir};
"""

SnakeMake FastQC From line 126 of rules/paired-end.smk

shell: """
# Get encoding of Phred Quality Scores
encoding=$(python {params.encoding} {input.R1})
echo "Detected Phred+${{encoding}} ASCII encoding"

bbtools bbmerge-auto in1={input.R1} in2={input.R2} qin=${{encoding}} \
    ihist={output} k=62 extend2=200 rem ecct -Xmx{params.memory}G
"""

SnakeMake From line 154 of rules/paired-end.smk

shell: """
fastq_screen --conf {params.fastq_screen_config} --outdir {params.outdir} \
    --threads {threads} --subset 1000000 \
    --aligner bowtie2 --force {input.file1} {input.file2}

fastq_screen --conf {params.fastq_screen_config2} --outdir {params.outdir2} \
    --threads {threads} --subset 1000000 \
    --aligner bowtie2 --force {input.file1} {input.file2}
"""

SnakeMake Bowtie 2 From line 201 of rules/paired-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Copy kraken2 db to /lscratch or temp 
# location to reduce filesystem strain
cp -rv {params.bacdb} ${{tmp}}/;
kdb_base=$(basename {params.bacdb})
kraken2 --db ${{tmp}}/${{kdb_base}} \
    --threads {threads} --report {output.krakentaxa} \
    --output {output.krakenout} \
    --gzip-compressed \
    --paired {input.fq1} {input.fq2}
# Generate Krona Report
cut -f2,3 {output.krakenout} | \
    ktImportTaxonomy - -o {output.kronahtml}
"""

SnakeMake kraken2 Krona From line 240 of rules/paired-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(
    zcat {input.file1} | \
    awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; \
    END {{print maxlen-1}}'
)

echo "sjdbOverhang for STAR: ${{readlength}}"

STAR --genomeDir {params.stardir} \
    --outFilterIntronMotifs {params.filterintronmotifs} \
    --outSAMstrandField {params.samstrandfield}  \
    --outFilterType {params.filtertype} \
    --outFilterMultimapNmax {params.filtermultimapnmax} \
    --alignSJoverhangMin {params.alignsjoverhangmin} \
    --alignSJDBoverhangMin {params.alignsjdboverhangmin} \
    --outFilterMismatchNmax {params.filtermismatchnmax} \
    --outFilterMismatchNoverLmax {params.filtermismatchnoverlmax} \
    --alignIntronMin {params.alignintronmin} \
    --alignIntronMax {params.alignintronmax} \
    --alignMatesGapMax {params.alignmatesgapmax} \
    --clip3pAdapterSeq {params.adapter1} {params.adapter2} \
    --readFilesIn {input.file1} {input.file2} \
    --readFilesCommand zcat \
    --runThreadN {threads} \
    --outFileNamePrefix {params.prefix}. \
    --outSAMunmapped {params.outsamunmapped} \
    --outWigType {params.wigtype} \
    --outWigStrand {params.wigstrand} \
    --twopassMode Basic \
    --sjdbGTFfile {params.gtffile} \
    --limitSjdbInsertNsj {params.nbjuncs} \
    --quantMode TranscriptomeSAM GeneCounts \
    --outSAMtype BAM Unsorted \
    --alignEndsProtrude 10 ConcordantPair \
    --peOverlapNbasesMin 10 \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --sjdbOverhang ${{readlength}}

# SAMtools sort (uses less memory than STAR SortedByCoordinate)
samtools sort -@ {threads} \
    -m 2G -T ${{tmp}}/SORTtmp_{wildcards.name} \
    -O bam {params.prefix}.Aligned.out.bam \
    > {output.out1}

rm {params.prefix}.Aligned.out.bam
mv {params.prefix}.Aligned.toTranscriptome.out.bam {workpath}/{bams_dir};
mv {params.prefix}.Log.final.out {workpath}/{log_dir}
"""

SnakeMake SAMtools STAR From line 331 of rules/paired-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(
    zcat {input.file1} | \
    awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; \
    END {{print maxlen-1}}'
)

echo "sjdbOverhang for STAR: ${{readlength}}"

STAR --genomeDir {params.stardir} \
    --outFilterIntronMotifs {params.filterintronmotifs} \
    --outSAMstrandField {params.samstrandfield} \
    --outFilterType {params.filtertype} \
    --outFilterMultimapNmax {params.filtermultimapnmax} \
    --alignSJoverhangMin {params.alignsjoverhangmin} \
    --alignSJDBoverhangMin {params.alignsjdboverhangmin} \
    --outFilterMismatchNmax {params.filtermismatchnmax} \
    --outFilterMismatchNoverLmax {params.filtermismatchnoverlmax} \
    --alignIntronMin {params.alignintronmin} \
    --alignIntronMax {params.alignintronmax} \
    --alignMatesGapMax {params.alignmatesgapmax} \
    --clip3pAdapterSeq {params.adapter1} {params.adapter2} \
    --readFilesIn {input.file1} {input.file2} \
    --readFilesCommand zcat \
    --runThreadN {threads} \
    --outFileNamePrefix {params.prefix}. \
    --outSAMtype BAM Unsorted \
    --alignEndsProtrude 10 ConcordantPair \
    --peOverlapNbasesMin 10 \
    --sjdbGTFfile {params.gtffile} \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --sjdbOverhang ${{readlength}}
"""

SnakeMake STAR From line 433 of rules/paired-end.smk

shell: """
cat {input.files} | \
    sort | \
    uniq | \
    awk -F \"\\t\" '{{if ($5>0 && $6==1) {{print}}}}'| \
    cut -f1-4 | sort | uniq | \
grep \"^chr\" | grep -v \"^chrM\" > {output.out1}
"""

SnakeMake From line 492 of rules/paired-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(
    zcat {input.file1} | \
    awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; \
    END {{print maxlen-1}}'
)

echo "sjdbOverhang for STAR: ${{readlength}}"

STAR --genomeDir {params.stardir} \
    --outFilterIntronMotifs {params.filterintronmotifs} \
    --outSAMstrandField {params.samstrandfield}  \
    --outFilterType {params.filtertype} \
    --outFilterMultimapNmax {params.filtermultimapnmax} \
    --alignSJoverhangMin {params.alignsjoverhangmin} \
    --alignSJDBoverhangMin {params.alignsjdboverhangmin} \
    --outFilterMismatchNmax {params.filtermismatchnmax} \
    --outFilterMismatchNoverLmax {params.filtermismatchnoverlmax} \
    --alignIntronMin {params.alignintronmin} \
    --alignIntronMax {params.alignintronmax} \
    --alignMatesGapMax {params.alignmatesgapmax} \
    --clip3pAdapterSeq {params.adapter1} {params.adapter2} \
    --readFilesIn {input.file1} {input.file2} \
    --readFilesCommand zcat \
    --runThreadN {threads} \
    --outFileNamePrefix {params.prefix}. \
    --outSAMunmapped {params.outsamunmapped} \
    --outWigType {params.wigtype} \
    --outWigStrand {params.wigstrand} \
    --sjdbFileChrStartEnd {input.tab} \
    --sjdbGTFfile {params.gtffile} \
    --limitSjdbInsertNsj {params.nbjuncs} \
    --quantMode TranscriptomeSAM GeneCounts \
    --outSAMtype BAM Unsorted \
    --alignEndsProtrude 10 ConcordantPair \
    --peOverlapNbasesMin 10 \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --sjdbOverhang ${{readlength}}

# SAMtools sort (uses less memory than STAR SortedByCoordinate)
samtools sort -@ {threads} \
    -m 2G -T ${{tmp}}/SORTtmp_{wildcards.name} \
    -O bam {params.prefix}.Aligned.out.bam \
    > {output.out1}

rm {params.prefix}.Aligned.out.bam
mv {params.prefix}.Aligned.toTranscriptome.out.bam {workpath}/{bams_dir};
mv {params.prefix}.Log.final.out {workpath}/{log_dir}
"""

SnakeMake SAMtools STAR From line 554 of rules/paired-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(
    zcat {input.R1} | \
    awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; \
    END {{print maxlen-1}}'
)

# Avoids inheriting $R_LIBS_SITE
# from local env variables
R_LIBS_SITE=/usr/local/lib/R/site-library

STAR --runThreadN {threads} \
    --sjdbGTFfile {params.gtffile} \
    --sjdbOverhang ${{readlength}} \
    --genomeDir {params.stardir} \
    --genomeLoad NoSharedMemory \
    --readFilesIn {input.R1} {input.R2} \
    --readFilesCommand zcat \
    --outStd BAM_Unsorted \
    --outSAMtype BAM Unsorted \
    --outSAMunmapped Within \
    --outFilterMultimapNmax 50 \
    --peOverlapNbasesMin 10 \
    --alignSplicedMateMapLminOverLmate 0.5 \
    --alignSJstitchMismatchNmax 5 -1 5 5 \
    --chimSegmentMin 10 \
    --chimOutType WithinBAM HardClip \
    --chimJunctionOverhangMin 10 \
    --chimScoreDropMax 30 \
    --chimScoreJunctionNonGTAG 0 \
    --chimScoreSeparation 1 \
    --chimSegmentReadGapMax 3 \
    --chimMultimapNmax 50 \
    --twopassMode Basic \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --outFileNamePrefix {params.prefix}. \
| tee ${{tmp}}/{params.chimericbam} | \
arriba -x /dev/stdin \
    -o {output.fusions} \
    -O {output.discarded} \
    -a {params.reffa} \
    -g {params.gtffile} \
    -b {input.blacklist} \

# Sorting and Indexing BAM files is required for Arriba's Visualization
samtools sort -@ {threads} \
    -m 2G -T ${{tmp}}/SORTtmp_{wildcards.name} \
    -O bam ${{tmp}}/{params.chimericbam} \
    > {output.bam}

samtools index {output.bam} {output.bai}
rm ${{tmp}}/{params.chimericbam}

# Generate Gene Fusions Visualization
draw_fusions.R \
    --fusions={output.fusions} \
    --alignments={output.bam} \
    --output={output.figure} \
    --annotation={params.gtffile} \
    --cytobands={input.cytoband} \
    --proteinDomains={input.protdomain}
"""

SnakeMake SAMtools STAR Arriba From line 652 of rules/paired-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Get strandedness to calculate Forward Probability
fp=$(tail -n1 {input.file2} | awk '{{if($NF > 0.75) print "0.0"; else if ($NF<0.25) print "1.0"; else print "0.5";}}')

echo "Forward Probability Passed to RSEM: $fp"
rsem-calculate-expression --no-bam-output --calc-ci --seed 12345  \
    --bam --paired-end -p {threads}  {input.file1} {params.rsemref} {params.prefix} --time \
    --temporary-folder ${{tmp}} --keep-intermediate-files --forward-prob=${{fp}} --estimate-rspd
"""

SnakeMake RSEM From line 750 of rules/paired-end.smk

shell: """
inner_distance.py -i {input.bam} -r {params.genemodel} \
    -k 10000000 -o {params.prefix}
"""

SnakeMake From line 786 of rules/paired-end.smk

shell: """
bash {params.bashscript} {input.bam} {params.outprefix}

# reverse files if method is not dUTP/NSR/NNSR ... ie, R1 in the direction of RNA strand.
strandinfo=`tail -n1 {input.strandinfo} | awk '{{print $NF}}'`
if [ `echo "$strandinfo < 0.25"|bc` -eq 1 ];then
    mv {output.fbw} {output.fbw}.tmp
    mv {output.rbw} {output.fbw}
    mv {output.fbw}.tmp {output.rbw}
fi
"""

SnakeMake From line 818 of rules/paired-end.smk

shell: """
multiqc --ignore '*/.singularity/*' -f -c {params.qcconfig} --interactive --outdir {params.outdir} {params.workdir}

# Parse RSeQC Inner Distance Maximas
echo -e "Sample\\tInner_Dist_Maxima" > {output.maximas}
for f in $(find {params.workdir} -iname '*.inner_distance_freq.txt'); do
    sample=$(basename "${{f}}");
    inner_dist_maxima=$(sort -k3,3nr "${{f}}" | awk -F '\\t' 'NR==1{{print $1}}');
    echo -e "${{sample}}\\t${{inner_dist_maxima}}";
done >> {output.maximas}

# Parse RSeQC Median TINs
echo -e "Sample\\tmedian_tin" > {output.medtins}
find {params.workdir} -name '*.star_rg_added.sorted.dmark.summary.txt' -exec cut -f1,3 {{}} \\; | \
    grep -v '^Bam_file' | \
    awk -F '\\t' '{{printf "%s\\t%.3f\\n", $1,$2}}' >> {output.medtins}

# Parse Flowcell and Lane information
echo -e "Sample\\tflowcell_lanes" > {output.fclanes}
find {params.workdir} -name '*.fastq.info.txt' -exec awk -F '\\t' -v OFS='\\t' 'NR==2 {{print $1,$5}}' {{}} \\; \
    >> {output.fclanes}

python3 {params.pyparser} {params.logfiles} {params.outdir}
"""

SnakeMake MultiQC From line 871 of rules/paired-end.smk

shell: """
# Avoids inheriting $R_LIBS_SITE
# from local env variables
R_LIBS_SITE=/usr/local/lib/R/site-library
# Generate RNA QC Dashboard
{params.rwrapper} \
    -m {params.rmarkdown} \
    -r {input.counts} \
    -t {input.tins} \
    -q {input.qc} \
    -o {params.odir} \
    -f RNA_Report.html
"""

SnakeMake From line 926 of rules/paired-end.smk

shell: """
mkdir -p {params.outdir}
fastQValidator --noeof --file {input.R1} > {output.out1}
"""

SnakeMake From line 28 of rules/single-end.smk

shell: """
fastqc {input.R1} -t {threads} -o {params.outdir};
"""

SnakeMake FastQC From line 55 of rules/single-end.smk

shell: """
cutadapt --nextseq-trim=2 --trim-n \
    -n 5 -O 5 -q {params.leadingquality},{params.trailingquality} \
    -m 16 -b file:{params.fastawithadaptersetd} -j {threads} \
    -o {output.outfq} {input.infq}
"""

SnakeMake Cutadapt From line 86 of rules/single-end.smk

shell: """
cutadapt --nextseq-trim=2 --trim-n \
    -n 5 -O 5 -q {params.leadingquality},{params.trailingquality} \
    -m {params.minlen} -b file:{params.fastawithadaptersetd} -j {threads} \
    -o {output.outfq} {input.infq}
"""

SnakeMake Cutadapt From line 118 of rules/single-end.smk

shell: """
fastqc {input} -t {threads} -o {params.outdir};
"""

SnakeMake FastQC From line 148 of rules/single-end.smk

shell: """
fastq_screen --conf {params.fastq_screen_config} --outdir {params.outdir} \
    --threads {threads} --subset 1000000 --aligner bowtie2 --force {input.file1}

fastq_screen --conf {params.fastq_screen_config2} --outdir {params.outdir2} \
    --threads {threads} --subset 1000000 --aligner bowtie2 --force {input.file1}
"""

SnakeMake Bowtie 2 From line 184 of rules/single-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Copy kraken2 db to /lscratch or temp 
# location to reduce filesytem strain
cp -rv {params.bacdb} ${{tmp}}/;
kdb_base=$(basename {params.bacdb})
kraken2 --db ${{tmp}}/${{kdb_base}} \
    --threads {threads} --report {output.krakentaxa} \
    --output {output.krakenout} \
    --gzip-compressed \
    {input.fq}
# Generate Krona Report
cut -f2,3 {output.krakenout} | \
    ktImportTaxonomy - -o {output.kronahtml}
"""

SnakeMake kraken2 Krona From line 220 of rules/single-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(zcat {input.file1} | awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; END {{print maxlen-1}}')
echo "sjdbOverhang for STAR: ${{readlength}}"

STAR --genomeDir {params.stardir} \
    --outFilterIntronMotifs {params.filterintronmotifs} \
    --outSAMstrandField {params.samstrandfield} \
    --outFilterType {params.filtertype} \
    --outFilterMultimapNmax {params.filtermultimapnmax} \
    --alignSJoverhangMin {params.alignsjoverhangmin} \
    --alignSJDBoverhangMin {params.alignsjdboverhangmin} \
    --outFilterMismatchNmax {params.filtermismatchnmax} \
    --outFilterMismatchNoverLmax {params.filtermismatchnoverlmax} \
    --alignIntronMin {params.alignintronmin} \
    --alignIntronMax {params.alignintronmax} \
    --alignMatesGapMax {params.alignmatesgapmax} \
    --clip3pAdapterSeq {params.adapter1} {params.adapter2} \
    --readFilesIn {input.file1} --readFilesCommand zcat \
    --runThreadN {threads} \
    --outFileNamePrefix {params.prefix}. \
    --outSAMunmapped {params.outsamunmapped} \
    --outWigType {params.wigtype} \
    --outWigStrand {params.wigstrand} \
    --twopassMode Basic \
    --sjdbGTFfile {params.gtffile} \
    --limitSjdbInsertNsj {params.nbjuncs} \
    --quantMode TranscriptomeSAM GeneCounts \
    --outSAMtype BAM SortedByCoordinate \
    --alignEndsProtrude 10 ConcordantPair \
    --peOverlapNbasesMin 10 \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --sjdbOverhang ${{readlength}}

mv {params.prefix}.Aligned.toTranscriptome.out.bam {workpath}/{bams_dir};
mv {params.prefix}.Log.final.out {workpath}/{log_dir}
"""

SnakeMake STAR From line 309 of rules/single-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(zcat {input.file1} | awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; END {{print maxlen-1}}')
echo "sjdbOverhang for STAR: ${{readlength}}"

STAR --genomeDir {params.stardir} \
    --outFilterMultimapNmax {params.filtermultimapnmax} \
    --alignSJDBoverhangMin 1000 \
    --outFilterScoreMinOverLread 0 \
    --outFilterMatchNminOverLread 0 \
    --outFilterMatchNmin 16 \
    --outFilterMismatchNoverLmax {params.filtermismatchnoverlmax} \
    --alignIntronMax 1 \
    --clip3pAdapterSeq {params.adapter1} {params.adapter2} \
    --readFilesIn {input.file1} --readFilesCommand zcat \
    --runThreadN {threads} \
    --outFileNamePrefix {params.prefix}. \
    --outSAMunmapped Within \
    --sjdbGTFfile {params.gtffile} \
    --limitSjdbInsertNsj {params.nbjuncs} \
    --quantMode TranscriptomeSAM GeneCounts \
    --outSAMtype BAM Unsorted \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --sjdbOverhang ${{readlength}}

# SAMtools sort (uses less memory than STAR SortedByCoordinate)
samtools sort -@ {threads} \
    -m 2G -T ${{tmp}}/SORTtmp_{wildcards.name} \
    -O bam {params.prefix}.Aligned.out.bam \
    > {output.out1}

rm {params.prefix}.Aligned.out.bam
mv {params.prefix}.Aligned.toTranscriptome.out.bam {workpath}/{bams_dir};
mv {params.prefix}.Log.final.out {workpath}/{log_dir}
"""

SnakeMake SAMtools STAR From line 403 of rules/single-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(zcat {input.file1} | awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; END {{print maxlen-1}}')
echo "sjdbOverhang for STAR: ${{readlength}}"

STAR --genomeDir {params.stardir} \
    --outFilterIntronMotifs {params.filterintronmotifs} \
    --outSAMstrandField {params.samstrandfield}  \
    --outFilterType {params.filtertype} \
    --outFilterMultimapNmax {params.filtermultimapnmax} \
    --alignSJoverhangMin {params.alignsjoverhangmin} \
    --alignSJDBoverhangMin {params.alignsjdboverhangmin} \
    --outFilterMismatchNmax {params.filtermismatchnmax} \
    --outFilterMismatchNoverLmax {params.filtermismatchnoverlmax} \
    --alignIntronMin {params.alignintronmin} \
    --alignIntronMax {params.alignintronmax} \
    --alignMatesGapMax {params.alignmatesgapmax} \
    --clip3pAdapterSeq {params.adapter1} {params.adapter2} \
    --readFilesIn {input.file1} --readFilesCommand zcat \
    --runThreadN {threads} \
    --outFileNamePrefix {params.prefix}. \
    --outSAMtype BAM Unsorted \
    --alignEndsProtrude 10 ConcordantPair \
    --peOverlapNbasesMin 10 \
    --sjdbGTFfile {params.gtffile} \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --sjdbOverhang ${{readlength}}
"""

SnakeMake STAR From line 489 of rules/single-end.smk

shell: """
cat {input.files} | \
    sort | uniq | \
    awk -F '\\t' '{{if ($5>0 && $6==1) {{print}}}}'| \
    cut -f1-4 | sort | uniq | \
grep \"^chr\"|grep -v \"^chrM\" > {output.out1}
"""

SnakeMake From line 542 of rules/single-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Optimal readlength for sjdbOverhang = max(ReadLength) - 1 [Default: 100]
readlength=$(zcat {input.file1} | awk -v maxlen=100 'NR%4==2 {{if (length($1) > maxlen+0) maxlen=length($1)}}; END {{print maxlen-1}}')
echo "sjdbOverhang for STAR: ${{readlength}}"

STAR --genomeDir {params.stardir} \
    --outFilterIntronMotifs {params.filterintronmotifs} \
    --outSAMstrandField {params.samstrandfield} \
    --outFilterType {params.filtertype} \
    --outFilterMultimapNmax {params.filtermultimapnmax} \
    --alignSJoverhangMin {params.alignsjoverhangmin} \
    --alignSJDBoverhangMin {params.alignsjdboverhangmin} \
    --outFilterMismatchNmax {params.filtermismatchnmax} \
    --outFilterMismatchNoverLmax {params.filtermismatchnoverlmax} \
    --alignIntronMin {params.alignintronmin} \
    --alignIntronMax {params.alignintronmax} \
    --alignMatesGapMax {params.alignmatesgapmax} \
    --clip3pAdapterSeq {params.adapter1} {params.adapter2} \
    --readFilesIn {input.file1} --readFilesCommand zcat \
    --runThreadN {threads} \
    --outFileNamePrefix {params.prefix}. \
    --outSAMunmapped {params.outsamunmapped} \
    --outWigType {params.wigtype} \
    --outWigStrand {params.wigstrand} \
    --sjdbFileChrStartEnd {input.tab} \
    --sjdbGTFfile {params.gtffile} \
    --limitSjdbInsertNsj {params.nbjuncs} \
    --quantMode TranscriptomeSAM GeneCounts \
    --outSAMtype BAM SortedByCoordinate \
    --alignEndsProtrude 10 ConcordantPair \
    --peOverlapNbasesMin 10 \
    --outTmpDir=${{tmp}}/STARtmp_{wildcards.name} \
    --sjdbOverhang ${{readlength}}

mv {params.prefix}.Aligned.toTranscriptome.out.bam {workpath}/{bams_dir};
mv {params.prefix}.Log.final.out {workpath}/{log_dir}
"""

SnakeMake STAR From line 602 of rules/single-end.smk

shell: """
# Setups temporary directory for
# intermediate files with built-in 
# mechanism for deletion on exit
if [ ! -d "{params.tmpdir}" ]; then mkdir -p "{params.tmpdir}"; fi
tmp=$(mktemp -d -p "{params.tmpdir}")
trap 'rm -rf "${{tmp}}"' EXIT

# Get strandedness to calculate Forward Probability
fp=$(tail -n1 {input.file2} | awk '{{if($NF > 0.75) print "0.0"; else if ($NF<0.25) print "1.0"; else print "0.5";}}')

echo "Forward Probability Passed to RSEM: $fp"
rsem-calculate-expression --no-bam-output --calc-ci --seed 12345 \
    --bam -p {threads}  {input.file1} {params.rsemref} {params.prefix} --time \
    --temporary-folder ${{tmp}} --keep-intermediate-files --forward-prob=${{fp}} --estimate-rspd
"""

SnakeMake RSEM From line 674 of rules/single-end.smk

shell: """
bash {params.bashscript} {input.bam} {params.outprefix}

# reverse files if method is not dUTP/NSR/NNSR ... ie, R1 in the direction of RNA strand.
strandinfo=`tail -n1 {input.strandinfo} | awk '{{print $NF}}'`
if [ `echo "$strandinfo < 0.25"|bc` -eq 1 ]; then
    mv {output.fbw} {output.fbw}.tmp
    mv {output.rbw} {output.fbw}
    mv {output.fbw}.tmp {output.rbw}
fi
"""

SnakeMake From line 718 of rules/single-end.smk

shell: """
multiqc --ignore '*/.singularity/*' -f -c {params.qcconfig} --interactive --outdir {params.outdir} {params.workdir}

# Parse RSeQC Median TINs
echo -e "Sample\\tmedian_tin" > {output.medtins}
find {params.workdir} -name '*.star_rg_added.sorted.dmark.summary.txt' -exec cut -f1,3 {{}} \\; | \
    grep -v '^Bam_file' | \
    awk -F '\\t' '{{printf "%s\\t%.3f\\n", $1,$2}}' >> {output.medtins}

# Parse Flowcell and Lane information
echo -e "Sample\\tflowcell_lanes" > {output.fclanes}
find {params.workdir} -name '*.fastq.info.txt' -exec awk -F '\\t' -v OFS='\\t' 'NR==2 {{print $1,$5}}' {{}} \\; \
    >> {output.fclanes}

python3 {params.pyparser} {params.logfiles} {params.outdir}
"""