-
Notifications
You must be signed in to change notification settings - Fork 0
/
Sensitivity_Test.R
150 lines (130 loc) · 5.49 KB
/
Sensitivity_Test.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
library(tidyverse)
library(SingleCellExperiment)
library(scp)
library(ggplot2)
mod_pep <- read_tsv("C10SVEC_singlecells_Peptide_intensities.tsv")
peptide1 <- mod_pep %>%
filter(grepl("MOUSE", Protein)) %>%
dplyr::select(`Modified.Sequence` | contains("_")) %>%
column_to_rownames(var="Modified.Sequence")
###################
mod_pep <- read_tsv("C10Treated_Peptide_intensities.tsv")
meta1 <- read.csv("C10Treated_NormCounts.csv", check.names = T) %>%
select(contains("Ctrl")| contains("Trtd")) %>%
colnames(.)
peptide2 <- mod_pep %>%
filter(grepl("MOUSE", Protein)) %>%
dplyr::select(`Modified.Sequence` | contains("Ctrl")| contains("Trtd")) %>%
column_to_rownames(var="Modified.Sequence") %>%
dplyr::select((meta1))
######################
mod_pep <- read_tsv("combined_modified_peptide_MaxLFQ.tsv")
meta2 <- read.csv("annotations.csv")
peptide3 <- mod_pep %>%
filter(grepl("HUMAN", Protein)) %>%
dplyr::select(`Modified.Sequence` | contains("_")) %>%
column_to_rownames(var="Modified.Sequence") %>%
select(contains(meta2$SampleID))
peptide1 <- peptide1 != 0
peptide2 <- peptide2 != 0
peptide3 <- peptide3 != 0
sampleGroup1 <- data.frame(SampleID = colnames(peptide1)) %>%
mutate(Group = case_when(grepl("_C10_", SampleID)~ "C10",
grepl("_SVEC_", SampleID)~ "SVEC"))
sampleGroup2 <- data.frame(SampleID = colnames(peptide2)) %>%
mutate(Group = "CDK1_Experiment")
sampleGroup3 <- data.frame(SampleID = colnames(peptide3)) %>%
mutate(Group = "Pancreatic_Cells")
sce <- SingleCellExperiment(assays = peptide1)
sce2 <- SingleCellExperiment(assays = peptide2)
sce3 <- SingleCellExperiment(assays = peptide3)
sim <- QFeatures(experiments = List(Assay1 = sce,
Assay2 = sce2,
Assay3 = sce3))
pep_Missing1 <- reportMissingValues(sim, "Assay1", by = sampleGroup1$Group)
pep_Missing2 <- reportMissingValues(sim, "Assay2", by = sampleGroup2$Group)
pep_Missing3 <- reportMissingValues(sim, "Assay3", by = sampleGroup3$Group)
csc1 <- cumulativeSensitivityCurve(
sim, "Assay1", by = sampleGroup1$Group,
)
csc2 <- cumulativeSensitivityCurve(
sim, "Assay2", by = sampleGroup2$Group,
)
csc3 <- cumulativeSensitivityCurve(
sim, "Assay3", by = sampleGroup3$Group,
)
ji1 <- jaccardIndex(sim, "Assay1", by = sampleGroup1$Group)
ji2 <- jaccardIndex(sim, "Assay2", by = sampleGroup2$Group)
ji3 <- jaccardIndex(sim, "Assay3", by = sampleGroup3$Group)
csc <- do.call("rbind", list(csc1,csc2,csc3))
csc %>%
ggplot()+
aes(x = SampleSize, y = Sensitivity, colour = by) +
geom_point()+
geom_hline(yintercept = pep_Missing1$LocalSensitivityMean,
linetype = "dashed",
color = c("#F8766D","#C77CFF"))+
geom_hline(yintercept = pep_Missing3$LocalSensitivityMean,
linetype = "dashed",
color = "#00BFC4")+
geom_hline(yintercept = pep_Missing2$LocalSensitivityMean,
linetype = "dashed",
color = "#7CAE00")+
scale_y_continuous(limits = c(0,26000))+
theme_bw(base_size = 20) +
theme(#legend.position = "none",
panel.background = element_rect(fill= 'white'),
axis.text.y=element_text(color = 'black', size = 20),
axis.text.x=element_text(angle = 0, vjust = 0.5, hjust= 0.5, color='black',size = 20),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_blank(),
axis.line = element_line())+
ylab("Sensitivity (Peptides ,n)")+
xlab("Number of Cells")
ji <- do.call("rbind", list(ji1,ji2,ji3))
ji %>%
ggplot()+
aes(y = jaccard, x = by, fill = by) +
geom_violin()+
scale_y_continuous(limits = c(0,1))+
theme_bw(base_size = 15) +
theme(legend.position = "none",
panel.background = element_rect(fill= 'white'),
axis.text.y=element_text(color = 'black', size = 12),
axis.text.x=element_text(angle = 0, vjust = 0.5, hjust= 0.5, color='black',size = 15),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_blank(),
axis.line = element_line())+
ylab("Jaccard Index")
data_completedness <- do.call("rbind", list(pep_Missing1,pep_Missing2,pep_Missing3)) %>%
mutate(Completeness = Completeness * 100) %>%
mutate(Type = rownames(.))
data_completedness %>%
ggplot()+
aes(x = Completeness, y = LocalSensitivityMean, color = Type,
group = Type,
size = NumberCells)+
geom_point(size =10,
shape = 5)+
geom_point()+
theme_bw()+
scale_x_continuous(limits = c(0,90))+
scale_y_continuous(limits = c(0,20000))+
geom_pointrange(aes(ymin=LocalSensitivityMean-LocalSensitivitySd,
ymax=LocalSensitivityMean+LocalSensitivitySd),
position=position_dodge(0.05))+
theme_bw(base_size = 20) +
theme(panel.background = element_rect(fill= 'white'),
axis.text.y=element_text(color = 'black', size = 20),
axis.text.x=element_text(angle = 0, vjust = 0.5, hjust= 0.5, color='black',size = 20),
panel.grid.minor = element_blank(),
panel.grid.major = element_blank(),
panel.border = element_blank(),
axis.line = element_line(),
text=element_text(family="Helvetica"))+
ylab("Local Sensitivty")+
xlab("Data completeness (%)")+
scale_size(range = c(0.27,1.2),
breaks = c(31,39,62,128))