adding strand specificity to rsem counts command; adding rule to conv…

…ert bams to forward and reverse bigwigs for easy visualization

Rules/initialqcrnaseq.snakefile

@@ -811,6 +811,7 @@ if pe=="yes":
    rule rsem:
       input: 
         file1=join(workpath,bams_dir,"{name}.p2.Aligned.toTranscriptome.out.bam"),
+        file2=join(workpath,rseqc_dir,"{name}.strand.info")
       output: 
         out1=join(workpath,degall_dir,"{name}.RSEM.genes.results"),
         out2=join(workpath,degall_dir,"{name}.RSEM.isoforms.results"),
@@ -829,9 +830,12 @@ if pe=="yes":
 if [ ! -d {params.outdir} ]; then mkdir {params.outdir}; fi
 cd {params.outdir}
 module load {params.rsemver}
-rsem-calculate-expression --no-bam-output --calc-ci --seed 12345  --bam --paired-end -p {threads}  {input.file1} {params.rsemref} {params.prefix} --time --temporary-folder /lscratch/$SLURM_JOBID --keep-intermediate-files
+fp=`awk '{{if($NF > 0.75) print "0.0"; else if ($NF<0.25) print "D1.0"; else print "0.5";}}' {input.file2}`
+echo $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 /lscratch/$SLURM_JOBID --keep-intermediate-files --forward-prob=$fp --estimate-rpsd
 """
 
+
 if se=="yes":
 
    rule rsem:
@@ -855,7 +859,40 @@ if se=="yes":
 if [ ! -d {params.outdir} ]; then mkdir {params.outdir}; fi
 cd {params.outdir}
 module load {params.rsemver}
-rsem-calculate-expression --no-bam-output --calc-ci --seed 12345  --bam -p {threads}  {input.file1} {params.rsemref} {params.prefix} --time --temporary-folder /lscratch/$SLURM_JOBID --keep-intermediate-files
+fp=`awk '{{if($NF > 0.75) print "0.0"; else if ($NF<0.25) print "D1.0"; else print "0.5";}}' {input.file2}`
+echo $fp
+rsem-calculate-expression --no-bam-output --calc-ci --seed 12345  --bam -p {threads}  {input.file1} {params.rsemref} {params.prefix} --time --temporary-folder /lscratch/$SLURM_JOBID --keep-intermediate-files --forward-prob=$fp --estimate-rpsd
+"""
+
+rule bam2bw:
+	input:
+		bam=join(workpath,bams_dir,"{name}.star_rg_added.sorted.dmark.bam"),
+		strandinfo=join(workpath,rseqc_dir,"{name}.strand.info")
+	output:
+		fbw=join(workpath,bams_dir,"{name}.fwd.bw"),
+		rbw=join(workpath,bams_dir,"{name}.rev.bw")
+	params:
+		rname='pl:bam2bw',
+		prefix="{name}",
+		deeptoolsver=config['bin'][pfamily]['tool_versions']['DEEPTOOLSVER']
+	threads: 32
+	shell:"""
+module load {params.deeptoolsver};
+
+bam={input.bam}
+# Forward strand
+bamCoverage -b $bam -o {params.prefix}.fwd.bw --filterRNAstrand forward --binSize 20 --smoothLength 40 -p 32
+
+# Reverse strand
+bamCoverage -b $bam -o {params.prefix}.rev.bw --filterRNAstrand reverse --binSize 20 --smoothLength 40 -p 32
+
+# reverse files if method is not dUTP/NSR/NNSR ... ie, R1 in the direction of RNA strand.
+fp=`awk '{{if($NF > 0.75) print "0.0"; else if ($NF<0.25) print "D1.0"; else print "0.5";}}' {input.strandinfo}`
+if [ $fp -lt 0.25 ];then
+mv {params.prefix}.fwd.bw {params.prefix}.fwd.bw.tmp
+mv {params.prefix}.rev.bw {params.prefix}.fwd.bw
+mv {params.prefix}.fwd.bw.tmp {params.prefix}.rev.bw
+fi
 """
 
 
@@ -929,7 +966,7 @@ rule rseqc:
    shell: """
 module load {params.rseqcver}
 cd {rseqc_dir}
-infer_experiment.py -r {params.bedref} -i {input.file1} > {output.out1}
+infer_experiment.py -r {params.bedref} -i {input.file1} -s 1000000 > {output.out1}
 read_distribution.py -i {input.file1} -r {params.bedref} > {output.out4}
 """
 

cluster.json

@@ -2,7 +2,7 @@
     "__default__": {
         "gres": "lscratch:256",
         "mem": "48g",
-        "partition": "norm",
+        "partition": "ccr,norm",
         "threads": "4",
         "time": "4-00:00:00"
     },

standard-bin.json

@@ -95,6 +95,7 @@
         "RSEQCVER" : "rseqc/2.6.4",
         "RVER" : "R/3.5",
         "SAMTOOLSVER" : "samtools/1.6",
+        "DEEPTOOLSVER": "deeptools/3.0.1",
         "STARVER": "STAR/2.5.2b",
         "SUBREADVER": "subread/1.5.2",
         "TRIMMOMATICVER": "trimmomatic/0.33",