RNA-Seq Data Analysis Workflow using Snakemake and Salmon for Gene Expression Quantification and Differential Expression Analysis

public 1yr ago 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

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 .

A simple Snakemake workflow for processing RNA-Seq using snakemake in conjunction with salmon . Data is assumed to be paired-end , and a two-treatment design is assumed, although multiple paired comparisons are also possible The steps implemented are:

Prepare the reference
- Download the fa and gtf files for the relevant build
- Prepare the set of decoy transcripts
- Index the reference using the set of decoy transcripts
Pre-process raw fq files
- Run QC and prepare a report on the raw fq files ( FastQC / ngsReports )
- Identify adapters ( bbtools )
- Trim samples ( AdapterRemoval )
Pre-process raw fq files trimmed samples
- Run QC and prepare a report on the raw fq files ( FastQC / ngsReports )
Align and quantify ( salmon quant )
Gene Level Differential Expression Analysis
- Normalise using cqn
- Compare groups using quasi-likelihood approaches and glmTreat
- Perform an alternative analysis using limma-voom

Required Setup

Prior to executing the workflow, please ensure the following steps have been performed

Place unprocessed fastq (or fastq.gz ) files in data/raw/fastq . This path is hard-wired into the workflow and cannot be changed. If samples have been run across multiple lanes, please merge prior to running this workflow as all samples are expected to be in a single pair of files.
Prepare a tab-separated file, usually samples.tsv and place this in the config directory. The column name sample is mandatory, however any other columns required in the analysis can be specified in config.yml . When entering filenames in the sample column, you can leave off any suffixes (e.g. fastq.gz) as these are specified in config.yml . File names are expected to finish in _R[12].fastq.gz and only the section preceding this should be used for the samples column
Edit config.yml as required. Fields cannot be changed

Snakemake implementation

The basic workflow is written snakemake and can be called using the following steps.

Firstly, set-up the required conda environments

snakemake \	--use-conda \	--conda-prefix '/home/steveped/mambaforge/envs/' \	--conda-create-envs-only \	--cores 1

Secondly, create and inspect the rulegraph

snakemake --rulegraph > workflow/rules/rulegraph.dot
dot -Tpdf workflow/rules/rulegraph.dot > workflow/rules/rulegraph.pdf

Finally, the workflow itself can be run using:

snakemake \	-p \	--use-conda \	--conda-prefix '/home/steveped/mambaforge/envs/' \	--notemp \	--keep-going \	--cores 16

All data (salmon quants, trimmed fastq etc) will be written to the data directory. Trimmed fastq and salmon quant aux_info files are marked as temporary and can be safely deleted after completion of the analysis
RMarkdown files will be added to the anaylsis directory using modules provided in the workflow
Compiled html files will be written to the docs folder, along with any topTable files

Note that this creates common environments able to be called by other workflows and is dependent on the user. For me, my global conda environments are stored in /home/steveped/mambaforge/envs/ . For other users, this path will need to be modified.

Code Snippets

shell:
    """
    Rscript --vanilla {input} &>> {log}
    if [[ {params.git} == "True" ]]; then
            TRIES={params.tries}
            while [[ -f .git/index.lock ]]
            do
                if [[ "$TRIES" == 0 ]]; then
                    echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
                    exit 1
                fi
                sleep {params.interval}
                ((TRIES--))
            done
            git add {output}
    fi
    """

SnakeMake From line 11 of rules/build_site.smk

shell:
    """
    ## Make sure existing files are not overwritten,
    ## unless explicitly requested
    if [[ ! -f {output.rmd} || {params.overwrite} == 'True' ]]; then
        cp {input.rmd} {output.rmd}
    fi

    R -e "rmarkdown::render_site('{output.rmd}')" &>> {log}
    if [[ {params.git} == "True" ]]; then
            TRIES={params.tries}
            while [[ -f .git/index.lock ]]
            do
                if [[ "$TRIES" == 0 ]]; then
                    echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
                    exit 1
                fi
                sleep {params.interval}
                ((TRIES--))
            done
            git add {output}
    fi
    """

SnakeMake From line 50 of rules/build_site.smk

shell:
    """
    ## Make sure existing files are not overwritten,
    ## unless explicitly requested
    if [[ ! -f {output.rmd} || {params.overwrite} == 'True' ]]; then
        cp {input.rmd} {output.rmd}
    fi

    R -e "rmarkdown::render_site('{output.rmd}')" &>> {log}
    if [[ {params.git} == "True" ]]; then
            TRIES={params.tries}
            while [[ -f .git/index.lock ]]
            do
                if [[ "$TRIES" == 0 ]]; then
                    echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
                    exit 1
                fi
                sleep {params.interval}
                ((TRIES--))
            done
            git add {output}
    fi
    """

SnakeMake From line 96 of rules/build_site.smk

shell:
    """
    ## Make sure existing files are not overwritten,
    ## unless explicitly requested
    if [[ ! -f {output.rmd} || {params.overwrite} == 'True' ]]; then
        cp {input.rmd} {output.rmd}
    fi

    R -e "rmarkdown::render_site('{output.rmd}')" &>> {log}
    if [[ {params.git} == "True" ]]; then
            TRIES={params.tries}
            while [[ -f .git/index.lock ]]
            do
                if [[ "$TRIES" == 0 ]]; then
                    echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
                    exit 1
                fi
                sleep {params.interval}
                ((TRIES--))
            done
            git add {output}
        fi
    """

SnakeMake From line 143 of rules/build_site.smk

shell:
    """
    ## Make sure existing files are not overwritten,
    ## unless explicitly requested
    if [[ ! -f {output.rmd} || {params.overwrite} == 'True' ]]; then
        cp {input.rmd} {output.rmd}
    fi

    R -e "rmarkdown::render_site('{output.rmd}')" &>> {log}
    if [[ {params.git} == "True" ]]; then
        TRIES={params.tries}
        while [[ -f .git/index.lock ]]
        do
            if [[ "$TRIES" == 0 ]]; then
                echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
                exit 1
            fi
            sleep {params.interval}
            ((TRIES--))
        done
        git add {output}
    fi
    """

SnakeMake From line 185 of rules/build_site.smk

shell:
    """
    ## Make sure existing files are not overwritten,
    ## unless explicitly requested
    if [[ ! -f {output.rmd} || {params.overwrite} == 'True' ]]; then
        Rscript --vanilla \
            {input.script} \
            {params.comp} \
            {output.rmd} &>>{log}
    fi

    if [[ {params.git} == "True" ]]; then
            TRIES={params.tries}
            while [[ -f .git/index.lock ]]
            do
                if [[ "$TRIES" == 0 ]]; then
                    echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
                    exit 1
                fi
                sleep {params.interval}
                ((TRIES--))
            done
            git add {output}
    fi
    """

SnakeMake From line 17 of rules/dge.smk

shell:
    """
    R -e "rmarkdown::render_site('{input.rmd}')" &>> {log}
    if [[ {params.git} == "True" ]]; then
            TRIES={params.tries}
            while [[ -f .git/index.lock ]]
            do
                if [[ "$TRIES" == 0 ]]; then
                    echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
                    exit 1
                fi
                sleep {params.interval}
                ((TRIES--))
            done
            git add {output}
    fi
    """

SnakeMake From line 65 of rules/dge.smk

shell:
    """
    # Write to a separate temp directory for each run to avoid I/O clashes
    TEMPDIR=$(mktemp -d -t fqcXXXXXXXXXX)
    fastqc \
      {params.extra} \
      -t {threads} \
      --outdir $TEMPDIR \
      {input.fq} &> {log}
    # Move the files
    mv $TEMPDIR/*html $(dirname {output.html})
    mv $TEMPDIR/*zip $(dirname {output.zip})
    # Clean up the temp directory
    rm -rf $TEMPDIR

    ## Update the git repo if one is active
    if [[ {params.git} == "True" ]]; then
        TRIES={params.tries}
        while [[ -f .git/index.lock ]]
        do
            if [[ "$TRIES" == 0 ]]; then
                echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
                exit 1
            fi
            sleep {params.interval}
            ((TRIES--))
        done
        git add {output}
    fi
    """

SnakeMake FastQC From line 18 of rules/process_raw.smk

    shell:
        """
        bbmerge.sh \
            threads={threads} \
            in1={input.r1} \
            in2={input.r2} \
            outa={output.fa} 2> {log}

        if [[ {params.git} == "True" ]]; then
			TRIES={params.tries}
			while [[ -f .git/index.lock ]]
			do
				if [[ "$TRIES" == 0 ]]; then
					echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
					exit 1
				fi
				sleep {params.interval}
				((TRIES--))
			done
			git add {output}
		fi
        """

SnakeMake From line 62 of rules/process_raw.smk

    shell:
        """
        AdapterRemoval \
            --adapter1 {params.adapter1} \
            --adapter2 {params.adapter2} \
            --file1 {input.r1} \
            --file2 {input.r2} \
            {params.extra} \
            --threads {threads} \
            --maxns {params.maxns} \
            --minquality {params.minqual} \
            --minlength {params.minlength} \
            --output1 {output.r1} \
            --output2 {output.r2} \
            --discarded /dev/null \
            --singleton /dev/null \
            --settings {output.log} &> {log}

        if [[ {params.git} == "True" ]]; then
			TRIES={params.tries}
			while [[ -f .git/index.lock ]]
			do
				if [[ "$TRIES" == 0 ]]; then
					echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
					exit 1
				fi
				sleep {params.interval}
				((TRIES--))
			done
			git add {output.log}
		fi
        """

SnakeMake AdapterRemoval From line 112 of rules/process_raw.smk

shell:
	"""
	## Download and compress to bgzip on the fly
	curl {params.url} 2> {log} | zcat | bgzip -c > {output}
	"""

SnakeMake From line 44 of rules/reference.smk

shell:
	"""
	samtools faidx \
	  --gzi-idx {output.gzi} \
	  --fai-idx {output.fai} \
	  {input}
	"""

SnakeMake SAMtools From line 58 of rules/reference.smk

shell:
	"""
	samtools faidx {input.fa} \
		1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 MT X Y | \
		gzip > {output}
	"""

SnakeMake SAMtools From line 76 of rules/reference.smk

shell:
	"""
	curl {params.url} 2> {log} | zcat | bgzip -c > {output}
	"""

SnakeMake From line 90 of rules/reference.smk

shell:
	"""
	samtools faidx \
	  --gzi-idx {output.gzi} \
	  --fai-idx {output.fai} \
	  {input}
	"""

SnakeMake SAMtools From line 103 of rules/reference.smk

shell:
	"""
	zcat {input.fa} | \
		egrep "{params.regexp}" | \
		cut -f1 -d\ | \
		sed -r 's/^>//g' > {output.ids}
	samtools faidx \
	  -r {output.ids} \
	  {input.fa} | \
	  bgzip -c > {output.fa}
	"""

SnakeMake SAMtools From line 124 of rules/reference.smk

shell:
	"""
	curl -o {output} {params.url} 2> {log}
	"""

SnakeMake From line 143 of rules/reference.smk

shell:
	"""
	grep "^>" <(gunzip -c {input}) | cut -d " " -f 1 > {output}
	sed -i.bak -e 's/>//g' {output}
	"""

SnakeMake From line 152 of rules/reference.smk

shell:
	"""
	cat {input.trans_fa} {input.ref_fa} > {output}
	"""

SnakeMake From line 164 of rules/reference.smk

shell:
	"""
	salmon index \
		-t {input.gentrome} \
		-d {input.decoys} \
		-p {threads} \
		-i {output} 2> {log}
	"""

SnakeMake Salmon From line 177 of rules/reference.smk

shell:
    """
    salmon quant \
        -i {input.index} \
        -l {params.lib} \
        -p {threads} \
        --numBootstraps {params.nBoot} \
        --numGibbsSamples {params.nGibbs} \
        --thinningFactor {params.thin} \
        {params.extra} \
        -o {params.dir} \
        -1 {input.r1} \
        -2 {input.r2} &>> {log}

    if [[ {params.git} == "True" ]]; then
        TRIES={params.tries}
        while [[ -f .git/index.lock ]]
        do
            if [[ "$TRIES" == 0 ]]; then
                echo "ERROR: Timeout while waiting for removal of git index.lock" &>> {log}
                exit 1
            fi
            sleep {params.interval}
            ((TRIES--))
        done
        git add {output.quant} {output.json} {output.meta}
    fi
    """