rMATs and IsoformSwitchAnalyeR workflow for multi-group multi-contrasts scenarios

public 1yr ago Version: v0.1.1 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 .

ASE (alternate splicing events) are identified and quantified using the ASCENT ( A lternate S pli C ing E ve N t T ools) pipeline. This workflow can be used for multi-group multi-contrasts scenarios.

ASCENT currently runs:

NOTE : More ASE tools will be added to ASCENT in the future as and when needed. Please contact Vishal Koparde to make requests to add specific ASE tool(s).

Usage:

 % bash /data/Ziegelbauer_lab/Pipelines/CCBR_ASE/v0.1.1/run
Pipeline Dir: /gpfs/gsfs12/users/Ziegelbauer_lab/Pipelines/CCBR_ASE/v0.1.1
Git Commit/Tag: a41b33007b6063b1d95744808f527a126fb1f0dc	v0.1.1
/gpfs/gsfs12/users/Ziegelbauer_lab/Pipelines/CCBR_ASE/v0.1.1/run
--> run CCBR Alternate Splicing Events Pipeline
USAGE:
 bash /data/Ziegelbauer_lab/Pipelines/CCBR_ASE/v0.1.1/run -m/--runmode=<RUNMODE> -w/--workdir=<WORKDIR>
Required Arguments:
1. RUNMODE: [Type: String] Valid options:
 *) init : initialize workdir
 *) run : run with slurm
 *) reset : DELETE workdir dir and re-init it
 *) dryrun : dry run snakemake to generate DAG
 *) unlock : unlock workdir if locked by snakemake
 *) runlocal : run without submitting to sbatch
2. WORKDIR: [Type: String]: Absolute or relative path to the output folder with write permissions.

DISCLAIMER :

IsoformSwitchAnalyzeR fails if replicates are absent.

References:

Code Snippets

    shell:"""
set -exuf -o pipefail
if [ ! -d {params.parentdir} ];then
    mkdir -p {params.parentdir}
fi
for i in {input}; do
    bn=$(basename $i)
    sn=$(echo $bn|awk -F".RSEM" '{{print $1}}')
    mkdir -p {params.parentdir}/${{sn}}
    ln $i {params.parentdir}/${{sn}}/${{bn}}
done
touch {output.dummy}
"""

SnakeMake From line 22 of rules/isa.smk

    shell:"""
set -exuf -o pipefail
c1=$(echo "{params.contrast}"|awk -F"_vs_" '{{print $1}}')
c2=$(echo "{params.contrast}"|awk -F"_vs_" '{{print $2}}')
Rscript {params.rscript} \
 -p {params.parentdir} \
 -g {params.gtf} \
 -s {params.samplestsv} \
 -c $c1 \
 -d $c2 \
 -o {params.outdir}
"""

SnakeMake From line 50 of rules/isa.smk

    shell:"""
set -exuf -o pipefail
Rscript {params.rscript} -d {params.resultsdir}
"""

SnakeMake From line 74 of rules/isa.smk

    shell:"""
set -euf -o pipefail
outdir=$(dirname {output.summary})
cp "{input.b1}" $outdir/
cp "{input.b2}" $outdir/
python ${{RMATS_SRC}}/rmats.py \
    --b1 "{input.b1}" \
    --b2 "{input.b2}" \
    --od "$outdir" \
    --gtf "{params.gtf}" \
    --readLength {params.rl} \
    --bi "{params.starindexdir}" \
    --nthread {threads} \
    -t "paired" \
    --novelSS \
    --libType "fr-secondstrand" \
    --tmp /lscratch/${{SLURM_JOB_ID}}/
"""

SnakeMake rmats From line 30 of rules/rmats.smk

    shell:"""
set -exuf -o pipefail
Rscript {params.rscript} -d {params.resultsdir}
"""

SnakeMake From line 60 of rules/rmats.smk

    shell:"""
set -euf -o pipefail
cd {params.rsemindexdir}
rsem-prepare-reference -p {threads} --gtf {input.gtf} {input.reffa} ref
rsem-generate-ngvector ref.transcripts.fa ref.transcripts
"""

SnakeMake RSEM From line 11 of rules/rsem.smk

    shell:"""
set -euf -o pipefail
cd {params.rsemindexdir}
genesgtf={input.gtf}
bn=$(basename $genesgtf)
gtfToGenePred -genePredExt -geneNameAsName2 $genesgtf ${{bn}}.tmp
awk -v OFS="\\t" '{{print $2,$4,$5,$1,"0",$3,$6,$7,"0",$8,$9,$10}}' ${{bn}}.tmp > {output.bed12}
rm -f ${{bn}}.tmp
"""

SnakeMake gtftogenepred From line 27 of rules/rsem.smk

    shell:"""
set -euf -o pipefail
mkdir -p {params.stardir}
STAR \
    --runThreadN {threads} \
    --runMode genomeGenerate \
    --genomeDir {params.stardir} \
    --genomeFastaFiles {input.reffa} \
    --outTmpDir /lscratch/${{SLURM_JOB_ID}}/tmp_{params.rl}
"""

SnakeMake STAR From line 13 of rules/star.smk

    shell:"""
set -euf -o pipefail
overhang=$(echo {params.rl}|awk '{{print $1-1}}')
cd {params.bamdir}
STAR \
    --chimSegmentMin 2 \
    --outFilterMismatchNmax 3 \
    --alignEndsType EndToEnd \
    --runThreadN {threads} \
    --outSAMstrandField intronMotif \
    --outSAMtype BAM SortedByCoordinate \
    --alignSJDBoverhangMin 6 \
    --alignIntronMax 299999 \
    --genomeDir {params.starindexdir} \
    --sjdbGTFfile {params.gtf} \
    --outFileNamePrefix "{params.sample}." \
    --readFilesIn {input.R1} {input.R2} \
    --readFilesCommand zcat \
    --outTmpDir=/lscratch/${{SLURM_JOB_ID}}/{params.sample} \
    --sjdbOverhang $overhang \
    --quantMode GeneCounts TranscriptomeSAM
"""

SnakeMake STAR From line 43 of rules/star.smk

    shell:"""
set -euf -o pipefail
cd {params.rsemdir}
# inter strandedness
samtools index -@{threads} {input.bam}
infer_experiment.py -r {input.bed12} -i {input.bam} -s 1000000 > {output.strandinfo}
# Get strandedness to calculate Forward Probability
fp=$(tail -n1 {output.strandinfo} | 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"
rsemindex=$(echo {input.bed12}|sed "s@.bed12@@g")
rsem-calculate-expression --no-bam-output --calc-ci --seed 12345  \
        --bam --paired-end -p {threads}  {input.tbam} $rsemindex {params.sample} --time \
        --temporary-folder /lscratch/$SLURM_JOBID --keep-intermediate-files --forward-prob=${{fp}} --estimate-rspd
mv {params.sample}.genes.results {output.gcounts}
mv {params.sample}.isoforms.results {output.tcounts}
"""    