Commit 25d1f03f authored by mirachirkova's avatar mirachirkova
Browse files

add heat maps

parent d20b4048
install.packages(c('knitr', 'tidyverse', 'hexbin', 'doParallel', 'ggplot2', 'ggbeeswarm', 'ggrepel', 'import', 'logging', 'amap',
install.packages(c('knitr', 'pheatmap', 'tidyverse', 'hexbin', 'doParallel', 'ggplot2', 'ggbeeswarm', 'ggrepel', 'import', 'logging', 'amap',
'grid', 'futile.logger', 'labeling'), dependencies = TRUE)
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
......
......@@ -27,4 +27,4 @@ STAC2
OR51E2
HPN
ACSM1
CLU
RPPH1
......@@ -67,7 +67,7 @@
"STAC2","-","-","-","+","-",1
"OR51E2","-","-","-","+","-",1
"ACSM1","-","-","-","+","-",1
"CLU","-","-","-","+","-",1
"RPPH1","-","-","-","+","-",1
"C2orf88","-","-","-","-","+",1
"APOBEC3C","-","-","-","-","+",1
"CA14","-","-","-","-","+",1
......
No preview for this file type
......@@ -2,4 +2,4 @@ DESeq 0.790475367665374
EBSeq 0.515586137502965
edgeR 0.522764214114926
voom 0.770472302946693
NOISeq 0.762570385003167
NOISeq 0.765683399425221
No preview for this file type
......@@ -123,10 +123,31 @@ loginfo('Save top genes crossing')
crossing(out)
loginfo('Saved')
condition <- coldata$condition
pdf(file.path(out, "heatmaps.pdf"))
# Visualize using heat maps
loginfo('Calculate heat maps for tools')
source("source/calculate_distmatrix_utils.R")
loginfo('Calculate DESeq heat map')
deseq_distmat <- heatmap_v (count, "DESeq", condition, out)
loginfo('DESeq heat map is done')
loginfo('Calculate EBSeq heat map')
ebseq_distmat <- heatmap_v (count, "EBSeq", condition, out)
loginfo('EBSeq heat map is done')
loginfo('Calculate edgeR heat map')
edger_distmat <- heatmap_v (count, "edgeR", condition, out)
loginfo('edgeR heat map is done')
loginfo('Calculate voom heat map')
voom_distmat <- heatmap_v (count, "voom", condition, out)
loginfo('voom heat map is done')
loginfo('Calculate NOISeq heat map')
noiseq_distmat <- heatmap_v (count, "NOISeq", condition, out)
loginfo('NOISeq heat map is done')
loginfo('Heat maps for tools are made')
dev.off()
loginfo('Calculate distance matrices for tools')
# Calculate dist matrices for tools
source("source/calculate_distmatrix_utils.R")
loginfo('Calculate DESeq distance matrix')
deseq_distmat <- calculate_distmatrix (count, file.path(out, "DESeq_sig.txt"))
loginfo('DESeq distance matrix is done')
......
......@@ -54,7 +54,7 @@ deseq2_top <- function(results, n) {
suppressMessages(library("DESeq2"))
suppressMessages(library("biomaRt"))
# Filter results by logFC > 1 or logFC < -1
# Filter results by logFC > 0.25 or logFC < -0.25
filtered_results <- results[!is.na(results$log2FoldChange) > 0 && abs(results$log2FoldChange) > 0.25, ]
# Extract top-n differentially expressed genes ordered by p-value
......
......@@ -23,6 +23,13 @@ if (!file.exists(out)) {
dir.create(out)
}
# Read annotation and make table with useful form
coldata <- read.table(file=cols, sep = ",", dec = ".")
colnames(coldata) <- c("names", "condition")
coldata <- coldata[-c(1), ]
condition <- coldata$condition
n <- 30
# Import functions from modules
source("source/DESeq.R")
source("source/EBSeq.R")
......@@ -36,34 +43,34 @@ edger_res <- readRDS(file = file.path(out, "edgeR_res.rds"))
voom_res <- readRDS(file = file.path(out, "voom_res.rds"))
noiseq_res <- readRDS(file = file.path(out, "NOISeq_res.rds"))
# Visualize results of differential expression
#loginfo('Visualize results of differential expression')
h_results <- read.table(file=file.path(out, "hobotnica_scores.txt"), sep = " ", dec = ".")
colnames(h_results) <- c("names", "scores")
h_results <- h_results[order(h_results$scores, decreasing = TRUE), ]
pdf(file.path(out, "heatmaps.pdf"))
# Visualize using heat maps
loginfo('Calculate heat maps for tools')
source("source/calculate_distmatrix_utils.R")
loginfo('Calculate DESeq heat map')
deseq_distmat <- heatmap_v (count, "DESeq", condition, out)
loginfo('DESeq heat map is done')
loginfo('Calculate EBSeq heat map')
ebseq_distmat <- heatmap_v (count, "EBSeq", condition, out)
loginfo('EBSeq heat map is done')
loginfo('Calculate edgeR heat map')
edger_distmat <- heatmap_v (count, "edgeR", condition, out)
loginfo('edgeR heat map is done')
loginfo('Calculate voom heat map')
voom_distmat <- heatmap_v (count, "voom", condition, out)
loginfo('voom heat map is done')
loginfo('Calculate NOISeq heat map')
noiseq_distmat <- heatmap_v (count, "NOISeq", condition, out)
loginfo('NOISeq heat map is done')
loginfo('Heat maps for tools are made')
dev.off()
colnames(cross_results) <- c("DESeq", "EBSeq", "edgeR", "NOISeq", "voom"
loginfo('The end!')
#pdf(file.path(out, "de_plots.pdf"))
library(gridExtra)
library(grid)
library(ggplot2)
library(lattice)
deseq_plot <- deseq2_v(deseq2_res, out)
ebseq_plot <- ebseq_v(ebseq_res, out)
edger_plot <- edger_v(edger_res, out)
voom_plot <- voom_v(voom_res, out)
noiseq_plot <- noiseq_v(noiseq_res, out)
plots <- data.frame(deseq_plot, ebseq_plot, edger_plot,
voom_plot, noiseq_plot)
colnames(cross_results) <- c("DESeq", "EBSeq", "edgeR", "NOISeq", "voom", "count")
#grid.arrange(b, c, ncol=2, top = textGrob('EnhancedVolcano',
# just = c('center'),
# gp = gpar(fontsize = 32)))
#deseq2_v(deseq2_res, out)
#ebseq_v(ebseq_res, out)
#edger_v(edger_res, out)
#voom_v(voom_res, out)
#noiseq_v(noiseq_res, out)
#dev.off()
loginfo('Visualization results of differential expression is done')
#loginfo('Visualization results of differential expression is done')
......@@ -16,4 +16,28 @@ calculate_distmatrix <- function(countMatrixFile, signatureFile) {
distMatrix <- Dist(t(cm_subset), method="kendall", nbproc=6)
return(distMatrix)
}
heatmap_v <- function(countMatrixFile, tool_name, condition, out) {
signatureFile <- file.path(out, paste0(tool_name, "_sig.txt"))
signature <- readLines(signatureFile)
cm <- read.table(countMatrixFile, header=T, sep=",")
cm <- read.table(count, header=T, sep=",")
cm <- data.frame(cm[, -1], row.names=cm[, 1])
cm_subset <- cm[signature, ]
my_sample_col <- data.frame(condition)
row.names(my_sample_col) <- colnames(cm_subset)
cm_subset <- cm_subset[, row.names(my_sample_col)]
cal_z_score <- function(x){
(x - min(x)) / (max(x) - min(x))
}
data_subset_norm <- t(apply(cm_subset, 1, cal_z_score))
library('pheatmap')
pheatmap(data_subset_norm, annotation_col = my_sample_col, annotation_names_col = FALSE, annotation_legend = TRUE,
annotation_names_row = FALSE, drop_levels = FALSE, show_rownames = T, show_colnames = F,
cluster_cols = F)
}
\ No newline at end of file
Markdown is supported
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment