-
Notifications
You must be signed in to change notification settings - Fork 0
Expand file tree
/
Copy pathfocusGenes.R
More file actions
209 lines (185 loc) · 11 KB
/
Copy pathfocusGenes.R
File metadata and controls
209 lines (185 loc) · 11 KB
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
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
FocusGenesIntegrationPlot = function(fileGenes, folderStoredPlots,
integrSiteObject, includeOnly,
aes_y_char, aes_shape_char, aes_color_char,# aes_size_char,
selectedCovariates = character(),
selectedFactorCovariates = character(),
covariatesInLabels = "time",
ylab = "Passage",
formatPlot = c(".png", ".svg", ".pdf")[1],
showLegend = FALSE,
transparency = 1,
jitter = 0,
colors = list(),
sizeInchesPlots = c(9, 2)) {
if(class(integrSiteObject) != "integrSiteObject") stop("error")
iSO = integrSiteObject
it = iSO$design$filesIntegrations
cv = iSO$design$covariatesIntegrations
if(length(it) == 0L | length(cv) == 0L) stop("error")
if(length(iSO$data$combined) != 0L) stop("error") # not yet implemented
if(!missing(includeOnly)) {
iSO$design$filesIntegrations = it = it[includeOnly]
iSO$design$covariatesIntegrations = cv = cv[includeOnly,]
iSO$design$id = 1:length(it)
}
#selectedCovariates = unique(c(aes_y_char, aes_shape_char, aes_color_char, selectedCovariates))#, aes_size_char)
if(!all(selectedCovariates %in% c("chr", "strand", "count", names(cv)))) stop("error")
#if(!all(selectedFactorCovariates %in% selectedCovariates)) stop("error")
fg = data.table :: fread(fileGenes, sep = "\t", header = T)
iSO = setGroups(iSO, group1 = 1:length(it), group2 = integer(0L))
iSO = setInfoSelectedCovariates(iSO, selectedCovariates = selectedCovariates, selCovAreFactors = selectedFactorCovariates)
iSO = loadIntegrations(iSO, excludeIntegrations = list())
iSO = loadLocationsGenes(iSO)
iSO = mergeData(iSO)
#iSO$design$selectedCovariates = selectedCovariates
#iSO$design$selCovAreFactors = sapply(selectedCovariates, function(x) x %in% selectedFactorCovariates)
#class(iSO)
#iSO = setFactorsCombinedData(iSO)
#iSO = mergeData(iSO, selectedCovariates[which(!(selectedCovariates %in% c("chr", "start", "end", "strand", "group")))])
#iSO = setFactorsCombinedData(iSO, variables = c(c("chr", "strand", "group"), selectedCovariates[which(selectedCovariates %in% selectedFactorCovariates)]))
selectedFactorCovariates = c("strand", selectedCovariates[selectedFactorCovariates])
if(length(selectedFactorCovariates) > 0L) for(i in 1:length(selectedFactorCovariates)) iSO$data$combined[[selectedFactorCovariates[i]]] = as.factor(iSO$data$combined[[selectedFactorCovariates[i]]])
iSO$data$combined$relativeAbundance = iSO$data$combined$count / iSO$data$combined$totalCount
aes_size_char = "relativeAbundance"
textSize = 16
#colorCell = list(`BM MSC` = "#0054ff", `Ad MSC` = "#008b0b", `HSPC` = "#090909")
i = 1L
for(i in 1:nrow(fg)) {
da = iSO$data$combined[chr == fg$chr[i],]
#da = iSO$data$combined[start >= fg$start[i] & end < fg$end[i],]
da = da[start >= fg$start[i] & end < fg$end[i],]
p = ggplot2 :: ggplot(data = da)
#p = p + geom_point(data = da, aes(x = pos, y = timePointMonths, shape = strand, color = cellType, size = relAbund ))
p = p + ggplot2 :: geom_point(ggplot2 :: aes(x = start,
y = .data[[aes_y_char]],
shape = .data[[aes_shape_char]],
color = .data[[aes_color_char]],
size = .data[[aes_size_char]]),
position = ggplot2 :: position_jitter(width = 0, height = jitter),
alpha = transparency)
if(length(colors)) p = p + ggplot2 :: scale_color_manual(values = as.list(colors))
p = p + ggplot2 :: scale_shape_manual(values = c("-"="<","+"=">"))
#scale_size_continuous(breaks = c(0.005,0.01,0.05,0.1,0.2,0.3))+
p = p + ggplot2 :: scale_size_continuous(range = c(4, 10))
# geom_segment(aes(x=startS,y=0.5,xend=endS,yend=0.5),size=2)+
p = p + ggplot2 :: scale_y_discrete(labels = unique(cv[, covariatesInLabels]))
p = p + ggplot2 :: xlim(fg$start[i],fg$end[i])
p = p + ggplot2 :: theme(text= ggplot2 :: element_text(family="sans"),
axis.title.x = ggplot2 :: element_text(size =textSize),
axis.title.y = ggplot2 :: element_text(size =textSize),
panel.background = ggplot2 :: element_rect(fill = NA),
panel.grid.major.y = ggplot2 :: element_line(colour = "grey90"),
panel.grid.minor.y = ggplot2 :: element_line(colour =NA),
axis.text.x = ggplot2 :: element_text(size=textSize),
axis.text.y = ggplot2 :: element_text(size=textSize))#,
#legend.position = "bottom")
# legend.position = "none",
# panel.grid.major.x = element_line(colour = "grey90"),
# panel.grid.minor.x = element_line(colour = "grey90"),
#axis.ticks.x=element_text(size=14),
p = p + ggplot2 :: xlab(fg$chr[i])
p = p + ggplot2 :: ylab(ylab)
p = p + ggplot2 :: guides(shape = ggplot2 :: guide_legend(override.aes=list(size = 3)),
size = ggplot2 :: guide_legend(override.aes=list(shape = ">"))#,
#color = ggplot2 :: guide_legend(override.aes=list(alpha = transparency, size = 3))
)
ggplot2 :: ggsave(filename = paste0(folderStoredPlots,
c("integrations_", "legend_")[1+showLegend],
paste(fg$geneName[i],fg$chr[i],fg$start[i],fg$end[i],"2", formatPlot ,sep="_",collapse = "")),
#plot = p, width = 9,height = c(2,7)[1+showLegend],units = "in")
plot = p, width = sizeInchesPlots[1], height = sizeInchesPlots[2], units = "in")
}
return(invisible(NULL))
}
FocusGenesTrackingPlot = function(fileGenes, folderStoredPlots, integrSiteObject, includeOnly,
formatPlot = c(".png", ".svg", ".pdf")[1],
covariatesInLabels = c("time"),
compressUniqueCovariates = TRUE,
sizeInchesPlots = c(6, 4),
colors = list()) {
if(class(integrSiteObject) != "integrSiteObject") stop("error")
iSO = integrSiteObject
it = iSO$design$filesIntegrations
cv = iSO$design$covariatesIntegrations
if(length(it) == 0L | length(cv) == 0L) stop("error")
if(length(iSO$data$combined) != 0L) stop("error") # not yet implemented
if(!missing(includeOnly)) {
iSO$design$filesIntegrations = it = it[includeOnly]
iSO$design$covariatesIntegrations = cv = cv[includeOnly,]
iSO$design$id = 1:length(it)
}
fg = data.table :: fread(fileGenes, sep = "\t", header = T)
if(compressUniqueCovariates) {
dp = duplicated(cv)
cv = unique(cv)
if(any(dp)) {
ss = which(!dp)
es = which(c(diff(dp), -1L) == -1L)
cv$series = purrr::map2(ss, es, ~ .x : .y)
}
else {
cv$series = 1:nrow(cv)
}
}
else {
cv$series = 1:nrow(cv)
}
iSO = setGroups(iSO, group1 = 1:length(it), group2 = integer(0L))
iSO = loadIntegrations(iSO, excludeIntegrations = list())
iSO = loadLocationsGenes(iSO)
iSO = mergeData(iSO)
iSO$data$combined$relativeAbundance = iSO$data$combined$count / iSO$data$combined$totalCount
ic = iSO$data$combined
textSize = 13
pchSize = 2
colorCell = list(`BM MSC` = "#0054ff", `Ad MSC` = "#008b0b", `HSPC` = "#090909")
minAlpha = .2
maxAlpha = .8
minTime = min(cv$time)
maxTime = max(cv$time)
i = 1L
for(i in 1:nrow(fg)) {
totAb = rep(0, nrow(cv))
relAb = replicate(n = nrow(cv), expr = numeric(), simplify = "list")
for(j in 1:nrow(cv)) {
#integrInGene_i = ic[which(ic$series == j & ic$chr == fg$chr[i] & ic$start >= fg$start[i] & ic$start < fg$end[i]),]
integrInGene_i = ic[which((ic$series %in% cv$series[[j]]) & (ic$chr == fg$chr[i]) & (ic$start >= fg$start[i]) & (ic$start < fg$end[i])),]
if(nrow(integrInGene_i) > 0L) {
relAb[[j]] = integrInGene_i$relativeAbundance
totAb[j] = sum(relAb[[j]])
}
}
lenAb = sapply(relAb, length)
selAb = rep(1:nrow(cv), times = lenAb)
#alpAb = unlist(lapply(split(totAb, cv$type), function(x) (x - min(x)) / (max(x) - min(x))))[order(cv$type)]
alpAb = sapply(cv$time, function(x) (x - minTime) / (maxTime - minTime))
#d1 = data.frame(series = 1:nrow(cv), total = totAb, type = cv$type, labels = paste(cv$type, cv$time), alpha = alpAb)
d1 = data.frame(series = 1:nrow(cv), total = totAb, type = cv$type,
labels = apply(cv[,covariatesInLabels], 1, function(x) paste(x, collapse = " ")),
alpha = alpAb)
d2 = data.frame(series = selAb, relative = unlist(relAb), type = cv$type[selAb])
p = ggplot2 :: ggplot(data = d1)
p = p + ggplot2 :: geom_col(ggplot2 :: aes(x = series, y = total, fill = type, alpha = alpha))
if(length(colors) > 0) p = p + ggplot2 :: scale_fill_manual(values = as.list(colors))
p = p + ggplot2 :: scale_alpha_continuous(range = c(minAlpha, maxAlpha))
p = p + ggplot2 :: scale_x_continuous(breaks = 1:nrow(cv), labels = d1$labels)
p = p + ggplot2 :: scale_y_continuous(labels = function(x) x * 100)
p = p + ggnewscale::new_scale_colour()
p = p + ggplot2 :: geom_point(data = d2, ggplot2 :: aes(x = series, y = relative), shape = 21, fill = "white", color = "black", position = ggplot2 :: position_jitter(w = 0.4, h = 0), size = pchSize)
p = p + ggplot2 :: xlab("")
p = p + ggplot2 :: ylab(paste(fg$geneName[i], "(%)"))
p = p + ggplot2 :: theme(text= ggplot2 :: element_text(family="sans"),
legend.position = "none",
axis.title.x = ggplot2 :: element_text(size =textSize),
axis.title.y = ggplot2 :: element_text(size =textSize),
panel.background = ggplot2 :: element_rect(fill = NA),
panel.grid.major.y = ggplot2 :: element_line(colour = "grey90"),
panel.grid.minor.y = ggplot2 :: element_line(colour = "grey90"),
panel.grid.major.x = ggplot2 :: element_blank(),
axis.text.x = ggplot2 :: element_text(angle = 60, size=textSize, hjust = 1),
axis.text.y = ggplot2 :: element_text(size=textSize))#,
#p
ggplot2 :: ggsave(filename = paste0(folderStoredPlots, "cumulative_", paste(fg$geneName[i],fg$chr[i],fg$start[i],fg$end[i],"2", formatPlot ,sep="_",collapse = "")), plot = p, width = sizeInchesPlots[1], height = sizeInchesPlots[2],units = "in")
}
return(invisible(NULL))
}