library(COTAN)
library(pROC)
options(parallelly.fork.enable = TRUE)
library(Seurat)
library(monocle3)
library(reticulate)
library(stringr)
library(dplyr)
library(tidyr)
library(ggplot2)
library(ggpubr)
<- "Results/FDR/"
dirOut if (!dir.exists(dirOut)) {
dir.create(dirOut)
}
<- "Data/MouseCortexFromLoom/FDR/MergedClusters_For_FDR/" dataSetDir
FDR analysis - Results - thresholds true 5% 20%
Preamble
Dataset composition
<- read.csv(file.path(dataSetDir,"Cells_Usage_DataFrame.csv"),
datasets_csv row.names = 1
)
1:3,] datasets_csv[
Group Collection E13.5.432 E13.5.187 E13.5.434
1 2_Clusters_even_near E13.5-434_+_E15.0-428 0 0 318
2 2_Clusters_even_near E15.0-432_+_E13.5-432 536 0 0
3 2_Clusters_even_near E15.0-508_+_E15.0-509 0 0 0
E13.5.184 E13.5.437 E13.5.510 E15.0.432 E15.0.509 E15.0.510 E15.0.508
1 0 0 0 0 0 0 0
2 0 0 0 536 0 0 0
3 0 0 0 0 397 0 397
E15.0.428 E15.0.434 E15.0.437 E17.5.516 E17.5.505
1 318 0 0 0 0
2 0 0 0 0 0
3 0 0 0 0 0
Define which genes are expressed
For each data set we need to define, independently from the DEA methods, which genes are specific for each cluster. So we need to define first which genes are expressed and which are not expressed. To do so we can take advantage from the fact that we have the original clusters from which the cells were sampled to create the artificial datasets. So looking to the original cluster we define as expressed all genes present in at least the 20% of cells and we define as not expressed the genes completely absent or expressed in less than 5% of cells.
Since these two thresholds can have a big influence on the tools performances we will test also others in other pages.
<- readRDS("Data/MouseCortexFromLoom/FDR/Results/GenePresence_PerCluster.RDS")
file.presence
for (file in list.files("Data/MouseCortexFromLoom/SingleClusterRawData/")) {
# print(file)
<- str_split(file,pattern = "_",simplify = T)[1]
Code <- str_split(Code,pattern = "e",simplify = T)[2]
Time <- str_split(Code,pattern = "e",simplify = T)[1]
Cluster <- str_remove(Cluster,pattern = "Cl")
Cluster <- paste0("E",Time,"-",Cluster)
Cluster <- "Absent"
file.presence[,Cluster] <- readRDS(file.path("Data/MouseCortexFromLoom/SingleClusterRawData/",
dataset.cl
file))<- rowSums(dataset.cl > 0)
number.cell.expressing <- round(0.05*dim(dataset.cl)[2],digits = 0)
AbsentThreshold <- round(0.2*dim(dataset.cl)[2],digits = 0)
PresenceThreshold names(number.cell.expressing[number.cell.expressing > AbsentThreshold]),Cluster] <- "Uncertain"
file.presence[
names(number.cell.expressing[number.cell.expressing >= PresenceThreshold]),Cluster] <- "Present"
file.presence[print(Cluster)
print(table(file.presence[,Cluster]))
}
[1] "E13.5-184"
Absent Present Uncertain
6719 3476 4500
[1] "E13.5-187"
Absent Present Uncertain
5930 4609 4156
[1] "E15.0-428"
Absent Present Uncertain
7470 2864 4361
[1] "E13.5-432"
Absent Present Uncertain
6847 3640 4208
[1] "E15.0-432"
Absent Present Uncertain
6951 3514 4230
[1] "E13.5-434"
Absent Present Uncertain
7191 3185 4319
[1] "E15.0-434"
Absent Present Uncertain
7554 3028 4113
[1] "E13.5-437"
Absent Present Uncertain
7024 3430 4241
[1] "E15.0-437"
Absent Present Uncertain
7090 3419 4186
[1] "E17.5-505"
Absent Present Uncertain
7141 3154 4400
[1] "E15.0-508"
Absent Present Uncertain
6637 3892 4166
[1] "E15.0-509"
Absent Present Uncertain
6495 4039 4161
[1] "E13.5-510"
Absent Present Uncertain
5961 4541 4193
[1] "E15.0-510"
Absent Present Uncertain
6082 4480 4133
[1] "E17.5-516"
Absent Present Uncertain
6931 3496 4268
2 Clusters even
2_Clusters_even_near
True vector
<-datasets_csv[datasets_csv$Group == "2_Clusters_even_near",]
subset.datasets_csv
<- NA
ground_truth_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
<- str_split(subset.datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- file.presence[,clusters]
file.presence.subset
#file.presence.subset <- as.matrix(file.presence.subset)
<- as.data.frame(matrix(nrow = nrow(file.presence.subset),
ground_truth ncol = ncol(file.presence.subset)))
rownames(ground_truth) <- rownames(file.presence.subset)
colnames(ground_truth) <- colnames(file.presence.subset)
== "Absent"] <- 0
ground_truth[file.presence.subset == "Present"] <- 1
ground_truth[file.presence.subset == "Uncertain"] <- 0.6
ground_truth[file.presence.subset
<- ground_truth
file.presence.subset
for (col in 1:ncol(ground_truth)) {
<- FALSE
ground_truth[,col] == 1 & rowMeans(file.presence.subset[,-col,drop = FALSE]) < 0.555 ,col] <- TRUE
ground_truth[file.presence.subset[,col]
}
$genes <- rownames(ground_truth)
ground_truth<- pivot_longer(ground_truth,
ground_truth cols = 1:(ncol(ground_truth)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
ground_truth$set_number <- ind
ground_truth<- rbind(ground_truth_tot, ground_truth)
ground_truth_tot
}<- ground_truth_tot[2:nrow(ground_truth_tot),]
ground_truth_tot
head(ground_truth_tot)
# A tibble: 6 × 5
genes clusters value data_set set_number
<chr> <chr> <lgl> <chr> <int>
1 Neil2 E13.5-434 FALSE 2_Clusters_even_near 1
2 Neil2 E15.0-428 FALSE 2_Clusters_even_near 1
3 Lamc1 E13.5-434 FALSE 2_Clusters_even_near 1
4 Lamc1 E15.0-428 FALSE 2_Clusters_even_near 1
5 Lama1 E13.5-434 FALSE 2_Clusters_even_near 1
6 Lama1 E15.0-428 FALSE 2_Clusters_even_near 1
length(unique(ground_truth_tot$genes))
[1] 14695
sum(ground_truth_tot$value)
[1] 8
ROC for COTAN
<- NA
onlyPositive.pVal.Cotan_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",
file.code $Collection[ind])
subset.datasets_csv<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset <- getClusterizationData(dataset,clName = "mergedClusters")[[2]]
deaCOTAN <- pValueFromDEA(deaCOTAN,
pvalCOTAN numCells = getNumCells(dataset),method = "none")
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization <- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- c(cl.names,
cl.names str_split(names(clusterization[clusterization == cl.val])[1],
pattern = "_",simplify = T)[1])
<- cl.names[!is.na(cl.names)]
cl.names
}colnames(deaCOTAN) <- cl.names
colnames(pvalCOTAN) <- cl.names
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- pvalCOTAN
onlyPositive.pVal.Cotan
for (cl in cl.names) {
print(cl)
#temp.DEA.CotanSign <- deaCOTAN[rownames(pvalCOTAN[pvalCOTAN[,cl] < 0.05,]) ,]
rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl] <- 1 #onlyPositive.pVal.Cotan[rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl]+1
onlyPositive.pVal.Cotan[
}
$genes <- rownames(onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan<- pivot_longer(onlyPositive.pVal.Cotan,
onlyPositive.pVal.Cotan values_to = "p_values",
cols = 1:(ncol(onlyPositive.pVal.Cotan)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
onlyPositive.pVal.Cotan$set_number <- ind
onlyPositive.pVal.Cotan<- rbind(onlyPositive.pVal.Cotan_tot, onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan_tot
}
[1] "E13.5-434"
[1] "E15.0-428"
[1] "E15.0-432"
[1] "E13.5-432"
[1] "E15.0-509"
[1] "E15.0-508"
<- onlyPositive.pVal.Cotan_tot[2:nrow(onlyPositive.pVal.Cotan_tot),]
onlyPositive.pVal.Cotan_tot
<- merge.data.frame(onlyPositive.pVal.Cotan_tot,
onlyPositive.pVal.Cotan_tot by = c("genes","clusters","data_set","set_number"),all.x = T,all.y = F)
ground_truth_tot,
#onlyPositive.pVal.Cotan_tot <- onlyPositive.pVal.Cotan_tot[order(onlyPositive.pVal.Cotan_tot$p_values,
# decreasing = F),]
# df <- as.data.frame(matrix(nrow = nrow(onlyPositive.pVal.Cotan_tot),ncol = 3))
# colnames(df) <- c("TPR","FPR","Method")
# df$Method <- "COTAN"
#
# Positive <- sum(onlyPositive.pVal.Cotan_tot$value)
# Negative <- sum(!onlyPositive.pVal.Cotan_tot$value)
#
# for (i in 1:nrow(onlyPositive.pVal.Cotan_tot)) {
# df[i,"TPR"] <- sum(onlyPositive.pVal.Cotan_tot[1:i,"value"])/Positive
# df[i,"FPR"] <- (i-sum(onlyPositive.pVal.Cotan_tot[1:i,"value"]))/Negative
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(onlyPositive.pVal.Cotan_tot$value)
onlyPositive.pVal.Cotan_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(onlyPositive.pVal.Cotan_tot$value, 1 - onlyPositive.pVal.Cotan_tot$p_values)
roc_resultCOTAN
# Plot the ROC curve
#plot(roc_resultCOTAN)
ROC for Seurat
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val)
roc_resultSeurat
# Plot the ROC curve
#plot(roc_resultSeurat)
ROC for Seurat scTransform
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_ScTransform_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_scTr
ROC for Seurat Bimod
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_Bimod_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_Bimod
ROC for Monocle
<- NA
deaMonocle_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"Monocle_DEA_genes.csv")),row.names = 1)
deaMonocle
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cell_group == cl.val,"cell_group"] <- cl.name
deaMonocle[deaMonocle
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "cell_group",
replacement = "clusters")
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "gene_id",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaMonocle$set_number <- ind
deaMonocle<- as.data.frame(deaMonocle)
deaMonocle <- rbind(deaMonocle_tot, deaMonocle)
deaMonocle_tot
}<- deaMonocle_tot[2:nrow(deaMonocle_tot),]
deaMonocle_tot
<- merge.data.frame(deaMonocle_tot,
deaMonocle_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaMonocle_tot$value)
deaMonocle_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaMonocle_tot$value, 1 - deaMonocle_tot$marker_test_p_value)
roc_resultMonocle
# Plot the ROC curve
#plot(roc_resultMonocle)
ROC from ScamPy
<- NA
deaScamPy_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"ScamPy_DEA_genes.csv")),
deaScamPy row.names = 1)
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$clusters == paste0("cl",cl.val),"clusters"] <- cl.name
deaScamPy[deaScamPy
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
$data_set <- subset.datasets_csv[ind,1]
deaScamPy$set_number <- ind
deaScamPy<- rbind(deaScamPy_tot, deaScamPy)
deaScamPy_tot
}<- deaScamPy_tot[2:nrow(deaScamPy_tot),]
deaScamPy_tot
<- merge.data.frame(deaScamPy_tot,
deaScamPy_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaScamPy_tot$value)
deaScamPy_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaScamPy_tot$value, 1 - deaScamPy_tot$pval)
roc_resultScamPy
# Plot the ROC curve
#plot(roc_resultScamPy)
Summary ROC for all methods
<- ggroc(list(COTAN=roc_resultCOTAN, Seurat=roc_resultSeurat,
TwoClusters_even_near Seurat_scTr = roc_resultSeurat_scTr, Seurat_Bimod = roc_resultSeurat_Bimod,
Monocle=roc_resultMonocle, ScamPy=roc_resultScamPy))
<- TwoClusters_even_near + xlab("FPR") + ylab("TPR")
TwoClusters_even_nearPL TwoClusters_even_nearPL
2_Clusters_even_medium
True vector
<-datasets_csv[datasets_csv$Group == "2_Clusters_even_medium",]
subset.datasets_csv
<- NA
ground_truth_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
<- str_split(subset.datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- file.presence[,clusters]
file.presence.subset
#file.presence.subset <- as.matrix(file.presence.subset)
<- as.data.frame(matrix(nrow = nrow(file.presence.subset),
ground_truth ncol = ncol(file.presence.subset)))
rownames(ground_truth) <- rownames(file.presence.subset)
colnames(ground_truth) <- colnames(file.presence.subset)
== "Absent"] <- 0
ground_truth[file.presence.subset == "Present"] <- 1
ground_truth[file.presence.subset == "Uncertain"] <- 0.6
ground_truth[file.presence.subset
<- ground_truth
file.presence.subset
for (col in 1:ncol(ground_truth)) {
<- FALSE
ground_truth[,col] == 1 & rowMeans(file.presence.subset[,-col,drop = FALSE]) < 0.555 ,col] <- TRUE
ground_truth[file.presence.subset[,col]
}
$genes <- rownames(ground_truth)
ground_truth<- pivot_longer(ground_truth,
ground_truth cols = 1:(ncol(ground_truth)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
ground_truth$set_number <- ind
ground_truth<- rbind(ground_truth_tot, ground_truth)
ground_truth_tot
}<- ground_truth_tot[2:nrow(ground_truth_tot),]
ground_truth_tot
head(ground_truth_tot)
# A tibble: 6 × 5
genes clusters value data_set set_number
<chr> <chr> <lgl> <chr> <int>
1 Neil2 E13.5-187 FALSE 2_Clusters_even_medium 1
2 Neil2 E13.5-184 FALSE 2_Clusters_even_medium 1
3 Lamc1 E13.5-187 FALSE 2_Clusters_even_medium 1
4 Lamc1 E13.5-184 FALSE 2_Clusters_even_medium 1
5 Lama1 E13.5-187 FALSE 2_Clusters_even_medium 1
6 Lama1 E13.5-184 FALSE 2_Clusters_even_medium 1
length(unique(ground_truth_tot$genes))
[1] 14695
sum(ground_truth_tot$value)
[1] 447
ROC for COTAN
<- NA
onlyPositive.pVal.Cotan_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",
file.code $Collection[ind])
subset.datasets_csv<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset <- getClusterizationData(dataset,clName = "mergedClusters")[[2]]
deaCOTAN <- pValueFromDEA(deaCOTAN,
pvalCOTAN numCells = getNumCells(dataset),method = "none")
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization <- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- c(cl.names,
cl.names str_split(names(clusterization[clusterization == cl.val])[1],
pattern = "_",simplify = T)[1])
<- cl.names[!is.na(cl.names)]
cl.names
}colnames(deaCOTAN) <- cl.names
colnames(pvalCOTAN) <- cl.names
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- pvalCOTAN
onlyPositive.pVal.Cotan
for (cl in cl.names) {
print(cl)
#temp.DEA.CotanSign <- deaCOTAN[rownames(pvalCOTAN[pvalCOTAN[,cl] < 0.05,]) ,]
rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl] <- 1 #onlyPositive.pVal.Cotan[rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl]+1
onlyPositive.pVal.Cotan[
}
$genes <- rownames(onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan<- pivot_longer(onlyPositive.pVal.Cotan,
onlyPositive.pVal.Cotan values_to = "p_values",
cols = 1:(ncol(onlyPositive.pVal.Cotan)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
onlyPositive.pVal.Cotan$set_number <- ind
onlyPositive.pVal.Cotan<- rbind(onlyPositive.pVal.Cotan_tot, onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan_tot
}
[1] "E13.5-187"
[1] "E13.5-184"
[1] "E17.5-516"
[1] "E15.0-434"
[1] "E15.0-508"
[1] "E15.0-437"
<- onlyPositive.pVal.Cotan_tot[2:nrow(onlyPositive.pVal.Cotan_tot),]
onlyPositive.pVal.Cotan_tot
<- merge.data.frame(onlyPositive.pVal.Cotan_tot,
onlyPositive.pVal.Cotan_tot by = c("genes","clusters","data_set","set_number"),all.x = T,all.y = F)
ground_truth_tot,
#onlyPositive.pVal.Cotan_tot <- onlyPositive.pVal.Cotan_tot[order(onlyPositive.pVal.Cotan_tot$p_values,
# decreasing = F),]
# df <- as.data.frame(matrix(nrow = nrow(onlyPositive.pVal.Cotan_tot),ncol = 3))
# colnames(df) <- c("TPR","FPR","Method")
# df$Method <- "COTAN"
#
# Positive <- sum(onlyPositive.pVal.Cotan_tot$value)
# Negative <- sum(!onlyPositive.pVal.Cotan_tot$value)
#
# for (i in 1:nrow(onlyPositive.pVal.Cotan_tot)) {
# df[i,"TPR"] <- sum(onlyPositive.pVal.Cotan_tot[1:i,"value"])/Positive
# df[i,"FPR"] <- (i-sum(onlyPositive.pVal.Cotan_tot[1:i,"value"]))/Negative
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(onlyPositive.pVal.Cotan_tot$value)
onlyPositive.pVal.Cotan_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(onlyPositive.pVal.Cotan_tot$value, 1 - onlyPositive.pVal.Cotan_tot$p_values)
roc_resultCOTAN
# Plot the ROC curve
#plot(roc_resultCOTAN)
ROC for Seurat
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val)
roc_resultSeurat
# Plot the ROC curve
#plot(roc_resultSeurat)
ROC for Seurat scTransform
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_ScTransform_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_scTr
ROC for Seurat Bimod
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_Bimod_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_Bimod
ROC for Monocle
<- NA
deaMonocle_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"Monocle_DEA_genes.csv")),row.names = 1)
deaMonocle
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cell_group == cl.val,"cell_group"] <- cl.name
deaMonocle[deaMonocle
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "cell_group",
replacement = "clusters")
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "gene_id",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaMonocle$set_number <- ind
deaMonocle<- as.data.frame(deaMonocle)
deaMonocle <- rbind(deaMonocle_tot, deaMonocle)
deaMonocle_tot
}<- deaMonocle_tot[2:nrow(deaMonocle_tot),]
deaMonocle_tot
<- merge.data.frame(deaMonocle_tot,
deaMonocle_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaMonocle_tot$value)
deaMonocle_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaMonocle_tot$value, 1 - deaMonocle_tot$marker_test_p_value)
roc_resultMonocle
# Plot the ROC curve
#plot(roc_resultMonocle)
ROC from ScamPy
<- NA
deaScamPy_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"ScamPy_DEA_genes.csv")),
deaScamPy row.names = 1)
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$clusters == paste0("cl",cl.val),"clusters"] <- cl.name
deaScamPy[deaScamPy
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
$data_set <- subset.datasets_csv[ind,1]
deaScamPy$set_number <- ind
deaScamPy<- rbind(deaScamPy_tot, deaScamPy)
deaScamPy_tot
}<- deaScamPy_tot[2:nrow(deaScamPy_tot),]
deaScamPy_tot
<- merge.data.frame(deaScamPy_tot,
deaScamPy_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaScamPy_tot$value)
deaScamPy_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaScamPy_tot$value, 1 - deaScamPy_tot$pval)
roc_resultScamPy
# Plot the ROC curve
#plot(roc_resultScamPy)
Summary ROC for all methods
<- ggroc(list(COTAN=roc_resultCOTAN, Seurat=roc_resultSeurat, Seurat_scTr = roc_resultSeurat_scTr, Seurat_Bimod = roc_resultSeurat_Bimod,
Two_Clusters_even_medium Monocle=roc_resultMonocle, ScamPy=roc_resultScamPy))
<- Two_Clusters_even_medium + xlab("FPR") + ylab("TPR")
Two_Clusters_even_mediumPL Two_Clusters_even_mediumPL
2_Clusters_even_far
True vector
<-datasets_csv[datasets_csv$Group == "2_Clusters_even_far",]
subset.datasets_csv
<- NA
ground_truth_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
<- str_split(subset.datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- file.presence[,clusters]
file.presence.subset
#file.presence.subset <- as.matrix(file.presence.subset)
<- as.data.frame(matrix(nrow = nrow(file.presence.subset),
ground_truth ncol = ncol(file.presence.subset)))
rownames(ground_truth) <- rownames(file.presence.subset)
colnames(ground_truth) <- colnames(file.presence.subset)
== "Absent"] <- 0
ground_truth[file.presence.subset == "Present"] <- 1
ground_truth[file.presence.subset == "Uncertain"] <- 0.6
ground_truth[file.presence.subset
<- ground_truth
file.presence.subset
for (col in 1:ncol(ground_truth)) {
<- FALSE
ground_truth[,col] == 1 & rowMeans(file.presence.subset[,-col,drop = FALSE]) < 0.555 ,col] <- TRUE
ground_truth[file.presence.subset[,col]
}
$genes <- rownames(ground_truth)
ground_truth<- pivot_longer(ground_truth,
ground_truth cols = 1:(ncol(ground_truth)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
ground_truth$set_number <- ind
ground_truth<- rbind(ground_truth_tot, ground_truth)
ground_truth_tot
}<- ground_truth_tot[2:nrow(ground_truth_tot),]
ground_truth_tot
head(ground_truth_tot)
# A tibble: 6 × 5
genes clusters value data_set set_number
<chr> <chr> <lgl> <chr> <int>
1 Neil2 E17.5-516 FALSE 2_Clusters_even_far 1
2 Neil2 E13.5-187 FALSE 2_Clusters_even_far 1
3 Lamc1 E17.5-516 FALSE 2_Clusters_even_far 1
4 Lamc1 E13.5-187 FALSE 2_Clusters_even_far 1
5 Lama1 E17.5-516 FALSE 2_Clusters_even_far 1
6 Lama1 E13.5-187 FALSE 2_Clusters_even_far 1
length(unique(ground_truth_tot$genes))
[1] 14695
sum(ground_truth_tot$value)
[1] 1359
ROC for COTAN
<- NA
onlyPositive.pVal.Cotan_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",
file.code $Collection[ind])
subset.datasets_csv<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset <- getClusterizationData(dataset,clName = "mergedClusters")[[2]]
deaCOTAN <- pValueFromDEA(deaCOTAN,
pvalCOTAN numCells = getNumCells(dataset),method = "none")
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization <- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- c(cl.names,
cl.names str_split(names(clusterization[clusterization == cl.val])[1],
pattern = "_",simplify = T)[1])
<- cl.names[!is.na(cl.names)]
cl.names
}colnames(deaCOTAN) <- cl.names
colnames(pvalCOTAN) <- cl.names
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- pvalCOTAN
onlyPositive.pVal.Cotan
for (cl in cl.names) {
print(cl)
#temp.DEA.CotanSign <- deaCOTAN[rownames(pvalCOTAN[pvalCOTAN[,cl] < 0.05,]) ,]
rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl] <- 1 #onlyPositive.pVal.Cotan[rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl]+1
onlyPositive.pVal.Cotan[
}
$genes <- rownames(onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan<- pivot_longer(onlyPositive.pVal.Cotan,
onlyPositive.pVal.Cotan values_to = "p_values",
cols = 1:(ncol(onlyPositive.pVal.Cotan)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
onlyPositive.pVal.Cotan$set_number <- ind
onlyPositive.pVal.Cotan<- rbind(onlyPositive.pVal.Cotan_tot, onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan_tot
}
[1] "E13.5-187"
[1] "E17.5-516"
[1] "E15.0-510"
[1] "E13.5-437"
[1] "E15.0-509"
[1] "E13.5-184"
<- onlyPositive.pVal.Cotan_tot[2:nrow(onlyPositive.pVal.Cotan_tot),]
onlyPositive.pVal.Cotan_tot
<- merge.data.frame(onlyPositive.pVal.Cotan_tot,
onlyPositive.pVal.Cotan_tot by = c("genes","clusters","data_set","set_number"),all.x = T,all.y = F)
ground_truth_tot,
#onlyPositive.pVal.Cotan_tot <- onlyPositive.pVal.Cotan_tot[order(onlyPositive.pVal.Cotan_tot$p_values,
# decreasing = F),]
# df <- as.data.frame(matrix(nrow = nrow(onlyPositive.pVal.Cotan_tot),ncol = 3))
# colnames(df) <- c("TPR","FPR","Method")
# df$Method <- "COTAN"
#
# Positive <- sum(onlyPositive.pVal.Cotan_tot$value)
# Negative <- sum(!onlyPositive.pVal.Cotan_tot$value)
#
# for (i in 1:nrow(onlyPositive.pVal.Cotan_tot)) {
# df[i,"TPR"] <- sum(onlyPositive.pVal.Cotan_tot[1:i,"value"])/Positive
# df[i,"FPR"] <- (i-sum(onlyPositive.pVal.Cotan_tot[1:i,"value"]))/Negative
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(onlyPositive.pVal.Cotan_tot$value)
onlyPositive.pVal.Cotan_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(onlyPositive.pVal.Cotan_tot$value, 1 - onlyPositive.pVal.Cotan_tot$p_values)
roc_resultCOTAN
# Plot the ROC curve
#plot(roc_resultCOTAN)
ROC for Seurat
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val)
roc_resultSeurat
# Plot the ROC curve
#plot(roc_resultSeurat)
ROC for Seurat scTransform
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_ScTransform_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_scTr
ROC for Seurat Bimod
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_Bimod_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_Bimod
ROC for Monocle
<- NA
deaMonocle_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"Monocle_DEA_genes.csv")),row.names = 1)
deaMonocle
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cell_group == cl.val,"cell_group"] <- cl.name
deaMonocle[deaMonocle
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "cell_group",
replacement = "clusters")
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "gene_id",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaMonocle$set_number <- ind
deaMonocle<- as.data.frame(deaMonocle)
deaMonocle <- rbind(deaMonocle_tot, deaMonocle)
deaMonocle_tot
}<- deaMonocle_tot[2:nrow(deaMonocle_tot),]
deaMonocle_tot
<- merge.data.frame(deaMonocle_tot,
deaMonocle_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaMonocle_tot$value)
deaMonocle_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaMonocle_tot$value, 1 - deaMonocle_tot$marker_test_p_value)
roc_resultMonocle
# Plot the ROC curve
#plot(roc_resultMonocle)
ROC from ScamPy
<- NA
deaScamPy_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"ScamPy_DEA_genes.csv")),
deaScamPy row.names = 1)
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$clusters == paste0("cl",cl.val),"clusters"] <- cl.name
deaScamPy[deaScamPy
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
$data_set <- subset.datasets_csv[ind,1]
deaScamPy$set_number <- ind
deaScamPy<- rbind(deaScamPy_tot, deaScamPy)
deaScamPy_tot
}<- deaScamPy_tot[2:nrow(deaScamPy_tot),]
deaScamPy_tot
<- merge.data.frame(deaScamPy_tot,
deaScamPy_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaScamPy_tot$value)
deaScamPy_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaScamPy_tot$value, 1 - deaScamPy_tot$pval)
roc_resultScamPy
# Plot the ROC curve
#plot(roc_resultScamPy)
Summary ROC for all methods
<- ggroc(list(COTAN=roc_resultCOTAN, Seurat=roc_resultSeurat, Seurat_scTr = roc_resultSeurat_scTr, Seurat_Bimod = roc_resultSeurat_Bimod,
Two_Clusters_even_far Monocle=roc_resultMonocle, ScamPy=roc_resultScamPy))
<- Two_Clusters_even_far + xlab("FPR") + ylab("TPR")
Two_Clusters_even_farPL Two_Clusters_even_farPL
2 clusters uneven
2_Clusters_uneven_near
True vector
<-datasets_csv[datasets_csv$Group == "2_Clusters_uneven_near",]
subset.datasets_csv
<- NA
ground_truth_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
<- str_split(subset.datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- file.presence[,clusters]
file.presence.subset
#file.presence.subset <- as.matrix(file.presence.subset)
<- as.data.frame(matrix(nrow = nrow(file.presence.subset),
ground_truth ncol = ncol(file.presence.subset)))
rownames(ground_truth) <- rownames(file.presence.subset)
colnames(ground_truth) <- colnames(file.presence.subset)
== "Absent"] <- 0
ground_truth[file.presence.subset == "Present"] <- 1
ground_truth[file.presence.subset == "Uncertain"] <- 0.6
ground_truth[file.presence.subset
<- ground_truth
file.presence.subset
for (col in 1:ncol(ground_truth)) {
<- FALSE
ground_truth[,col] == 1 & rowMeans(file.presence.subset[,-col,drop = FALSE]) < 0.555 ,col] <- TRUE
ground_truth[file.presence.subset[,col]
}
$genes <- rownames(ground_truth)
ground_truth<- pivot_longer(ground_truth,
ground_truth cols = 1:(ncol(ground_truth)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
ground_truth$set_number <- ind
ground_truth<- rbind(ground_truth_tot, ground_truth)
ground_truth_tot
}<- ground_truth_tot[2:nrow(ground_truth_tot),]
ground_truth_tot
head(ground_truth_tot)
# A tibble: 6 × 5
genes clusters value data_set set_number
<chr> <chr> <lgl> <chr> <int>
1 Neil2 E13.5-434 FALSE 2_Clusters_uneven_near 1
2 Neil2 E15.0-428 FALSE 2_Clusters_uneven_near 1
3 Lamc1 E13.5-434 FALSE 2_Clusters_uneven_near 1
4 Lamc1 E15.0-428 FALSE 2_Clusters_uneven_near 1
5 Lama1 E13.5-434 FALSE 2_Clusters_uneven_near 1
6 Lama1 E15.0-428 FALSE 2_Clusters_uneven_near 1
length(unique(ground_truth_tot$genes))
[1] 14695
sum(ground_truth_tot$value)
[1] 8
ROC for COTAN
<- NA
onlyPositive.pVal.Cotan_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",
file.code $Collection[ind])
subset.datasets_csv<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset <- getClusterizationData(dataset,clName = "mergedClusters")[[2]]
deaCOTAN <- pValueFromDEA(deaCOTAN,
pvalCOTAN numCells = getNumCells(dataset),method = "none")
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization <- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- c(cl.names,
cl.names str_split(names(clusterization[clusterization == cl.val])[1],
pattern = "_",simplify = T)[1])
<- cl.names[!is.na(cl.names)]
cl.names
}colnames(deaCOTAN) <- cl.names
colnames(pvalCOTAN) <- cl.names
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- pvalCOTAN
onlyPositive.pVal.Cotan
for (cl in cl.names) {
print(cl)
#temp.DEA.CotanSign <- deaCOTAN[rownames(pvalCOTAN[pvalCOTAN[,cl] < 0.05,]) ,]
rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl] <- 1 #onlyPositive.pVal.Cotan[rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl]+1
onlyPositive.pVal.Cotan[
}
$genes <- rownames(onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan<- pivot_longer(onlyPositive.pVal.Cotan,
onlyPositive.pVal.Cotan values_to = "p_values",
cols = 1:(ncol(onlyPositive.pVal.Cotan)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
onlyPositive.pVal.Cotan$set_number <- ind
onlyPositive.pVal.Cotan<- rbind(onlyPositive.pVal.Cotan_tot, onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan_tot
}
[1] "E13.5-434"
[1] "E15.0-428"
[1] "E15.0-432"
[1] "E13.5-432"
[1] "E15.0-509"
[1] "E15.0-508"
<- onlyPositive.pVal.Cotan_tot[2:nrow(onlyPositive.pVal.Cotan_tot),]
onlyPositive.pVal.Cotan_tot
<- merge.data.frame(onlyPositive.pVal.Cotan_tot,
onlyPositive.pVal.Cotan_tot by = c("genes","clusters","data_set","set_number"),all.x = T,all.y = F)
ground_truth_tot,
#onlyPositive.pVal.Cotan_tot <- onlyPositive.pVal.Cotan_tot[order(onlyPositive.pVal.Cotan_tot$p_values,
# decreasing = F),]
# df <- as.data.frame(matrix(nrow = nrow(onlyPositive.pVal.Cotan_tot),ncol = 3))
# colnames(df) <- c("TPR","FPR","Method")
# df$Method <- "COTAN"
#
# Positive <- sum(onlyPositive.pVal.Cotan_tot$value)
# Negative <- sum(!onlyPositive.pVal.Cotan_tot$value)
#
# for (i in 1:nrow(onlyPositive.pVal.Cotan_tot)) {
# df[i,"TPR"] <- sum(onlyPositive.pVal.Cotan_tot[1:i,"value"])/Positive
# df[i,"FPR"] <- (i-sum(onlyPositive.pVal.Cotan_tot[1:i,"value"]))/Negative
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(onlyPositive.pVal.Cotan_tot$value)
onlyPositive.pVal.Cotan_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(onlyPositive.pVal.Cotan_tot$value, 1 - onlyPositive.pVal.Cotan_tot$p_values)
roc_resultCOTAN
# Plot the ROC curve
#plot(roc_resultCOTAN)
ROC for Seurat
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val)
roc_resultSeurat
# Plot the ROC curve
plot(roc_resultSeurat)
ROC for Seurat scTransform
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_ScTransform_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_scTr
ROC for Seurat Bimod
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_Bimod_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_Bimod
ROC for Monocle
<- NA
deaMonocle_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"Monocle_DEA_genes.csv")),row.names = 1)
deaMonocle
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cell_group == cl.val,"cell_group"] <- cl.name
deaMonocle[deaMonocle
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "cell_group",
replacement = "clusters")
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "gene_id",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaMonocle$set_number <- ind
deaMonocle<- as.data.frame(deaMonocle)
deaMonocle <- rbind(deaMonocle_tot, deaMonocle)
deaMonocle_tot
}<- deaMonocle_tot[2:nrow(deaMonocle_tot),]
deaMonocle_tot
<- merge.data.frame(deaMonocle_tot,
deaMonocle_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaMonocle_tot$value)
deaMonocle_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaMonocle_tot$value, 1 - deaMonocle_tot$marker_test_p_value)
roc_resultMonocle
# Plot the ROC curve
#plot(roc_resultMonocle)
ROC from ScamPy
<- NA
deaScamPy_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"ScamPy_DEA_genes.csv")),
deaScamPy row.names = 1)
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$clusters == paste0("cl",cl.val),"clusters"] <- cl.name
deaScamPy[deaScamPy
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
$data_set <- subset.datasets_csv[ind,1]
deaScamPy$set_number <- ind
deaScamPy<- rbind(deaScamPy_tot, deaScamPy)
deaScamPy_tot
}<- deaScamPy_tot[2:nrow(deaScamPy_tot),]
deaScamPy_tot
<- merge.data.frame(deaScamPy_tot,
deaScamPy_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaScamPy_tot$value)
deaScamPy_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaScamPy_tot$value, 1 - deaScamPy_tot$pval)
roc_resultScamPy
# Plot the ROC curve
#plot(roc_resultScamPy)
Summary ROC for all methods
<- ggroc(list(COTAN=roc_resultCOTAN, Seurat=roc_resultSeurat, Seurat_scTr = roc_resultSeurat_scTr, Seurat_Bimod = roc_resultSeurat_Bimod,
Two_Clusters_uneven_near Monocle=roc_resultMonocle, ScamPy=roc_resultScamPy))
<- Two_Clusters_uneven_near + xlab("FPR") + ylab("TPR")
Two_Clusters_uneven_nearPL Two_Clusters_uneven_nearPL
2_Clusters_uneven_medium
True vector
<-datasets_csv[datasets_csv$Group == "2_Clusters_uneven_medium",]
subset.datasets_csv
<- NA
ground_truth_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
<- str_split(subset.datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- file.presence[,clusters]
file.presence.subset
#file.presence.subset <- as.matrix(file.presence.subset)
<- as.data.frame(matrix(nrow = nrow(file.presence.subset),
ground_truth ncol = ncol(file.presence.subset)))
rownames(ground_truth) <- rownames(file.presence.subset)
colnames(ground_truth) <- colnames(file.presence.subset)
== "Absent"] <- 0
ground_truth[file.presence.subset == "Present"] <- 1
ground_truth[file.presence.subset == "Uncertain"] <- 0.6
ground_truth[file.presence.subset
<- ground_truth
file.presence.subset
for (col in 1:ncol(ground_truth)) {
<- FALSE
ground_truth[,col] == 1 & rowMeans(file.presence.subset[,-col,drop = FALSE]) < 0.555 ,col] <- TRUE
ground_truth[file.presence.subset[,col]
}
$genes <- rownames(ground_truth)
ground_truth<- pivot_longer(ground_truth,
ground_truth cols = 1:(ncol(ground_truth)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
ground_truth$set_number <- ind
ground_truth<- rbind(ground_truth_tot, ground_truth)
ground_truth_tot
}<- ground_truth_tot[2:nrow(ground_truth_tot),]
ground_truth_tot
head(ground_truth_tot)
# A tibble: 6 × 5
genes clusters value data_set set_number
<chr> <chr> <lgl> <chr> <int>
1 Neil2 E13.5-187 FALSE 2_Clusters_uneven_medium 1
2 Neil2 E13.5-184 FALSE 2_Clusters_uneven_medium 1
3 Lamc1 E13.5-187 FALSE 2_Clusters_uneven_medium 1
4 Lamc1 E13.5-184 FALSE 2_Clusters_uneven_medium 1
5 Lama1 E13.5-187 FALSE 2_Clusters_uneven_medium 1
6 Lama1 E13.5-184 FALSE 2_Clusters_uneven_medium 1
length(unique(ground_truth_tot$genes))
[1] 14695
sum(ground_truth_tot$value)
[1] 447
ROC for COTAN
<- NA
onlyPositive.pVal.Cotan_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",
file.code $Collection[ind])
subset.datasets_csv<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset <- getClusterizationData(dataset,clName = "mergedClusters")[[2]]
deaCOTAN <- pValueFromDEA(deaCOTAN,
pvalCOTAN numCells = getNumCells(dataset),method = "none")
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization <- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- c(cl.names,
cl.names str_split(names(clusterization[clusterization == cl.val])[1],
pattern = "_",simplify = T)[1])
<- cl.names[!is.na(cl.names)]
cl.names
}colnames(deaCOTAN) <- cl.names
colnames(pvalCOTAN) <- cl.names
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- pvalCOTAN
onlyPositive.pVal.Cotan
for (cl in cl.names) {
print(cl)
#temp.DEA.CotanSign <- deaCOTAN[rownames(pvalCOTAN[pvalCOTAN[,cl] < 0.05,]) ,]
rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl] <- 1 #onlyPositive.pVal.Cotan[rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl]+1
onlyPositive.pVal.Cotan[
}
$genes <- rownames(onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan<- pivot_longer(onlyPositive.pVal.Cotan,
onlyPositive.pVal.Cotan values_to = "p_values",
cols = 1:(ncol(onlyPositive.pVal.Cotan)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
onlyPositive.pVal.Cotan$set_number <- ind
onlyPositive.pVal.Cotan<- rbind(onlyPositive.pVal.Cotan_tot, onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan_tot
}
[1] "E13.5-187"
[1] "E13.5-184"
[1] "E17.5-516"
[1] "E15.0-434"
[1] "E15.0-508"
[1] "E15.0-437"
<- onlyPositive.pVal.Cotan_tot[2:nrow(onlyPositive.pVal.Cotan_tot),]
onlyPositive.pVal.Cotan_tot
<- merge.data.frame(onlyPositive.pVal.Cotan_tot,
onlyPositive.pVal.Cotan_tot by = c("genes","clusters","data_set","set_number"),all.x = T,all.y = F)
ground_truth_tot,
#onlyPositive.pVal.Cotan_tot <- onlyPositive.pVal.Cotan_tot[order(onlyPositive.pVal.Cotan_tot$p_values,
# decreasing = F),]
# df <- as.data.frame(matrix(nrow = nrow(onlyPositive.pVal.Cotan_tot),ncol = 3))
# colnames(df) <- c("TPR","FPR","Method")
# df$Method <- "COTAN"
#
# Positive <- sum(onlyPositive.pVal.Cotan_tot$value)
# Negative <- sum(!onlyPositive.pVal.Cotan_tot$value)
#
# for (i in 1:nrow(onlyPositive.pVal.Cotan_tot)) {
# df[i,"TPR"] <- sum(onlyPositive.pVal.Cotan_tot[1:i,"value"])/Positive
# df[i,"FPR"] <- (i-sum(onlyPositive.pVal.Cotan_tot[1:i,"value"]))/Negative
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(onlyPositive.pVal.Cotan_tot$value)
onlyPositive.pVal.Cotan_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(onlyPositive.pVal.Cotan_tot$value, 1 - onlyPositive.pVal.Cotan_tot$p_values)
roc_resultCOTAN
# Plot the ROC curve
#plot(roc_resultCOTAN)
ROC for Seurat
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val)
roc_resultSeurat
# Plot the ROC curve
#plot(roc_resultSeurat)
ROC for Seurat scTransform
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_ScTransform_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_scTr
ROC for Seurat Bimod
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_Bimod_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_Bimod
ROC for Monocle
<- NA
deaMonocle_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"Monocle_DEA_genes.csv")),row.names = 1)
deaMonocle
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cell_group == cl.val,"cell_group"] <- cl.name
deaMonocle[deaMonocle
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "cell_group",
replacement = "clusters")
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "gene_id",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaMonocle$set_number <- ind
deaMonocle<- as.data.frame(deaMonocle)
deaMonocle <- rbind(deaMonocle_tot, deaMonocle)
deaMonocle_tot
}<- deaMonocle_tot[2:nrow(deaMonocle_tot),]
deaMonocle_tot
<- merge.data.frame(deaMonocle_tot,
deaMonocle_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaMonocle_tot$value)
deaMonocle_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaMonocle_tot$value, 1 - deaMonocle_tot$marker_test_p_value)
roc_resultMonocle
# Plot the ROC curve
#plot(roc_resultMonocle)
ROC from ScamPy
<- NA
deaScamPy_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"ScamPy_DEA_genes.csv")),
deaScamPy row.names = 1)
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$clusters == paste0("cl",cl.val),"clusters"] <- cl.name
deaScamPy[deaScamPy
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
$data_set <- subset.datasets_csv[ind,1]
deaScamPy$set_number <- ind
deaScamPy<- rbind(deaScamPy_tot, deaScamPy)
deaScamPy_tot
}<- deaScamPy_tot[2:nrow(deaScamPy_tot),]
deaScamPy_tot
<- merge.data.frame(deaScamPy_tot,
deaScamPy_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaScamPy_tot$value)
deaScamPy_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaScamPy_tot$value, 1 - deaScamPy_tot$pval)
roc_resultScamPy
# Plot the ROC curve
#plot(roc_resultScamPy)
Summary ROC for all methods
<- ggroc(list(COTAN=roc_resultCOTAN, Seurat=roc_resultSeurat, Seurat_scTr = roc_resultSeurat_scTr, Seurat_Bimod = roc_resultSeurat_Bimod,
Two_Clusters_uneven_medium Monocle=roc_resultMonocle, ScamPy=roc_resultScamPy))
<- Two_Clusters_uneven_medium + xlab("FPR") + ylab("TPR")
Two_Clusters_uneven_mediumPL
Two_Clusters_uneven_mediumPL
2_Clusters_uneven_far
True vector
<-datasets_csv[datasets_csv$Group == "2_Clusters_uneven_far",]
subset.datasets_csv
<- NA
ground_truth_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
<- str_split(subset.datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- file.presence[,clusters]
file.presence.subset
#file.presence.subset <- as.matrix(file.presence.subset)
<- as.data.frame(matrix(nrow = nrow(file.presence.subset),
ground_truth ncol = ncol(file.presence.subset)))
rownames(ground_truth) <- rownames(file.presence.subset)
colnames(ground_truth) <- colnames(file.presence.subset)
== "Absent"] <- 0
ground_truth[file.presence.subset == "Present"] <- 1
ground_truth[file.presence.subset == "Uncertain"] <- 0.6
ground_truth[file.presence.subset
<- ground_truth
file.presence.subset
for (col in 1:ncol(ground_truth)) {
<- FALSE
ground_truth[,col] == 1 & rowMeans(file.presence.subset[,-col,drop = FALSE]) < 0.555 ,col] <- TRUE
ground_truth[file.presence.subset[,col]
}
$genes <- rownames(ground_truth)
ground_truth<- pivot_longer(ground_truth,
ground_truth cols = 1:(ncol(ground_truth)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
ground_truth$set_number <- ind
ground_truth<- rbind(ground_truth_tot, ground_truth)
ground_truth_tot
}<- ground_truth_tot[2:nrow(ground_truth_tot),]
ground_truth_tot
head(ground_truth_tot)
# A tibble: 6 × 5
genes clusters value data_set set_number
<chr> <chr> <lgl> <chr> <int>
1 Neil2 E17.5-516 FALSE 2_Clusters_uneven_far 1
2 Neil2 E13.5-187 FALSE 2_Clusters_uneven_far 1
3 Lamc1 E17.5-516 FALSE 2_Clusters_uneven_far 1
4 Lamc1 E13.5-187 FALSE 2_Clusters_uneven_far 1
5 Lama1 E17.5-516 FALSE 2_Clusters_uneven_far 1
6 Lama1 E13.5-187 FALSE 2_Clusters_uneven_far 1
length(unique(ground_truth_tot$genes))
[1] 14695
sum(ground_truth_tot$value)
[1] 1359
ROC for COTAN
<- NA
onlyPositive.pVal.Cotan_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",
file.code $Collection[ind])
subset.datasets_csv<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset <- getClusterizationData(dataset,clName = "mergedClusters")[[2]]
deaCOTAN <- pValueFromDEA(deaCOTAN,
pvalCOTAN numCells = getNumCells(dataset),method = "none")
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization <- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- c(cl.names,
cl.names str_split(names(clusterization[clusterization == cl.val])[1],
pattern = "_",simplify = T)[1])
<- cl.names[!is.na(cl.names)]
cl.names
}colnames(deaCOTAN) <- cl.names
colnames(pvalCOTAN) <- cl.names
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- pvalCOTAN
onlyPositive.pVal.Cotan
for (cl in cl.names) {
print(cl)
#temp.DEA.CotanSign <- deaCOTAN[rownames(pvalCOTAN[pvalCOTAN[,cl] < 0.05,]) ,]
rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl] <- 1 #onlyPositive.pVal.Cotan[rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl]+1
onlyPositive.pVal.Cotan[
}
$genes <- rownames(onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan<- pivot_longer(onlyPositive.pVal.Cotan,
onlyPositive.pVal.Cotan values_to = "p_values",
cols = 1:(ncol(onlyPositive.pVal.Cotan)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
onlyPositive.pVal.Cotan$set_number <- ind
onlyPositive.pVal.Cotan<- rbind(onlyPositive.pVal.Cotan_tot, onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan_tot
}
[1] "E13.5-187"
[1] "E17.5-516"
[1] "E15.0-510"
[1] "E13.5-437"
[1] "E15.0-509"
[1] "E13.5-184"
<- onlyPositive.pVal.Cotan_tot[2:nrow(onlyPositive.pVal.Cotan_tot),]
onlyPositive.pVal.Cotan_tot
<- merge.data.frame(onlyPositive.pVal.Cotan_tot,
onlyPositive.pVal.Cotan_tot by = c("genes","clusters","data_set","set_number"),all.x = T,all.y = F)
ground_truth_tot,
#onlyPositive.pVal.Cotan_tot <- onlyPositive.pVal.Cotan_tot[order(onlyPositive.pVal.Cotan_tot$p_values,
# decreasing = F),]
# df <- as.data.frame(matrix(nrow = nrow(onlyPositive.pVal.Cotan_tot),ncol = 3))
# colnames(df) <- c("TPR","FPR","Method")
# df$Method <- "COTAN"
#
# Positive <- sum(onlyPositive.pVal.Cotan_tot$value)
# Negative <- sum(!onlyPositive.pVal.Cotan_tot$value)
#
# for (i in 1:nrow(onlyPositive.pVal.Cotan_tot)) {
# df[i,"TPR"] <- sum(onlyPositive.pVal.Cotan_tot[1:i,"value"])/Positive
# df[i,"FPR"] <- (i-sum(onlyPositive.pVal.Cotan_tot[1:i,"value"]))/Negative
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(onlyPositive.pVal.Cotan_tot$value)
onlyPositive.pVal.Cotan_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(onlyPositive.pVal.Cotan_tot$value, 1 - onlyPositive.pVal.Cotan_tot$p_values)
roc_resultCOTAN
# Plot the ROC curve
#plot(roc_resultCOTAN)
ROC for Seurat
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val)
roc_resultSeurat
# Plot the ROC curve
#plot(roc_resultSeurat)
ROC for Seurat scTransform
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_ScTransform_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_scTr
ROC for Seurat Bimod
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_Bimod_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_Bimod
ROC for Monocle
<- NA
deaMonocle_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"Monocle_DEA_genes.csv")),row.names = 1)
deaMonocle
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cell_group == cl.val,"cell_group"] <- cl.name
deaMonocle[deaMonocle
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "cell_group",
replacement = "clusters")
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "gene_id",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaMonocle$set_number <- ind
deaMonocle<- as.data.frame(deaMonocle)
deaMonocle <- rbind(deaMonocle_tot, deaMonocle)
deaMonocle_tot
}<- deaMonocle_tot[2:nrow(deaMonocle_tot),]
deaMonocle_tot
<- merge.data.frame(deaMonocle_tot,
deaMonocle_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaMonocle_tot$value)
deaMonocle_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaMonocle_tot$value, 1 - deaMonocle_tot$marker_test_p_value)
roc_resultMonocle
# Plot the ROC curve
#plot(roc_resultMonocle)
ROC from ScamPy
<- NA
deaScamPy_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"ScamPy_DEA_genes.csv")),
deaScamPy row.names = 1)
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$clusters == paste0("cl",cl.val),"clusters"] <- cl.name
deaScamPy[deaScamPy
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
$data_set <- subset.datasets_csv[ind,1]
deaScamPy$set_number <- ind
deaScamPy<- rbind(deaScamPy_tot, deaScamPy)
deaScamPy_tot
}<- deaScamPy_tot[2:nrow(deaScamPy_tot),]
deaScamPy_tot
<- merge.data.frame(deaScamPy_tot,
deaScamPy_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaScamPy_tot$value)
deaScamPy_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaScamPy_tot$value, 1 - deaScamPy_tot$pval)
roc_resultScamPy
# Plot the ROC curve
#plot(roc_resultScamPy)
Summary ROC for all methods
<- ggroc(list(COTAN=roc_resultCOTAN, Seurat=roc_resultSeurat, Seurat_scTr = roc_resultSeurat_scTr, Seurat_Bimod = roc_resultSeurat_Bimod,
Two_Clusters_uneven_far Monocle=roc_resultMonocle, ScamPy=roc_resultScamPy))
<- Two_Clusters_uneven_far + xlab("FPR") + ylab("TPR")
Two_Clusters_uneven_farPL
Two_Clusters_uneven_farPL
3 clusters
3_Clusters_even
True vector
<-datasets_csv[datasets_csv$Group == "3_Clusters_even",]
subset.datasets_csv
<- NA
ground_truth_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
<- str_split(subset.datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- file.presence[,clusters]
file.presence.subset
#file.presence.subset <- as.matrix(file.presence.subset)
<- as.data.frame(matrix(nrow = nrow(file.presence.subset),
ground_truth ncol = ncol(file.presence.subset)))
rownames(ground_truth) <- rownames(file.presence.subset)
colnames(ground_truth) <- colnames(file.presence.subset)
== "Absent"] <- 0
ground_truth[file.presence.subset == "Present"] <- 1
ground_truth[file.presence.subset == "Uncertain"] <- 0.6
ground_truth[file.presence.subset
<- ground_truth
file.presence.subset
for (col in 1:ncol(ground_truth)) {
<- FALSE
ground_truth[,col] == 1 & rowMeans(file.presence.subset[,-col,drop = FALSE]) < 0.555 ,col] <- TRUE
ground_truth[file.presence.subset[,col]
}
$genes <- rownames(ground_truth)
ground_truth<- pivot_longer(ground_truth,
ground_truth cols = 1:(ncol(ground_truth)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
ground_truth$set_number <- ind
ground_truth<- rbind(ground_truth_tot, ground_truth)
ground_truth_tot
}<- ground_truth_tot[2:nrow(ground_truth_tot),]
ground_truth_tot
head(ground_truth_tot)
# A tibble: 6 × 5
genes clusters value data_set set_number
<chr> <chr> <lgl> <chr> <int>
1 Neil2 E15.0-437 FALSE 3_Clusters_even 1
2 Neil2 E13.5-510 FALSE 3_Clusters_even 1
3 Neil2 E13.5-437 FALSE 3_Clusters_even 1
4 Lamc1 E15.0-437 FALSE 3_Clusters_even 1
5 Lamc1 E13.5-510 FALSE 3_Clusters_even 1
6 Lamc1 E13.5-437 FALSE 3_Clusters_even 1
length(unique(ground_truth_tot$genes))
[1] 14695
sum(ground_truth_tot$value)
[1] 1305
ROC for COTAN
<- NA
onlyPositive.pVal.Cotan_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",
file.code $Collection[ind])
subset.datasets_csv<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset <- getClusterizationData(dataset,clName = "mergedClusters")[[2]]
deaCOTAN <- pValueFromDEA(deaCOTAN,
pvalCOTAN numCells = getNumCells(dataset),method = "none")
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization <- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- c(cl.names,
cl.names str_split(names(clusterization[clusterization == cl.val])[1],
pattern = "_",simplify = T)[1])
<- cl.names[!is.na(cl.names)]
cl.names
}colnames(deaCOTAN) <- cl.names
colnames(pvalCOTAN) <- cl.names
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- pvalCOTAN
onlyPositive.pVal.Cotan
for (cl in cl.names) {
print(cl)
#temp.DEA.CotanSign <- deaCOTAN[rownames(pvalCOTAN[pvalCOTAN[,cl] < 0.05,]) ,]
rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl] <- 1 #onlyPositive.pVal.Cotan[rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl]+1
onlyPositive.pVal.Cotan[
}
$genes <- rownames(onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan<- pivot_longer(onlyPositive.pVal.Cotan,
onlyPositive.pVal.Cotan values_to = "p_values",
cols = 1:(ncol(onlyPositive.pVal.Cotan)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
onlyPositive.pVal.Cotan$set_number <- ind
onlyPositive.pVal.Cotan<- rbind(onlyPositive.pVal.Cotan_tot, onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan_tot
}
[1] "E13.5-437"
[1] "E15.0-437"
[1] "E13.5-510"
[1] "E17.5-516"
[1] "E13.5-437"
[1] "E17.5-505"
[1] "E15.0-510"
[1] "E15.0-428"
[1] "E13.5-510"
<- onlyPositive.pVal.Cotan_tot[2:nrow(onlyPositive.pVal.Cotan_tot),]
onlyPositive.pVal.Cotan_tot
<- merge.data.frame(onlyPositive.pVal.Cotan_tot,
onlyPositive.pVal.Cotan_tot by = c("genes","clusters","data_set","set_number"),all.x = T,all.y = F)
ground_truth_tot,
#onlyPositive.pVal.Cotan_tot <- onlyPositive.pVal.Cotan_tot[order(onlyPositive.pVal.Cotan_tot$p_values,
# decreasing = F),]
# df <- as.data.frame(matrix(nrow = nrow(onlyPositive.pVal.Cotan_tot),ncol = 3))
# colnames(df) <- c("TPR","FPR","Method")
# df$Method <- "COTAN"
#
# Positive <- sum(onlyPositive.pVal.Cotan_tot$value)
# Negative <- sum(!onlyPositive.pVal.Cotan_tot$value)
#
# for (i in 1:nrow(onlyPositive.pVal.Cotan_tot)) {
# df[i,"TPR"] <- sum(onlyPositive.pVal.Cotan_tot[1:i,"value"])/Positive
# df[i,"FPR"] <- (i-sum(onlyPositive.pVal.Cotan_tot[1:i,"value"]))/Negative
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(onlyPositive.pVal.Cotan_tot$value)
onlyPositive.pVal.Cotan_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(onlyPositive.pVal.Cotan_tot$value, 1 - onlyPositive.pVal.Cotan_tot$p_values)
roc_resultCOTAN
# Plot the ROC curve
#plot(roc_resultCOTAN)
ROC for Seurat
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val)
roc_resultSeurat
# Plot the ROC curve
#plot(roc_resultSeurat)
ROC for Seurat scTransform
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_ScTransform_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_scTr
ROC for Seurat Bimod
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_Bimod_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_Bimod
ROC for Monocle
<- NA
deaMonocle_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"Monocle_DEA_genes.csv")),row.names = 1)
deaMonocle
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cell_group == cl.val,"cell_group"] <- cl.name
deaMonocle[deaMonocle
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "cell_group",
replacement = "clusters")
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "gene_id",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaMonocle$set_number <- ind
deaMonocle<- as.data.frame(deaMonocle)
deaMonocle <- rbind(deaMonocle_tot, deaMonocle)
deaMonocle_tot
}<- deaMonocle_tot[2:nrow(deaMonocle_tot),]
deaMonocle_tot
<- merge.data.frame(deaMonocle_tot,
deaMonocle_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaMonocle_tot$value)
deaMonocle_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaMonocle_tot$value, 1 - deaMonocle_tot$marker_test_p_value)
roc_resultMonocle
# Plot the ROC curve
#plot(roc_resultMonocle)
ROC from ScamPy
<- NA
deaScamPy_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"ScamPy_DEA_genes.csv")),
deaScamPy row.names = 1)
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$clusters == paste0("cl",cl.val),"clusters"] <- cl.name
deaScamPy[deaScamPy
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
$data_set <- subset.datasets_csv[ind,1]
deaScamPy$set_number <- ind
deaScamPy<- rbind(deaScamPy_tot, deaScamPy)
deaScamPy_tot
}<- deaScamPy_tot[2:nrow(deaScamPy_tot),]
deaScamPy_tot
<- merge.data.frame(deaScamPy_tot,
deaScamPy_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaScamPy_tot$value)
deaScamPy_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaScamPy_tot$value, 1 - deaScamPy_tot$pval)
roc_resultScamPy
# Plot the ROC curve
#plot(roc_resultScamPy)
Summary ROC for all methods
<- ggroc(list(COTAN=roc_resultCOTAN, Seurat=roc_resultSeurat, Seurat_scTr = roc_resultSeurat_scTr, Seurat_Bimod = roc_resultSeurat_Bimod,
Three_Clusters_even Monocle=roc_resultMonocle, ScamPy=roc_resultScamPy))
<- Three_Clusters_even +
Three_Clusters_evenPL xlab("FPR") + ylab("TPR")+theme(legend.position="none")
Three_Clusters_even
3_Clusters_uneven
True vector
<-datasets_csv[datasets_csv$Group == "3_Clusters_uneven",]
subset.datasets_csv
<- NA
ground_truth_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
<- str_split(subset.datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- file.presence[,clusters]
file.presence.subset
#file.presence.subset <- as.matrix(file.presence.subset)
<- as.data.frame(matrix(nrow = nrow(file.presence.subset),
ground_truth ncol = ncol(file.presence.subset)))
rownames(ground_truth) <- rownames(file.presence.subset)
colnames(ground_truth) <- colnames(file.presence.subset)
== "Absent"] <- 0
ground_truth[file.presence.subset == "Present"] <- 1
ground_truth[file.presence.subset == "Uncertain"] <- 0.6
ground_truth[file.presence.subset
<- ground_truth
file.presence.subset
for (col in 1:ncol(ground_truth)) {
<- FALSE
ground_truth[,col] == 1 & rowMeans(file.presence.subset[,-col,drop = FALSE]) < 0.555 ,col] <- TRUE
ground_truth[file.presence.subset[,col]
}
$genes <- rownames(ground_truth)
ground_truth<- pivot_longer(ground_truth,
ground_truth cols = 1:(ncol(ground_truth)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
ground_truth$set_number <- ind
ground_truth<- rbind(ground_truth_tot, ground_truth)
ground_truth_tot
}<- ground_truth_tot[2:nrow(ground_truth_tot),]
ground_truth_tot
head(ground_truth_tot)
# A tibble: 6 × 5
genes clusters value data_set set_number
<chr> <chr> <lgl> <chr> <int>
1 Neil2 E15.0-428 FALSE 3_Clusters_uneven 1
2 Neil2 E13.5-434 FALSE 3_Clusters_uneven 1
3 Neil2 E15.0-510 FALSE 3_Clusters_uneven 1
4 Lamc1 E15.0-428 FALSE 3_Clusters_uneven 1
5 Lamc1 E13.5-434 FALSE 3_Clusters_uneven 1
6 Lamc1 E15.0-510 FALSE 3_Clusters_uneven 1
length(unique(ground_truth_tot$genes))
[1] 14695
sum(ground_truth_tot$value)
[1] 1577
ROC for COTAN
<- NA
onlyPositive.pVal.Cotan_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",
file.code $Collection[ind])
subset.datasets_csv<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset <- getClusterizationData(dataset,clName = "mergedClusters")[[2]]
deaCOTAN <- pValueFromDEA(deaCOTAN,
pvalCOTAN numCells = getNumCells(dataset),method = "none")
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization <- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- c(cl.names,
cl.names str_split(names(clusterization[clusterization == cl.val])[1],
pattern = "_",simplify = T)[1])
<- cl.names[!is.na(cl.names)]
cl.names
}colnames(deaCOTAN) <- cl.names
colnames(pvalCOTAN) <- cl.names
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- pvalCOTAN
onlyPositive.pVal.Cotan
for (cl in cl.names) {
print(cl)
#temp.DEA.CotanSign <- deaCOTAN[rownames(pvalCOTAN[pvalCOTAN[,cl] < 0.05,]) ,]
rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl] <- 1 #onlyPositive.pVal.Cotan[rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl]+1
onlyPositive.pVal.Cotan[
}
$genes <- rownames(onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan<- pivot_longer(onlyPositive.pVal.Cotan,
onlyPositive.pVal.Cotan values_to = "p_values",
cols = 1:(ncol(onlyPositive.pVal.Cotan)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
onlyPositive.pVal.Cotan$set_number <- ind
onlyPositive.pVal.Cotan<- rbind(onlyPositive.pVal.Cotan_tot, onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan_tot
}
[1] "E15.0-510"
[1] "E13.5-434"
[1] "E15.0-428"
[1] "E15.0-432"
[1] "E13.5-432"
[1] "E13.5-187"
[1] "E15.0-509"
[1] "E15.0-508"
[1] "E13.5-184"
<- onlyPositive.pVal.Cotan_tot[2:nrow(onlyPositive.pVal.Cotan_tot),]
onlyPositive.pVal.Cotan_tot
<- merge.data.frame(onlyPositive.pVal.Cotan_tot,
onlyPositive.pVal.Cotan_tot by = c("genes","clusters","data_set","set_number"),all.x = T,all.y = F)
ground_truth_tot,
#onlyPositive.pVal.Cotan_tot <- onlyPositive.pVal.Cotan_tot[order(onlyPositive.pVal.Cotan_tot$p_values,
# decreasing = F),]
# df <- as.data.frame(matrix(nrow = nrow(onlyPositive.pVal.Cotan_tot),ncol = 3))
# colnames(df) <- c("TPR","FPR","Method")
# df$Method <- "COTAN"
#
# Positive <- sum(onlyPositive.pVal.Cotan_tot$value)
# Negative <- sum(!onlyPositive.pVal.Cotan_tot$value)
#
# for (i in 1:nrow(onlyPositive.pVal.Cotan_tot)) {
# df[i,"TPR"] <- sum(onlyPositive.pVal.Cotan_tot[1:i,"value"])/Positive
# df[i,"FPR"] <- (i-sum(onlyPositive.pVal.Cotan_tot[1:i,"value"]))/Negative
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(onlyPositive.pVal.Cotan_tot$value)
onlyPositive.pVal.Cotan_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(onlyPositive.pVal.Cotan_tot$value, 1 - onlyPositive.pVal.Cotan_tot$p_values)
roc_resultCOTAN
# Plot the ROC curve
#plot(roc_resultCOTAN)
ROC for Seurat
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val)
roc_resultSeurat
# Plot the ROC curve
#plot(roc_resultSeurat)
ROC for Seurat scTransform
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_ScTransform_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_scTr
ROC for Seurat Bimod
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_Bimod_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_Bimod
ROC for Monocle
<- NA
deaMonocle_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"Monocle_DEA_genes.csv")),row.names = 1)
deaMonocle
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cell_group == cl.val,"cell_group"] <- cl.name
deaMonocle[deaMonocle
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "cell_group",
replacement = "clusters")
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "gene_id",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaMonocle$set_number <- ind
deaMonocle<- as.data.frame(deaMonocle)
deaMonocle <- rbind(deaMonocle_tot, deaMonocle)
deaMonocle_tot
}<- deaMonocle_tot[2:nrow(deaMonocle_tot),]
deaMonocle_tot
<- merge.data.frame(deaMonocle_tot,
deaMonocle_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaMonocle_tot$value)
deaMonocle_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaMonocle_tot$value, 1 - deaMonocle_tot$marker_test_p_value)
roc_resultMonocle
# Plot the ROC curve
#plot(roc_resultMonocle)
ROC from ScamPy
<- NA
deaScamPy_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"ScamPy_DEA_genes.csv")),
deaScamPy row.names = 1)
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$clusters == paste0("cl",cl.val),"clusters"] <- cl.name
deaScamPy[deaScamPy
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
$data_set <- subset.datasets_csv[ind,1]
deaScamPy$set_number <- ind
deaScamPy<- rbind(deaScamPy_tot, deaScamPy)
deaScamPy_tot
}<- deaScamPy_tot[2:nrow(deaScamPy_tot),]
deaScamPy_tot
<- merge.data.frame(deaScamPy_tot,
deaScamPy_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaScamPy_tot$value)
deaScamPy_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaScamPy_tot$value, 1 - deaScamPy_tot$pval)
roc_resultScamPy
# Plot the ROC curve
#plot(roc_resultScamPy)
Summary ROC for all methods
<- ggroc(list(COTAN=roc_resultCOTAN, Seurat=roc_resultSeurat, Seurat_scTr = roc_resultSeurat_scTr, Seurat_Bimod = roc_resultSeurat_Bimod,
Three_Clusters_uneven Monocle=roc_resultMonocle, ScamPy=roc_resultScamPy))
<- Three_Clusters_uneven + xlab("FPR") + ylab("TPR")+theme(legend.position="none")
Three_Clusters_unevenPL
Three_Clusters_uneven
5 clusters
5_Clusters_uneven
True vector
<-datasets_csv[datasets_csv$Group == "5_Clusters_uneven",]
subset.datasets_csv
<- NA
ground_truth_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
<- str_split(subset.datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- file.presence[,clusters]
file.presence.subset
#file.presence.subset <- as.matrix(file.presence.subset)
<- as.data.frame(matrix(nrow = nrow(file.presence.subset),
ground_truth ncol = ncol(file.presence.subset)))
rownames(ground_truth) <- rownames(file.presence.subset)
colnames(ground_truth) <- colnames(file.presence.subset)
== "Absent"] <- 0
ground_truth[file.presence.subset == "Present"] <- 1
ground_truth[file.presence.subset == "Uncertain"] <- 0.6
ground_truth[file.presence.subset
<- ground_truth
file.presence.subset
for (col in 1:ncol(ground_truth)) {
<- FALSE
ground_truth[,col] == 1 & rowMeans(file.presence.subset[,-col,drop = FALSE]) < 0.555 ,col] <- TRUE
ground_truth[file.presence.subset[,col]
}
$genes <- rownames(ground_truth)
ground_truth<- pivot_longer(ground_truth,
ground_truth cols = 1:(ncol(ground_truth)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
ground_truth$set_number <- ind
ground_truth<- rbind(ground_truth_tot, ground_truth)
ground_truth_tot
}<- ground_truth_tot[2:nrow(ground_truth_tot),]
ground_truth_tot
head(ground_truth_tot)
# A tibble: 6 × 5
genes clusters value data_set set_number
<chr> <chr> <lgl> <chr> <int>
1 Neil2 E13.5-510 FALSE 5_Clusters_uneven 1
2 Neil2 E15.0-437 FALSE 5_Clusters_uneven 1
3 Neil2 E15.0-510 FALSE 5_Clusters_uneven 1
4 Neil2 E13.5-432 FALSE 5_Clusters_uneven 1
5 Neil2 E13.5-437 FALSE 5_Clusters_uneven 1
6 Lamc1 E13.5-510 FALSE 5_Clusters_uneven 1
length(unique(ground_truth_tot$genes))
[1] 14695
sum(ground_truth_tot$value)
[1] 1862
ROC for COTAN
<- NA
onlyPositive.pVal.Cotan_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",
file.code $Collection[ind])
subset.datasets_csv<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset <- getClusterizationData(dataset,clName = "mergedClusters")[[2]]
deaCOTAN <- pValueFromDEA(deaCOTAN,
pvalCOTAN numCells = getNumCells(dataset),method = "none")
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization <- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- c(cl.names,
cl.names str_split(names(clusterization[clusterization == cl.val])[1],
pattern = "_",simplify = T)[1])
<- cl.names[!is.na(cl.names)]
cl.names
}colnames(deaCOTAN) <- cl.names
colnames(pvalCOTAN) <- cl.names
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
<- pvalCOTAN
onlyPositive.pVal.Cotan
for (cl in cl.names) {
print(cl)
#temp.DEA.CotanSign <- deaCOTAN[rownames(pvalCOTAN[pvalCOTAN[,cl] < 0.05,]) ,]
rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl] <- 1 #onlyPositive.pVal.Cotan[rownames(deaCOTAN[deaCOTAN[,cl] < 0,]),cl]+1
onlyPositive.pVal.Cotan[
}
$genes <- rownames(onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan<- pivot_longer(onlyPositive.pVal.Cotan,
onlyPositive.pVal.Cotan values_to = "p_values",
cols = 1:(ncol(onlyPositive.pVal.Cotan)-1),
names_to = "clusters")
$data_set <- subset.datasets_csv[ind,1]
onlyPositive.pVal.Cotan$set_number <- ind
onlyPositive.pVal.Cotan<- rbind(onlyPositive.pVal.Cotan_tot, onlyPositive.pVal.Cotan)
onlyPositive.pVal.Cotan_tot
}
[1] "E13.5-432"
[1] "E15.0-510"
[1] "E13.5-437"
[1] "E15.0-437"
[1] "E13.5-510"
[1] "E13.5-434"
[1] "E15.0-428"
[1] "E13.5-184"
[1] "E15.0-434"
[1] "E17.5-505"
[1] "E15.0-432"
[1] "E13.5-432"
[1] "E15.0-509"
[1] "E15.0-508"
[1] "E13.5-187"
<- onlyPositive.pVal.Cotan_tot[2:nrow(onlyPositive.pVal.Cotan_tot),]
onlyPositive.pVal.Cotan_tot
<- merge.data.frame(onlyPositive.pVal.Cotan_tot,
onlyPositive.pVal.Cotan_tot by = c("genes","clusters","data_set","set_number"),all.x = T,all.y = F)
ground_truth_tot,
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(onlyPositive.pVal.Cotan_tot$value)
onlyPositive.pVal.Cotan_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(onlyPositive.pVal.Cotan_tot$value, 1 - onlyPositive.pVal.Cotan_tot$p_values)
roc_resultCOTAN
# Plot the ROC curve
#plot(roc_resultCOTAN)
ROC for Seurat
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val)
roc_resultSeurat
# Plot the ROC curve
#plot(roc_resultSeurat)
ROC for Seurat scTransform
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_ScTransform_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_scTr
ROC for Seurat Bimod
<- NA
deaSeurat_tot
for (ind in 1:dim(subset.datasets_csv)[1]) {
#print(ind)
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code
<- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization
<- read.csv(file.path(dirOut,paste0(file.code,"Seurat_DEA_Bimod_genes.csv")), row.names = 1)
deaSeurat
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cluster == cl.val,]$cluster <- cl.name
deaSeurat[deaSeurat
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "cluster",
replacement = "clusters")
colnames(deaSeurat) <- str_replace(colnames(deaSeurat),
pattern = "gene",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaSeurat$set_number <- ind
deaSeurat<- rbind(deaSeurat_tot, deaSeurat)
deaSeurat_tot
}<- deaSeurat_tot[2:nrow(deaSeurat_tot),]
deaSeurat_tot
<- merge.data.frame(deaSeurat_tot,
deaSeurat_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaSeurat_tot$value)
deaSeurat_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaSeurat_tot$value, 1 - deaSeurat_tot$p_val) roc_resultSeurat_Bimod
ROC for Monocle
<- NA
deaMonocle_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"Monocle_DEA_genes.csv")),row.names = 1)
deaMonocle
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$cell_group == cl.val,"cell_group"] <- cl.name
deaMonocle[deaMonocle
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "cell_group",
replacement = "clusters")
colnames(deaMonocle) <- str_replace(colnames(deaMonocle),
pattern = "gene_id",
replacement = "genes")
$data_set <- subset.datasets_csv[ind,1]
deaMonocle$set_number <- ind
deaMonocle<- as.data.frame(deaMonocle)
deaMonocle <- rbind(deaMonocle_tot, deaMonocle)
deaMonocle_tot
}<- deaMonocle_tot[2:nrow(deaMonocle_tot),]
deaMonocle_tot
<- merge.data.frame(deaMonocle_tot,
deaMonocle_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaMonocle_tot$value)
deaMonocle_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaMonocle_tot$value, 1 - deaMonocle_tot$marker_test_p_value)
roc_resultMonocle
# Plot the ROC curve
#plot(roc_resultMonocle)
ROC from ScamPy
<- NA
deaScamPy_tot for (ind in 1:dim(subset.datasets_csv)[1]) {
<- paste0(subset.datasets_csv$Group[ind],"_",subset.datasets_csv$Collection[ind])
file.code <- readRDS(file = file.path(dataSetDir,paste0(file.code,".RDS")))
dataset
<- getClusterizationData(dataset, clName = "mergedClusters")[[1]]
clusterization #print(file.code)
<- read.csv(file.path(dirOut,paste0(file.code,"ScamPy_DEA_genes.csv")),
deaScamPy row.names = 1)
<- NA
cl.names
for (cl.val in unique(clusterization)) {
#print(cl.val)
<- str_split(names(clusterization[clusterization == cl.val])[1],
cl.name pattern = "_",simplify = T)[1]
<- c(cl.names,cl.name)
cl.names <- cl.names[!is.na(cl.names)]
cl.names
$clusters == paste0("cl",cl.val),"clusters"] <- cl.name
deaScamPy[deaScamPy
}
<- str_split(datasets_csv$Collection[ind],pattern = "_[+]_",simplify = T)
clusters
$data_set <- subset.datasets_csv[ind,1]
deaScamPy$set_number <- ind
deaScamPy<- rbind(deaScamPy_tot, deaScamPy)
deaScamPy_tot
}<- deaScamPy_tot[2:nrow(deaScamPy_tot),]
deaScamPy_tot
<- merge.data.frame(deaScamPy_tot,
deaScamPy_tot
ground_truth_tot,by = c("genes","clusters","data_set","set_number"),
all.x = T,all.y = F)
# Convert TRUE/FALSE to 1/0
$value <- as.numeric(deaScamPy_tot$value)
deaScamPy_tot
# Compute the ROC curve - note that we invert the p-values with 1 - p_values
<- roc(deaScamPy_tot$value, 1 - deaScamPy_tot$pval)
roc_resultScamPy
# Plot the ROC curve
#plot(roc_resultScamPy)
Summary ROC for all methods
<- ggroc(list(COTAN=roc_resultCOTAN, Seurat=roc_resultSeurat, Seurat_scTr = roc_resultSeurat_scTr, Seurat_Bimod = roc_resultSeurat_Bimod,
Five_Clusters Monocle=roc_resultMonocle, ScamPy=roc_resultScamPy))
<- Five_Clusters + xlab("FPR") + ylab("TPR")
Five_Clusters_unevenPL
Five_Clusters_unevenPL
Global Summary
2 Clusters
ggarrange(TwoClusters_even_nearPL,Two_Clusters_even_mediumPL, Two_Clusters_even_farPL,Two_Clusters_uneven_nearPL, Two_Clusters_uneven_mediumPL,Two_Clusters_uneven_farPL,
labels = c("Even_Near", "Even_Medium", "Even_Far", "Uneven_Near","Uneven_Medium","Uneven_Far"),
ncol = 3, nrow = 2, common.legend = T, legend = "bottom")
3 and 5 Clusters
ggarrange(Three_Clusters_evenPL,Three_Clusters_unevenPL, NULL, Five_Clusters_unevenPL,
labels = c("3_Even", "3_Uneven", "", "5_Uneven"),
ncol = 2, nrow = 2, common.legend = T, legend = "bottom")