-
Notifications
You must be signed in to change notification settings - Fork 3
/
Figure3_Analyze_Markers.R
161 lines (138 loc) · 6.53 KB
/
Figure3_Analyze_Markers.R
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
source("~/MAGIC/CTP_Functions.R")
source("~/MAGIC/Colour_Scheme.R")
require("scater")
library(methods)
# Matching files
file_pairs <- rbind(
c("Bladder_10X_filtered.rds", "Bladder_FACS_filtered.rds"),
c("Kidney_10X_filtered.rds", "Kidney_FACS_filtered.rds"),
c("Lung_10X_filtered.rds", "Lung_FACS_filtered.rds"),
c("Mammary_10X_filtered.rds", "Mammary_FACS_filtered.rds"),
c("Marrow_10X_filtered.rds", "Marrow_FACS_filtered.rds"),
c("Tongue_10X_filtered.rds", "Tongue_FACS_filtered.rds"),
c("Muscle_10X_filtered.rds", "Muscle_FACS_filtered.rds"))
rownames(file_pairs) <- c("Bladder", "Kidney", "Lung", "Mammary", "Marrow", "Tongue", "Muscle")
# Command line argument: AUC threshold
#args <- commandArgs(trailingOnly=TRUE)
#auc_threshold <- as.numeric(args[1])
p_val_cors <- matrix(0, nrow=nrow(file_pairs), ncol=length(Method_Legend))
rownames(p_val_cors) <- rownames(file_pairs);
cross_dataset_outfile <- paste("TM", "Data_Consistency.rds", sep="_");
cross_method_outfile <- paste("TM", "Method_Consistency.rds", sep="_");
allDataOUT <- vector();
allMethOUT <- vector();
for (i in 1:nrow(file_pairs)) {
drop <- readRDS(file_pairs[i,1])
facs <- readRDS(file_pairs[i,2])
# Only cell-types that exist in both datasets
types <- intersect( as.character(unique(drop$cell_type1)), as.character(unique(facs$cell_type1)) )
drop <- drop[, drop$cell_type1 %in% types]
facs <- facs[, facs$cell_type1 %in% types]
# Exclude any permuted genes if present
# if ("Permuted" %in% colnames(rowData(drop))) {
# perm <- c( rownames(drop)[rowData(drop)$Permuted], rownames(facs)[rowData(facs)$Permuted] )
# drop <- drop[!rownames(drop) %in% perm, ]
# facs <- facs[!rownames(facs) %in% perm, ]
# }
genes <- intersect(rownames(drop), rownames(facs));
drop <- drop[match(genes, rownames(drop)),]
facs <- facs[match(genes, rownames(facs)),]
# Exclude low detection genes
exclude <- Matrix::rowSums(drop@assays[["counts"]] > 0) < ncol(drop@assays[["counts"]])*0.01 |
Matrix::rowSums(facs@assays[["counts"]] > 0) < ncol(facs@assays[["counts"]])*0.01
drop <- drop[!exclude,]
facs <- facs[!exclude,]
# Setup output
methods <- names(assays(drop));
methods_nice_names <- Method_Colours[match(methods, Method_Colours[,1]),2]
methods_nice_names <- names(Method_Legend)[match(methods_nice_names, Method_Legend)]
methods <- methods[!duplicated(methods_nice_names)]
methods_nice_names <- methods_nice_names[!duplicated(methods_nice_names)]
xDataOUT <- vector() #matrix(nrow=length(methods), ncol=5)
#rownames(xDataOUT) <- methods_nice_names
#colnames(xDataOUT) <- c("Ndrop", "Nfacs", "Nboth", "Nconsist", "perc_consist")
Method_Markers_Drop <- list();
Method_Markers_Facs <- list();
tissue <- rownames(file_pairs)[i]
#tissue_outfile <- paste(tissue, "_", auc_threshold, "_consistency_results.rds", sep="");
tissue_outfile <- paste(tissue, "_consistency_results.rds", sep="");
# Identify Markers
for (m in 1:length(methods)) {
meth <- methods[m];
# Get markers
mark_drop <- complex_markers(as.matrix(assays(drop)[[meth]]), drop$cell_type1, n_max=1)
saveRDS(mark_drop, paste(tissue, meth, "10X","markers.rds", sep="_"))
mark_facs <- complex_markers(as.matrix(assays(facs)[[meth]]), factor(facs$cell_type1), n_max=1)
saveRDS(mark_facs, paste(tissue, meth, "FACS","markers.rds", sep="_"))
# p-value correlation
tmp_facs <- mark_facs[rownames(mark_facs) %in% rownames(mark_drop),]
tmp_drop <- mark_drop[match(rownames(tmp_facs), rownames(mark_drop)),]
c <- cor(tmp_facs$p.value, tmp_drop$p.value, method="spearman")
p_val_cors[i,m] <- c
colnames(p_val_cors) <- remap_names(methods)
dataout <- vector();
for (auc_threshold in c(0, 0.5, 0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.85, 0.9)){
# apply thresholds
mark_drop_gd <- mark_drop[mark_drop$AUC > auc_threshold & mark_drop$q.value < 0.05,]
mark_facs_gd <- mark_facs[mark_facs$AUC > auc_threshold & mark_facs$q.value < 0.05,]
if (auc_threshold == 0) {
Method_Markers_Drop[[meth]] <- cbind(
rownames(mark_drop_gd), mark_drop_gd$Group)
Method_Markers_Facs[[meth]] <- cbind(
rownames(mark_facs_gd), mark_facs_gd$Group)
}
Ndrop <- nrow(mark_drop_gd)
Nfacs <- nrow(mark_facs_gd)
# Consistent
mark_drop_gd <- mark_drop_gd[rownames(mark_drop_gd) %in% rownames(mark_facs_gd), ]
mark_facs_gd <- mark_facs_gd[match(rownames(mark_drop_gd), rownames(mark_facs_gd)), ]
Nboth <- nrow(mark_drop_gd);
mark_drop_gd$Group <- as.character(mark_drop_gd$Group)
mark_facs_gd$Group <- as.character(mark_facs_gd$Group)
perc_agreement <- sum(mark_drop_gd$Group == mark_facs_gd$Group)/nrow(mark_drop_gd);
dataout <- rbind(dataout, c(auc_threshold, Ndrop, Nfacs, Nboth, sum(mark_drop_gd$Group == mark_facs_gd$Group), perc_agreement))
}
dataout <- data.frame(tissue=as.character(tissue),method=as.character(meth), dataout, stringsAsFactors=FALSE)
colnames(dataout) <- c("tissue","method", "auc", "ndrop", "nfacs", "nboth", "nagree", "percagree")
xDataOUT <- rbind(xDataOUT, dataout);
}
# Cross method consistency
DropxMethOUT <- matrix(nrow=length(methods), ncol=length(methods))
rownames(DropxMethOUT) <- methods_nice_names
colnames(DropxMethOUT) <- methods_nice_names
FACSxMethOUT <- matrix(nrow=length(methods), ncol=length(methods))
rownames(FACSxMethOUT) <- methods_nice_names
colnames(FACSxMethOUT) <- methods_nice_names
for (r in 1:length(methods)) {
for (j in 1:length(methods)) {
meth1 <- methods[r];
meth2 <- methods[j];
markers1 <- Method_Markers_Drop[[meth1]]
markers2 <- Method_Markers_Drop[[meth2]]
# Markers by both methods
markers1 <- markers1[markers1[,1] %in% markers2[,1],]
markers2 <- markers2[match(markers1[,1], markers2[,1]),]
# % same group
perc_agree1 <- sum( markers1[,2] == markers2[,2])/nrow(markers1);
DropxMethOUT[r,j] <- perc_agree1
markers1 <- Method_Markers_Facs[[meth1]]
markers2 <- Method_Markers_Facs[[meth2]]
# Markers by both methods
markers1 <- markers1[markers1[,1] %in% markers2[,1],]
markers2 <- markers2[match(markers1[,1], markers2[,1]),]
# % same group
perc_agree2 <- sum( markers1[,2] == markers2[,2])/nrow(markers1);
FACSxMethOUT[r,j] <- perc_agree2
}}
# saveRDS(list(crossData=xDataOUT, DropxMethOUT=DropxMethOUT, FACSxMethOUT=FACSxMethOUT), file=tissue_outfile )
allDataOUT <- rbind(allDataOUT, xDataOUT)
if (is.null(dim(allMethOUT))) {
allMethOUT <- DropxMethOUT+FACSxMethOUT
} else {
allMethOUT <- allMethOUT+DropxMethOUT+FACSxMethOUT
}
}
allMethOUT <- allMethOUT/(nrow(file_pairs)*2)
#saveRDS(allMethOUT, file=cross_method_outfile)
#saveRDS(allDataOUT, file=cross_dataset_outfile)
#saveRDS(p_val_cors, file="TM_mark_cors.rds")