add the deseq2 clustering heatmap

Author: Lin Yang

modules/differential_expression/differential_expression_cohort.snakefile

@@ -57,6 +57,7 @@ rule deseq2_differential_genes:
         treatment = treatment,
         control = control,
         meta = config["metasheet"],
+        pgenes = "static/deseq2/hg_pcoding.csv",
         path = "set +eu;source activate %s" % config['stat_root']
     message:
         "Running DESeq2 on the samples"
@@ -66,7 +67,7 @@ rule deseq2_differential_genes:
     shell:
     	"{params.path}; Rscript src/differentialexpr/DESeq2.R --input {params.filelist} --type salmon \
         --batch {params.batch} --meta {params.meta} --tx2gene {params.tx_annot} \
-        --condition {params.condition} --treatment {params.treatment} --control {params.control} --outpath {params.out_path}"
+        --condition {params.condition} --pcoding {params.pgenes} --treatment {params.treatment} --control {params.control} --outpath {params.out_path}"
 
 
 

src/differentialexpr/DESeq2.R

@@ -24,7 +24,9 @@ option_list = list(
   make_option(c("-r", "--treatment"), type="character", default="./",
               help="Treatment", metavar="character"),
   make_option(c("-c", "--control"), type="character", default="./",
-              help="Control", metavar="character")
+              help="Control", metavar="character"),
+  make_option(c("-p", "--pcoding"), type="character", default="./",
+              help="proding coding gene list", metavar="character")
 )
 
 
@@ -130,9 +132,10 @@ print(class(dds))
 print (paste("Comparing ",opt$treatment , " VS ", opt$control, sep = ""))
 res <- results(dds, contrast = c("Condition",c(opt$treatment,opt$control)))
 
+clustering_heatmap(dds, res, opt$pcoding, opt$outpath, opt$treatment, opt$control)
 
 res_final <- as.data.frame(res)
 res_final$Gene_name <- rownames(res_final)
 res_final <- res_final[c(7, 1:6)]
 res_final$`-log10(padj)` <- -log10(res_final$padj)
-write.table(res_final,file = paste(opt$outpath,opt$condition,'_',opt$treatment,'_vs_',opt$control,'_DESeq2.txt',sep = ""), quote = FALSE,sep = "\t")
+write.table(res_final,file = paste(opt$outpath,opt$condition,'_',opt$treatment,'_vs_',opt$control,'_DESeq2.txt',sep = ""), quote = FALSE,sep = "\t")

src/differentialexpr/clusteringheatmap.R

@@ -0,0 +1,108 @@
+library(vegan)
+library(pheatmap)
+
+clustering_heatmap <- function(dds, res, pgenes, outpath, treatment, control) {
+  
+  print("loding the human protein coding genes...")
+  pgenes <- read.csv(pgenes)	
+  deseq_result <- res
+  deseq_result <- as.data.frame(deseq_result)
+  
+  deseq_result <- deseq_result[rownames(deseq_result) %in% pgenes[[1]],]
+  deseq_result <- na.omit(deseq_result)
+  
+  #deseq_result <-deseq_result[order(as.numeric(deseq_result[,"stat"])),]
+  
+  #get the normalized count matrix 
+  cdata <- as.data.frame(counts(dds, normalized = TRUE))
+  cdata <- cdata[rownames(cdata) %in% rownames(deseq_result),]
+  
+  #cluster the sample based on phenotype 
+  coldata <- as.data.frame(colData(dds))
+  coldata <- coldata[order(coldata$Condition, decreasing = TRUE),]
+  
+  cdata <- cdata[rownames(coldata)]
+  
+  print("Calculating the expression density...")
+  div <- diversity(cdata, index = "invsimpson")
+  
+  cdata <- cbind(cdata, div)
+  
+  deseq_result <- merge(deseq_result, cdata[c("div")], by = 0, all = FALSE)
+  
+  #calculate the median of diversity for up-regulated genes and down-regulated genes
+  up_div <- median(deseq_result[deseq_result$log2FoldChange > 0 & deseq_result$pvalue <= 0.05,]$div)
+  down_div <- median(deseq_result[deseq_result$log2FoldChange < 0 & deseq_result$pvalue <= 0.05,]$div)
+  
+  #rank the genes
+  deseq_result <-deseq_result[order(as.numeric(deseq_result[,"stat"])),]
+  
+  down_50 <- deseq_result[deseq_result$div > down_div,]$Row.names[1:50]
+  up_50 <- tail(deseq_result[deseq_result$div > up_div,], 50)$Row.names
+  
+  # #cluster the sample based on phenotype 
+  # coldata <- as.data.frame(colData(dds))
+  # coldata <- coldata[order(coldata$Condition, decreasing = TRUE),]
+  
+  
+  df <- cdata[c(down_50, up_50),]
+  #df <- df[rownames(coldata)]
+  
+  df <- df[-ncol(df)]
+  #calculate the z-score across the samples 
+  a <- t(scale(t(df)))
+  
+  pdf(paste0(outpath,'heatmap_',treatment,'_vs_',control,'.pdf'), width = 5, height = 20)
+  pheatmap(a, cluster_cols = FALSE, cluster_rows = FALSE, annotation = coldata["Condition"])
+  dev.off()
+}
+#deseq_result <- res
+#deseq_result <- as.data.frame(deseq_result)
+#
+#deseq_result <- deseq_result[rownames(deseq_result) %in% pgenes$symbol,]
+#deseq_result <- na.omit(deseq_result)
+#
+##deseq_result <-deseq_result[order(as.numeric(deseq_result[,"stat"])),]
+#
+##get the normalized count matrix 
+#cdata <- as.data.frame(counts(dds, normalized = TRUE))
+#cdata <- cdata[rownames(cdata) %in% rownames(deseq_result),]
+#
+##cluster the sample based on phenotype 
+#coldata <- as.data.frame(colData(dds))
+#coldata <- coldata[order(coldata$Condition, decreasing = TRUE),]
+#
+#cdata <- cdata[rownames(coldata)]
+#
+#div <- diversity(cdata, index = "invsimpson")
+#
+#cdata <- cbind(cdata, div)
+#
+#deseq_result <- merge(deseq_result, cdata[c("div")], by = 0, all = FALSE)
+#
+##calculate the median of diversity for up-regulated genes and down-regulated genes
+#up_div <- median(deseq_result[deseq_result$log2FoldChange > 0 & deseq_result$pvalue <= 0.05,]$div)
+#down_div <- median(deseq_result[deseq_result$log2FoldChange < 0 & deseq_result$pvalue <= 0.05,]$div)
+#
+##rank the genes
+#deseq_result <-deseq_result[order(as.numeric(deseq_result[,"stat"])),]
+#
+#down_50 <- deseq_result[deseq_result$div > down_div,]$Row.names[1:50]
+#up_50 <- tail(deseq_result[deseq_result$div > up_div,], 50)$Row.names
+#
+## #cluster the sample based on phenotype 
+## coldata <- as.data.frame(colData(dds))
+## coldata <- coldata[order(coldata$Condition, decreasing = TRUE),]
+#
+#
+#df <- cdata[c(down_50, up_50),]
+##df <- df[rownames(coldata)]
+#
+#df <- df[-ncol(df)]
+##calculate the z-score across the samples 
+#a <- t(scale(t(df)))
+#
+#pdf("/Users/linyang/Documents/Rplot08.pdf", width = 5, height = 20)
+#pheatmap(a, cluster_cols = FALSE, cluster_rows = FALSE, annotation = coldata["Condition"])
+#dev.off()
+#