Venn Diagram showing overlap of DE genes

Results-template/Scripts/limma_edgeR_DESeq2_venn.R

@@ -0,0 +1,81 @@
+#Example Usage: Rscript limma_edgeR_DESeq2_venn.R --contrast 'group1-group2' --limma limma_DEG_{group1}-{group2}_all_genes.txt --edgeR edgeR_DEG_{group1}-{group2}_all_genes.txt --DESeq2 DESeq2_deg_{group1}_vs_{group2}.txt
+
+rm(list=ls())
+library(Vennerable)
+library(argparse)
+library(dplyr)
+
+filter_by_fc_fdr <-function(df,fdr_cutoff,upreg_fc_cutoff,downreg_fc_cutoff) {
+  df=filter(df, fdr!='NA' & fc!='NA')
+  keep1=df$fdr<=fdr_cutoff
+  keep2=df$fc>=upreg_fc_cutoff
+  keep3=df$fc<=downreg_fc_cutoff
+  keep <- (keep1) & (keep2 | keep3)
+  return(df[keep,])
+}
+
+
+parser <- ArgumentParser()
+
+parser$add_argument("-c", "--contrast", type="character", required=TRUE,
+                    help="contrast for differential expression")
+
+parser$add_argument("-l", "--limma", type="character", required=TRUE,
+                    help="reformatted limma output filename ")
+
+parser$add_argument("-e", "--edgeR", type="character", required=TRUE,
+                    help="reformatted edgeR output filename ")
+
+parser$add_argument("-D", "--DESeq2", type="character", required=TRUE,
+                    help="reformatted DESeq2 output filename ")
+
+parser$add_argument("-f", "--fdr_cutoff", type="double", default=0.05,
+                    help="FDR cutoff")
+
+parser$add_argument("-u", "--upreg_fc_cutoff", type="double", default=1.5,
+                    help="upregulated fold change cutoff")
+
+parser$add_argument("-d", "--downreg_fc_cutoff", type="double", default=-1.5,
+                    help="upregulated fold change cutoff")
+
+args <- parser$parse_args()
+
+limma=read.table(args$limma,header = TRUE)
+edger=read.table(args$edgeR,header=TRUE)
+deseq2=read.table(args$DESeq2,header=TRUE)
+
+fdr_cutoff=args$fdr_cutoff
+upreg_fc_cutoff=args$upreg_fc_cutoff
+downreg_fc_cutoff=args$downreg_fc_cutoff
+
+
+limma_filt=filter_by_fc_fdr(limma,fdr_cutoff,upreg_fc_cutoff,downreg_fc_cutoff)
+edger_filt=filter_by_fc_fdr(edger,fdr_cutoff,upreg_fc_cutoff,downreg_fc_cutoff)
+deseq2_filt=filter_by_fc_fdr(deseq2,fdr_cutoff,upreg_fc_cutoff,downreg_fc_cutoff)
+
+dim(limma_filt)
+limma_filt
+dim(edger_filt)
+edger_filt
+dim(deseq2_filt)
+deseq2_filt
+
+x.genes=limma_filt$ensid_gene
+y.genes=edger_filt$ensid_gene
+z.genes=deseq2_filt$ensid_gene
+
+vennD=Venn(SetNames = c("limma","DESeq2","edgeR"), 
+           Weight=c(`100`=length(setdiff(x.genes,union(y.genes,z.genes))),
+                    `010`=length(setdiff(y.genes,union(x.genes,z.genes))),
+                    `110`=length(setdiff(intersect(x.genes,y.genes),z.genes)),
+                    `001`=length(setdiff(z.genes,union(x.genes,y.genes))),
+                    `101`=length(setdiff(intersect(x.genes,z.genes),y.genes)),
+                    `011`=length(setdiff(intersect(z.genes,y.genes),x.genes)),
+                    `111`=length(intersect(intersect(x.genes,y.genes),z.genes))))
+
+outfile = paste("limma_edgeR_DESeq2_",args$contrast,"_vennDiagram.png",sep="")
+
+png(outfile)
+plot(vennD, doWeights = FALSE, type = "circles")
+dev.off()
+

Rules/rnaseq.snakefile

@@ -17,12 +17,21 @@ def check_writeaccess(filename):
   if not os.access(filename,os.W_OK):
     exit("File: %s exists, but cannot be read!"%(filename))
 
+def createConstrasts(cList):
+  contrastsList = []
+  for i in range(0, len(cList)-1, 2):
+    contrastsList.append("-".join(cList[i:i+2]))
+  return contrastsList
+
+
 
 configfile: "run.json"
 
 samples=config['project']['groups']['rsamps']
 
-workpath = config['project']['workpath']
+workpath=config['project']['workpath']
+
+contrastsList = createConstrasts(config['project']['contrasts']['rcontrasts'])
 
 se=""
 pe=""
@@ -74,6 +83,7 @@ if config['project']['DEG'] == "yes" and config['project']['TRIM'] == "yes":
       join(workpath,subreadg_dir,"edgeR_prcomp.png"),
       join(workpath,subreadg_dir,"limma_MDS.png"),
       join(workpath,subreadg_dir,"PcaReport.html"),
+      expand(join(workpath,subreadg_dir,"limma_edgeR_DESeq2_{con}_vennDiagram.png"),con=contrastsList),
 
       #Subread-junctions
 
@@ -82,6 +92,7 @@ if config['project']['DEG'] == "yes" and config['project']['TRIM'] == "yes":
       join(workpath,subreadj_dir,"edgeR_prcomp.png"),
       join(workpath,subreadj_dir,"limma_MDS.png"),
       join(workpath,subreadj_dir,"PcaReport.html"),
+      expand(join(workpath,subreadj_dir,"limma_edgeR_DESeq2_{con}_vennDiagram.png"),con=contrastsList),
 
       #Subread-gene-junctions
 
@@ -90,6 +101,7 @@ if config['project']['DEG'] == "yes" and config['project']['TRIM'] == "yes":
       join(workpath,subreadgj_dir,"edgeR_prcomp.png"),
       join(workpath,subreadgj_dir,"limma_MDS.png"),
       join(workpath,subreadgj_dir,"PcaReport.html"),
+      expand(join(workpath,subreadgj_dir,"limma_edgeR_DESeq2_{con}_vennDiagram.png"),con=contrastsList),
 
       #RSEM
 
@@ -98,13 +110,17 @@ if config['project']['DEG'] == "yes" and config['project']['TRIM'] == "yes":
       join(workpath,rsemg_dir,"edgeR_prcomp.png"),
       join(workpath,rsemg_dir,"limma_MDS.png"),
       join(workpath,rsemg_dir,"PcaReport.html"),
+      expand(join(workpath,rsemg_dir,"limma_edgeR_DESeq2_{con}_vennDiagram.png"),con=contrastsList),
       expand(join(workpath,rsemg_dir,"{name}.RSEM.genes.results"),name=samples),
       join(workpath,rsemg_dir,"RSEM.genes.FPKM.all_samples.txt"),
       join(workpath,rsemi_dir,"RSEM.isoforms.FPKM.all_samples.txt"),
 
       expand(join(workpath,star_dir,"{name}.star.count.overlap.txt"),name=samples),
       join(workpath,star_dir,"RawCountFileOverlap.txt"),
       join(workpath,star_dir,"RawCountFileStar.txt"),
+
+
+
 
 
 elif config['project']['DEG'] == "no" and config['project']['TRIM'] == "yes":
@@ -490,6 +506,30 @@ module load {params.rver}
 Rscript limmacall.R '{params.outdir}' '{input.file1}' '{input.file2}' '{params.contrasts}' '{params.refer}' '{params.projectId}' '{params.projDesc}' '{params.gtffile}' '{params.dtype}' '{params.karyobeds}'
 """
 
+rule vennDiagram:
+  input:
+    join(workpath,"DEG_{dtype}","limma_MDS.png"),
+    join(workpath,"DEG_{dtype}","edgeR_prcomp.png"),
+    join(workpath,"DEG_{dtype}","DESeq2_PCA.png"),
+    limmafile=join(workpath,"DEG_{dtype}","limma_DEG_{group1}-{group2}_all_genes.txt"),
+    deseq2file=join(workpath,"DEG_{dtype}","edgeR_DEG_{group1}-{group2}_all_genes.txt"),
+    edgeRfile=join(workpath,"DEG_{dtype}","DESeq2_deg_{group1}_vs_{group2}.txt"),
+  output: 
+    join(workpath,"DEG_{dtype}","limma_edgeR_DESeq2_{group1}-{group2}_vennDiagram.png")
+  params:
+    rname='pl:vennDiagram',
+    batch='--mem=12g --time=2:00:00',
+    outdir=join(workpath,"DEG_{dtype}"),
+    dtype="{dtype}",
+    rver=config['bin'][pfamily]['tool_versions']['RVER'],
+    rscript=join(workpath,"Scripts","limma_edgeR_DESeq2_venn.R"),
+  shell: """
+cp {params.rscript} {params.outdir}
+cd {params.outdir}
+module load {params.rver}
+Rscript {params.rscript} --contrast '{wildcards.group1}-{wildcards.group2}' --limma '{input.limmafile}' --edgeR '{input.edgeRfile}' --DESeq2 '{input.deseq2file}'
+"""
+
 rule pca:
   input: 
     file1=join(workpath,star_dir,"sampletable.txt"),
@@ -575,3 +615,4 @@ module load {params.pythonver}
 module load {params.rsemver}
 python {params.script1} '{params.outdir}' '{input.igfiles}' '{input.samtab}' '{params.contrasts}' '{params.rsemref}' 'gene' '{params.annotate}'
 """
+