Skip to content
This repository has been archived by the owner on Oct 23, 2023. It is now read-only.

Commit

Permalink
Updating Ebseq rules
Browse files Browse the repository at this point in the history
  • Loading branch information
@skchronicles
skchronicles committed Jun 11, 2019
1 parent 16ecfca commit 012965f
Showing 1 changed file with 28 additions and 23 deletions.
51 changes: 28 additions & 23 deletions Rules/rnaseq.snakefile
100755 → 100644
View file
Original file line number Diff line number Diff line change
Expand Up @@ -36,10 +36,14 @@ contrastsList = createConstrasts(config['project']['contrasts']['rcontrasts'])
cpm_cutoff_list = config['project']['contrasts']['rcontrasts_cpm_cutoff']
mincounts = config['project']['contrasts']['rcontrasts_min_counts']

# Creating DEG contrasts with CPM and Minimum Sample thresholds
contrasts_w_cpm_cutoff_list = []
print("Creating Contrasts with the following information")
groups_i = 0
for i in zip(contrastsList,cpm_cutoff_list,mincounts):
groups_i += 1
contrasts_w_cpm_cutoff_list.append(str(i[0]) + "_" + str(i[1]) + "_" + str(i[2]))
print(contrasts_w_cpm_cutoff_list, contrastsList)
print("\t{}. Groups: {}\tCPM:{}\tMinSamples:{}".format(groups_i, i[0], i[1], i[2]))

se=""
pe=""
Expand Down Expand Up @@ -67,7 +71,7 @@ rseqc_dir="RSeQC"

debug = expand(join(workpath,"DEG_{con}_{cpm}_{minsamples}", "DESeq2_DEG_{con}_all_genes.txt"), zip, con=contrastsList, cpm=cpm_cutoff_list,minsamples=mincounts) + expand(join(workpath,"DEG_{con}_{cpm}_{minsamples}", "limma_DEG_{con}_all_genes.txt"), zip, con=contrastsList, cpm=cpm_cutoff_list,minsamples=mincounts) + expand(join(workpath,"DEG_{con}_{cpm}_{minsamples}", "edgeR_DEG_{con}_all_genes.txt"), zip, con=contrastsList, cpm=cpm_cutoff_list,minsamples=mincounts),

print(debug)
#print(debug)

if config['project']['DEG'] == "yes" and config['project']['TRIM'] == "yes":
rule all:
Expand All @@ -78,23 +82,23 @@ if config['project']['DEG'] == "yes" and config['project']['TRIM'] == "yes":
#Initialize DEG (Filter raw counts matrix by cpm and min sample thresholds)
expand(join(workpath,"DEG_{deg_dir}", "RawCountFile_RSEM_genes_filtered.txt"), deg_dir=contrasts_w_cpm_cutoff_list),
expand(join(workpath,"DEG_{deg_dir}", "RawCountFile_RSEM_genes.txt"), deg_dir=contrasts_w_cpm_cutoff_list),

#EBSeq
expand(join(workpath,"DEG_{deg_dir}", "EBSeq_isoform_completed.txt"), deg_dir=contrasts_w_cpm_cutoff_list),
expand(join(workpath,"DEG_{deg_dir}", "EBSeq_gene_completed.txt"), deg_dir=contrasts_w_cpm_cutoff_list),


#DEG (limma, DESeq2, edgeR)
expand(join(workpath,"DEG_{deg_dir}", "DESeq2_PCA.png"), deg_dir=contrasts_w_cpm_cutoff_list),
expand(join(workpath,"DEG_{deg_dir}", "edgeR_prcomp.png"), deg_dir=contrasts_w_cpm_cutoff_list),
expand(join(workpath,"DEG_{deg_dir}", "limma_MDS.png"), deg_dir=contrasts_w_cpm_cutoff_list),

#EBSeq (Genes and Isoforms)
expand(join(workpath,"DEG_{con}_{cpm}_{minsamples}", "EBSeq_isoform_completed.txt"), zip, con=contrastsList, cpm=cpm_cutoff_list,minsamples=mincounts),
expand(join(workpath,"DEG_{con}_{cpm}_{minsamples}", "EBSeq_gene_completed.txt"), zip, con=contrastsList, cpm=cpm_cutoff_list,minsamples=mincounts),

expand(join(workpath,"DEG_{con}_{cpm}_{minsamples}", "DESeq2_DEG_{con}_all_genes.txt"), zip, con=contrastsList, cpm=cpm_cutoff_list,minsamples=mincounts),
expand(join(workpath,"DEG_{con}_{cpm}_{minsamples}", "limma_DEG_{con}_all_genes.txt"), zip, con=contrastsList, cpm=cpm_cutoff_list,minsamples=mincounts),
expand(join(workpath,"DEG_{con}_{cpm}_{minsamples}", "edgeR_DEG_{con}_all_genes.txt"), zip, con=contrastsList, cpm=cpm_cutoff_list,minsamples=mincounts),

#PCA (limma, DESeq2, edgeR)
expand(join(workpath,"DEG_{deg_dir}","PcaReport.html"), deg_dir=contrasts_w_cpm_cutoff_list),

#Venn Diagram (Overlap of sig genes between limma, DESeq2, edgeR)
expand(join(workpath,"DEG_{deg_dir}","limma_edgeR_DESeq2_vennDiagram.png"), deg_dir=contrasts_w_cpm_cutoff_list),

Expand Down Expand Up @@ -128,17 +132,16 @@ Rscript {params.rscript} {params.outdir} {wildcards.group1} {wildcards.group2} {

rule EBSeq_isoform:
input:
samtab=join(workpath,star_dir,"sampletable.txt"),
samtab2s=expand(join(workpath,"DEG_{deg_dir}","sampletable.txt"),deg_dir=contrasts_w_cpm_cutoff_list),
igfiles=expand(join(workpath,"DEG_ALL","{name}.RSEM.isoforms.results"),name=samples),
sampletable=join(workpath,"DEG_{con}_{cpm}_{minsamples}","sampletable.txt"),
isoformCounts=expand(join(workpath,"DEG_ALL","{name}.RSEM.isoforms.results"),name=samples),
output:
join(workpath,"DEG_{deg_dir}","EBSeq_isoform_completed.txt"),
join(workpath,"DEG_{con}_{cpm}_{minsamples}","EBSeq_isoform_completed.txt"),
params:
rname='pl:EBSeq_isoform',
batch='--mem=128g --cpus-per-task=8 --time=10:00:00',
workpath=workpath,
outdir=join(workpath,"DEG_{deg_dir}"),
contrasts=" ".join(config['project']['contrasts']['rcontrasts']),
outdir=join(workpath,"DEG_{con}_{cpm}_{minsamples}"),
contrasts= lambda w: str(w.con).replace("-", " "),
rsemref=config['references'][pfamily]['RSEMREF'],
rsem=config['bin'][pfamily]['RSEM'],
annotate=config['references'][pfamily]['ANNOTATEISOFORMS'],
Expand All @@ -149,20 +152,20 @@ rule EBSeq_isoform:
cd {params.outdir}
module load {params.pythonver}
module load {params.rsemver}
python {params.script1} '{params.outdir}' '{input.igfiles}' '{input.samtab}' '{params.contrasts}' '{params.rsemref}' 'isoform' '{params.annotate}'
python {params.script1} '{params.outdir}' '{input.isoformCounts}' '{input.sampletable}' '{params.contrasts}' '{params.rsemref}' 'isoform' '{params.annotate}'
"""

rule EBSeq_gene:
input:
samtab=join(workpath,star_dir,"sampletable.txt"),
igfiles=expand(join(workpath,"DEG_ALL","{name}.RSEM.genes.results"), name=samples),
sampletable=join(workpath,"DEG_{con}_{cpm}_{minsamples}","sampletable.txt"),
geneCounts=expand(join(workpath,"DEG_ALL","{name}.RSEM.genes.results"), name=samples),
output:
join(workpath,"DEG_{deg_dir}","EBSeq_gene_completed.txt"),
join(workpath,"DEG_{con}_{cpm}_{minsamples}","EBSeq_gene_completed.txt"),
params:
rname='pl:EBSeq_gene',
batch='--mem=128g --cpus-per-task=8 --time=10:00:00',
outdir=join(workpath,"DEG_{deg_dir}"),
contrasts=" ".join(config['project']['contrasts']['rcontrasts']),
outdir=join(workpath,"DEG_{con}_{cpm}_{minsamples}"),
contrasts= lambda w: str(w.con).replace("-", " "),
rsemref=config['references'][pfamily]['RSEMREF'],
rsem=config['bin'][pfamily]['RSEM'],
annotate=config['references'][pfamily]['ANNOTATE'],
Expand All @@ -173,7 +176,7 @@ rule EBSeq_gene:
cd {params.outdir}
module load {params.pythonver}
module load {params.rsemver}
python {params.script1} '{params.outdir}' '{input.igfiles}' '{input.samtab}' '{params.contrasts}' '{params.rsemref}' 'gene' '{params.annotate}'
python {params.script1} '{params.outdir}' '{input.geneCounts}' '{input.sampletable}' '{params.contrasts}' '{params.rsemref}' 'gene' '{params.annotate}'
"""


Expand Down Expand Up @@ -298,7 +301,7 @@ rule vennDiagram:
limmafile=join(workpath, "DEG_{con}_{cpm}_{minsamples}","limma_DEG_{con}_all_genes.txt"),
deseqfile=join(workpath, "DEG_{con}_{cpm}_{minsamples}","edgeR_DEG_{con}_all_genes.txt"),
edgeRfile=join(workpath, "DEG_{con}_{cpm}_{minsamples}","DESeq2_DEG_{con}_all_genes.txt"),
output:
output:
join(workpath,"DEG_{con}_{cpm}_{minsamples}","limma_edgeR_DESeq2_vennDiagram.png")
params:
rname='pl:vennDiagram',
Expand All @@ -313,3 +316,5 @@ module load {params.rver}
Rscript {params.rscript} --limma '{input.limmafile}' --edgeR '{input.edgeRfile}' --DESeq2 '{input.deseqfile}'
"""



0 comments on commit 012965f

Please sign in to comment.