A Snakemake workflow to process paired-end RNA-Seq data.

public public 1yr ago 0 bookmarks
View Workflow
Report
snakernaseq
snakernaseq
snakernaseq
View Workflow

A Snakemake workflow to process paired-end RNA-Seq data.

Steps:

The workflow consists following steps:

  • Quality control of the raw and/or trimmed data (FastQC, MultiQC)

  • Adapter trimming w/ trim_galore (Optional)

  • Contamination check and decontamination (Optional)

  • Alignment to the reference genome (hisat2, STAR)

  • Alignment quality control with RSeQC, QualiMap

  • Transcript/gene quantification (StringTie, featureCounts, RSEM)

  • Alignment-free transcript quantification (kallisto/salmon)

Future additions:

  • Differential gene expression analysis (deseq2)

  • Machine learning-based mapping uncertainty analysis (GeneQC)

dag

Code Snippets

17
18
19
20
shell:
    """
    fastqc --outdir {params.prefix} -t {threads} -f fastq {input}
    """
SnakeMake FastQC From line 17 of rules/01a_fastqc.smk 
 9
10
wrapper:
    "0.35.0/bio/multiqc"
SnakeMake MultiQC From line 9 of rules/01b_multiqc.smk 
13
14
wrapper:
    "0.35.2/bio/trim_galore/pe"
SnakeMake From line 13 of rules/02_trim.smk 
30
31
32
33
shell:
    """
    fastqc -t {threads} {input.r1} {input.r2} -q -f fastq -o {params.prefix}
    """
SnakeMake FastQC From line 30 of rules/02_trim.smk 
43
44
wrapper:
    "0.49.0/bio/multiqc"
SnakeMake MultiQC From line 43 of rules/02_trim.smk 
54
55
wrapper:
    "0.49.0/bio/multiqc"
SnakeMake MultiQC From line 54 of rules/02_trim.smk 
11
12
13
14
shell:
    '''
    bwa mem -t {threads} {input.rrna} {input.r1} {input.r2} > {output}
    '''
SnakeMake BWA From line 11 of rules/03a_rrna_check.smk 
24
25
26
27
shell:
    '''
    samtools view -@ {threads} -bS -o {output} {input.sam}
    '''
SnakeMake SAMtools From line 24 of rules/03a_rrna_check.smk 
37
38
39
40
shell:
    '''
    samtools flagstat {input.bam} --threads {threads} > {output}
    '''
SnakeMake SAMtools From line 37 of rules/03a_rrna_check.smk 
50
51
52
53
shell:
    '''
    samtools stats {input.bam} > {output}
    '''
SnakeMake SAMtools From line 50 of rules/03a_rrna_check.smk 
64
65
wrapper:
    "0.48.0/bio/multiqc"
SnakeMake MultiQC From line 64 of rules/03a_rrna_check.smk 
18
19
20
21
shell:
    """
    bbsplit.sh -Xmx120g threads={threads} in1={input.r1} in2={input.r2} ref_rrna={input.rrna} path={params.index} basename={params.out_rrna}/out_%.fq outu1={output.out1} outu2={output.out2} refstats={output.stats} 2> {log}
    """
SnakeMake From line 18 of rules/03b_rrna_clean.smk 
37
38
39
40
shell:
    """
    fastqc -t {threads} {input.r1} {input.r2} -q -f fastq -o {params.prefix}
    """
SnakeMake FastQC From line 37 of rules/03b_rrna_clean.smk 
50
51
wrapper:
    "0.47.0/bio/multiqc"
SnakeMake MultiQC From line 50 of rules/03b_rrna_clean.smk 
13
14
15
16
shell:
    '''
    bwa mem -t {threads} {input.rrna} {input.r1} {input.r2} > {output}
    '''
SnakeMake BWA From line 13 of rules/03c_rrna_cleaned_check.smk 
26
27
28
29
shell:
    '''
    samtools view -@ {threads} -bS -o {output} {input.sam}
    '''
SnakeMake SAMtools From line 26 of rules/03c_rrna_cleaned_check.smk 
39
40
41
42
shell:
    '''
    samtools stats {input.bam} -@ {threads} > {output}
    '''
SnakeMake SAMtools From line 39 of rules/03c_rrna_cleaned_check.smk 
52
53
wrapper:
    "0.48.0/bio/multiqc"
SnakeMake MultiQC From line 52 of rules/03c_rrna_cleaned_check.smk 
8
9
shell:
    "gffread {input.anno} -T -o {output.gtf}"
SnakeMake gffread From line 8 of rules/04a_star.smk 
22
23
wrapper:
    "0.49.0/bio/star/index"
SnakeMake From line 22 of rules/04a_star.smk 
41
42
wrapper:
    "0.49.0/bio/star/align"
SnakeMake From line 41 of rules/04a_star.smk 
49
50
51
52
shell:
    """
    cat {input} | awk '($5 > 0 && $7 > 2 && $6==0)' | cut -f1-6 | sort | uniq > {output}
    """
SnakeMake From line 49 of rules/04a_star.smk 
73
74
wrapper:
    "0.49.0/bio/star/align"
SnakeMake From line 73 of rules/04a_star.smk 
11
12
script:
    "../scripts/gtf2bed.py"
SnakeMake From line 11 of rules/04b_rseqc.smk 
38
39
40
shell:
    "junction_annotation.py {params.extra} -i {input.bam} -r {input.bed} -o {params.prefix} "
    "> {log[0]} 2>&1"
SnakeMake From line 38 of rules/04b_rseqc.smk 
51
52
shell:
    "bam_stat.py -i {input} > {output} 2> {log}"
SnakeMake From line 51 of rules/04b_rseqc.smk 
65
66
shell:
    "infer_experiment.py -r {input.bed} -i {input.bam} > {output} 2> {log}"
SnakeMake From line 65 of rules/04b_rseqc.smk 
78
79
shell:
    "read_distribution.py -r {input.bed} -i {input.bam} > {output} 2> {log}"
SnakeMake From line 78 of rules/04b_rseqc.smk 
92
93
shell:
    "read_GC.py -i {input} -o {params.prefix} > {log} 2>&1"
SnakeMake From line 92 of rules/04b_rseqc.smk 
126
127
wrapper:
    "0.49.0/bio/multiqc"
SnakeMake MultiQC From line 126 of rules/04b_rseqc.smk 
11
12
13
14
shell:
    '''
    samtools sort -@ {threads} -n -o {output} -T {params.prfx} {input.bam}
    '''
SnakeMake SAMtools From line 11 of rules/04c_qualimap.smk 
26
27
28
29
shell:
    '''
    qualimap rnaseq -bam {input.sorted_bam} -gtf {params.gtf} --outdir {params.outdir} --sorted --paired
    '''
SnakeMake QualiMap From line 26 of rules/04c_qualimap.smk 
38
39
wrapper:
    "0.49.0/bio/multiqc"
SnakeMake MultiQC From line 38 of rules/04c_qualimap.smk 
12
13
14
15
16
17
shell:
    """
    rsem-prepare-reference \
    -p {threads} \
    --gff3 {input.gff} {input.fasta} {params.ref_name}
    """
SnakeMake RSEM From line 12 of rules/05a_rsem.smk 
31
32
33
34
35
36
37
38
shell:
    """
    rsem-calculate-expression \
    --no-qualities \
    -p {threads} \
    --strandedness reverse \
    --alignments --paired-end {input.bam} {params.genomedir} {params.prefix}
    """
SnakeMake RSEM From line 31 of rules/05a_rsem.smk 
47
48
wrapper:
    "0.49.0/bio/multiqc"
SnakeMake MultiQC From line 47 of rules/05a_rsem.smk 
13
14
15
16
shell:
    """
    htseq-count -m intersection-nonempty --stranded=reverse --idattr gene_id -r pos -f bam {input} {params.gtf} > {output} 2> {log}
    """
SnakeMake htseqcount From line 13 of rules/05b_htseq.smk 
29
30
31
32
shell:
    """
    multiqc -m htseq {params.prefix} --filename {output} 2> {log}
    """
SnakeMake MultiQC HTSeq From line 29 of rules/05b_htseq.smk 
12
13
14
15
16
17
shell:
    """
    grep "^>" {input.ref} | cut -d " " -f 1 > {output.decoy}
    sed -i.bak -e 's/>//g' {output.decoy}
    cat {input.tcp} {input.ref} > {output.gent}
    """
SnakeMake From line 12 of rules/06a_salmon.smk 
30
31
32
33
shell:
    """
    salmon index -p {threads} -t {input.gent} -d {input.decoy} -i {output} &> {log}
    """
SnakeMake Salmon From line 30 of rules/06a_salmon.smk 
49
50
51
52
shell:
    """
    salmon quant -i {input.index} -l A -1 {input.r1} -2 {input.r2} -o {output} --validateMappings --gcBias --seqBias --writeUnmappedNames --writeMappings={output.mappings} -p {threads} --numBootstraps 100
    """
SnakeMake Quant Salmon From line 49 of rules/06a_salmon.smk 
61
62
wrapper:
    "0.49.0/bio/multiqc"
SnakeMake MultiQC From line 61 of rules/06a_salmon.smk 
73
74
75
76
shell:
    '''
    gffread -w {output} -g {input.ref} {input.gtf}
    '''
SnakeMake gffread From line 73 of rules/06a_salmon.smk 
91
92
93
94
shell:
    '''
    salmon quant -t {input.tcp} -l A -a {input.bam} -o {output} --gcBias --seqBias --writeUnmappedNames -p {threads} -g {input.gtf} --numBootstraps 100
    '''
SnakeMake Quant Salmon From line 91 of rules/06a_salmon.smk 
106
107
wrapper:
    "0.47.0/bio/multiqc"
SnakeMake MultiQC From line 106 of rules/06a_salmon.smk 
11
12
shell:
    "kallisto index -i {output} {input} 2> {log}"
SnakeMake kallisto From line 11 of rules/06b_kallisto.smk 
29
30
31
shell:
    "kallisto quant --threads {threads} -i {input.idx} -o {output} --gtf {input.gtf} "
    "{params.extra} {input.r1} {input.r2} 2> {log}"
SnakeMake Quant kallisto From line 29 of rules/06b_kallisto.smk 
40
41
wrapper:
    "0.49.0/bio/multiqc"
SnakeMake MultiQC From line 40 of rules/06b_kallisto.smk 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
import gffutils

db = gffutils.create_db(snakemake.input[0],
                        dbfn=snakemake.output.db,
                        force=True,
                        keep_order=True,
                        merge_strategy='merge',
                        sort_attribute_values=True,
                        disable_infer_genes=True,
                        disable_infer_transcripts=True)

with open(snakemake.output.bed, 'w') as outfileobj:
    for tx in db.features_of_type('transcript', order_by='start'):
        bed = [s.strip() for s in db.bed12(tx).split('\t')]
        bed[3] = tx.id
        outfileobj.write('{}\n'.format('\t'.join(bed)))
Python From line 3 of scripts/gtf2bed.py 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
__author__ = "Julian de Ruiter"
__copyright__ = "Copyright 2017, Julian de Ruiter"
__email__ = "[email protected]"
__license__ = "MIT"


from os import path

from snakemake.shell import shell


input_dirs = set(path.dirname(fp) for fp in snakemake.input)
output_dir = path.dirname(snakemake.output[0])
output_name = path.basename(snakemake.output[0])
log = snakemake.log_fmt_shell(stdout=True, stderr=True)

shell(
    "multiqc"
    " {snakemake.params}"
    " --force"
    " -o {output_dir}"
    " -n {output_name}"
    " {input_dirs}"
    " {log}")
Python Snakemake MultiQC From line 3 of multiqc/wrapper.py 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
__author__ = "Kerrin Mendler"
__copyright__ = "Copyright 2018, Kerrin Mendler"
__email__ = "[email protected]"
__license__ = "MIT"


from snakemake.shell import shell
import os.path


log = snakemake.log_fmt_shell()

# Check that two input files were supplied
n = len(snakemake.input)
assert n == 2, "Input must contain 2 files. Given: %r." % n

# Don't run with `--fastqc` flag
if "--fastqc" in snakemake.params.get("extra", ""):
    raise ValueError("The trim_galore Snakemake wrapper cannot "
                       "be run with the `--fastqc` flag. Please "
                       "remove the flag from extra params. "
                       "You can use the fastqc Snakemake wrapper on "
                       "the input and output files instead.")

# Check that four output files were supplied
m = len(snakemake.output)
assert m == 4, "Output must contain 4 files. Given: %r." % m

# Check that all output files are in the same directory
out_dir = os.path.dirname(snakemake.output[0])
for file_path in snakemake.output[1:]:
    assert out_dir == os.path.dirname(file_path), \
        "trim_galore can only output files to a single directory." \
        " Please indicate only one directory for the output files."

shell(
    "(trim_galore"
    " {snakemake.params.extra}"
    " --paired"
    " -o {out_dir}"
    " {snakemake.input})"
    " {log}")
Python Snakemake FastQC Trim_Galore From line 3 of pe/wrapper.py 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
__author__ = "Julian de Ruiter"
__copyright__ = "Copyright 2017, Julian de Ruiter"
__email__ = "[email protected]"
__license__ = "MIT"


from os import path

from snakemake.shell import shell


input_dirs = set(path.dirname(fp) for fp in snakemake.input)
output_dir = path.dirname(snakemake.output[0])
output_name = path.basename(snakemake.output[0])
log = snakemake.log_fmt_shell(stdout=True, stderr=True)

shell(
    "multiqc"
    " {snakemake.params}"
    " --force"
    " -o {output_dir}"
    " -n {output_name}"
    " {input_dirs}"
    " {log}"
)
Python Snakemake MultiQC From line 3 of multiqc/wrapper.py 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
__author__ = "Julian de Ruiter"
__copyright__ = "Copyright 2017, Julian de Ruiter"
__email__ = "[email protected]"
__license__ = "MIT"


from os import path

from snakemake.shell import shell


input_dirs = set(path.dirname(fp) for fp in snakemake.input)
output_dir = path.dirname(snakemake.output[0])
output_name = path.basename(snakemake.output[0])
log = snakemake.log_fmt_shell(stdout=True, stderr=True)

shell(
    "multiqc"
    " {snakemake.params}"
    " --force"
    " -o {output_dir}"
    " -n {output_name}"
    " {input_dirs}"
    " {log}"
)
Python Snakemake MultiQC From line 3 of multiqc/wrapper.py 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
__author__ = "Julian de Ruiter"
__copyright__ = "Copyright 2017, Julian de Ruiter"
__email__ = "[email protected]"
__license__ = "MIT"


from os import path

from snakemake.shell import shell


input_dirs = set(path.dirname(fp) for fp in snakemake.input)
output_dir = path.dirname(snakemake.output[0])
output_name = path.basename(snakemake.output[0])
log = snakemake.log_fmt_shell(stdout=True, stderr=True)

shell(
    "multiqc"
    " {snakemake.params}"
    " --force"
    " -o {output_dir}"
    " -n {output_name}"
    " {input_dirs}"
    " {log}"
)
Python Snakemake MultiQC From line 3 of multiqc/wrapper.py 
 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
__author__ = "Johannes Köster"
__copyright__ = "Copyright 2016, Johannes Köster"
__email__ = "[email protected]"
__license__ = "MIT"


import os
from snakemake.shell import shell

extra = snakemake.params.get("extra", "")
log = snakemake.log_fmt_shell(stdout=True, stderr=True)

fq1 = snakemake.input.get("fq1")
assert fq1 is not None, "input-> fq1 is a required input parameter"
fq1 = (
    [snakemake.input.fq1]
    if isinstance(snakemake.input.fq1, str)
    else snakemake.input.fq1
)
fq2 = snakemake.input.get("fq2")
if fq2:
    fq2 = (
        [snakemake.input.fq2]
        if isinstance(snakemake.input.fq2, str)
        else snakemake.input.fq2
    )
    assert len(fq1) == len(
        fq2
    ), "input-> equal number of files required for fq1 and fq2"
input_str_fq1 = ",".join(fq1)
input_str_fq2 = ",".join(fq2) if fq2 is not None else ""
input_str = " ".join([input_str_fq1, input_str_fq2])

if fq1[0].endswith(".gz"):
    readcmd = "--readFilesCommand zcat"
else:
    readcmd = ""

outprefix = os.path.dirname(snakemake.output[0]) + "/"

shell(
    "STAR "
    "{extra} "
    "--runThreadN {snakemake.threads} "
    "--genomeDir {snakemake.params.index} "
    "--readFilesIn {input_str} "
    "{readcmd} "
    "--outFileNamePrefix {outprefix} "
    "--outStd Log "
    "{log}"
)
Python Snakemake STAR From line 1 of align/wrapper.py 
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
__author__ = "Thibault Dayris"
__copyright__ = "Copyright 2019, Dayris Thibault"
__email__ = "[email protected]"
__license__ = "MIT"

from snakemake.shell import shell
from snakemake.utils import makedirs

log = snakemake.log_fmt_shell(stdout=True, stderr=True)

extra = snakemake.params.get("extra", "")
sjdb_overhang = snakemake.params.get("sjdbOverhang", "100")

gtf = snakemake.input.get("gtf")
if gtf is not None:
    gtf = "--sjdbGTFfile " + gtf
    sjdb_overhang = "--sjdbOverhang " + sjdb_overhang
else:
    gtf = sjdb_overhang = ""

makedirs(snakemake.output)

shell(
    "STAR "  # Tool
    "--runMode genomeGenerate "  # Indexation mode
    "{extra} "  # Optional parameters
    "--runThreadN {snakemake.threads} "  # Number of threads
    "--genomeDir {snakemake.output} "  # Path to output
    "--genomeFastaFiles {snakemake.input.fasta} "  # Path to fasta files
    "{sjdb_overhang} "  # Read-len - 1
    "{gtf} "  # Highly recommended GTF
    "{log}"  # Logging
)
Python Snakemake STAR From line 3 of index/wrapper.py
ShowHide 41 more snippets with no or duplicated tags.

Login to post a comment if you would like to share your experience with this workflow.

Do you know this workflow well? If so, you can request seller status , and start supporting this workflow.

Free

Created: 1yr ago
Updated: 1yr ago
Maitainers: public
URL: https://github.com/akcorut/SnakeRNASeq
Name: snakernaseq
Version: 1
Badge:
workflow icon

Insert copied code into your website to add a link to this workflow.

Downloaded: 0
Copyright: Public Domain
License: None
  • Future updates

Related Workflows

psychip_snakemake
Show Details View Workflow
ENCODE pipeline for histone marks developed for the psychENCODE project
public

public

psychip pipeline is an improved version of the ENCODE pipeline for histone marks developed for the psychENCODE project. The o...
raw sequence reads Alignment Sequence alignment report macs2 ucsc-bedclip bedGraphToBigWig BEDTools BWA Picard SAMtools Snakemake
Free
74
ncov_2
Show Details View Workflow
Near-real time tracking of SARS-CoV-2 in Connecticut
public

public

Repository containing scripts to perform near-real time tracking of SARS-CoV-2 in Connecticut using genomic data. This pipeli...
JSON nextclade Augur Biopython FOCUS Pandas Snakemake bs4 epiweeks geopy matplotlib numpy pycountry pycountry-convert uszipcode
Free
28
cellranger-snakemake-gke
Show Details View Workflow
snakemake workflow to run cellranger on a given bucket using gke.
public

public

A Snakemake workflow for running cellranger on a given bucket using Google Kubernetes Engine. The usage of this workflow ...
macs2 ucsc-bedclip bedGraphToBigWig BEDTools BWA Picard SAMtools Snakemake
Free
74
atlas
Show Details View Workflow
ATLAS - Three commands to start analyzing your metagenome data
public

public

Metagenome-atlas is a easy-to-use metagenomic pipeline based on snakemake. It handles all steps from QC, Assembly, Binning, t...
raw sequence reads Genome assembly Annotation track checkm2 gunc prodigal snakemake-wrapper-utils MEGAHIT Atlas BBMap Biopython BioRuby Bwa-mem2 cd-hit CheckM DAS Diamond eggNOG-mapper v2 MetaBAT 2 Minimap2 MMseqs MultiQC Pandas Picard pyfastx SAMtools SemiBin Snakemake SPAdes SqueezeMeta TADpole VAMB CONCOCT ete3 gtdbtk h5py networkx numpy plotly psutil utils metagenomics
Free
175
rna-seq-star-deseq2
Show Details View Workflow
RNA-seq workflow using STAR and DESeq2
public

public

This workflow performs a differential gene expression analysis with STAR and Deseq2. The usage of this workflow is described ...
raw sequence reads Gene expression profiling ID list snakemake-wrapper-utils biomaRt PCAtools BiocParallel BWA Cutadapt DESeq2 GFFutils MultiQC Pandas SAMtools Snakemake STAR cli stringr tidyverse RNA-Seq
Free
41
dna-seq-gatk-variant-calling
Show Details View Workflow
This Snakemake pipeline implements the GATK best-practices workflow
public

public

This Snakemake pipeline implements the GATK best-practices workflow for calling small germline variants. The usage of thi...
VCF raw sequence reads Variant calling genetic variants gatk rust-bio-tools snakemake-wrapper-utils tabix BCFtools BWA FastQC MultiQC Pandas Picard SAMtools Snakemake Trimmomatic Variant Effect Predictor (VEP) common matplotlib numpy seaborn DNA
Free
57