|
0
|
1 args <- commandArgs(trailingOnly = TRUE)
|
|
|
2 options(scipen=999)
|
|
|
3
|
|
|
4 inFile = args[1]
|
|
|
5 outDir = args[2]
|
|
|
6 logfile = args[3]
|
|
|
7 min_freq = as.numeric(args[4])
|
|
|
8 min_cells = as.numeric(args[5])
|
|
|
9 mergeOn = args[6]
|
|
|
10
|
|
|
11 cat("<html><table><tr><td>Starting analysis</td></tr>", file=logfile, append=F)
|
|
|
12
|
|
|
13 library(ggplot2)
|
|
|
14 library(reshape2)
|
|
|
15 library(data.table)
|
|
|
16 library(grid)
|
|
|
17 library(parallel)
|
|
|
18 #require(xtable)
|
|
|
19 cat("<tr><td>Reading input</td></tr>", file=logfile, append=T)
|
|
|
20 dat = read.table(inFile, header=T, sep="\t", dec=".", fill=T, stringsAsFactors=F)
|
|
|
21 dat = dat[,c("Patient", "Receptor", "Sample", "Cell_Count", "Clone_Molecule_Count_From_Spikes", "Log10_Frequency", "Total_Read_Count", "J_Segment_Major_Gene", "V_Segment_Major_Gene", "CDR3_Sense_Sequence", "Clone_Sequence")]
|
|
|
22 dat$dsPerM = 0
|
|
|
23 dat = dat[!is.na(dat$Patient),]
|
|
|
24 dat$Related_to_leukemia_clone = F
|
|
|
25
|
|
|
26 setwd(outDir)
|
|
|
27 cat("<tr><td>Selecting first V/J Genes</td></tr>", file=logfile, append=T)
|
|
|
28 dat$V_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$V_Segment_Major_Gene), "; "), "[[", 1)))
|
|
|
29 dat$J_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$J_Segment_Major_Gene), "; "), "[[", 1)))
|
|
|
30
|
|
|
31 cat("<tr><td>Calculating Frequency</td></tr>", file=logfile, append=T)
|
|
|
32
|
|
|
33 dat$Frequency = ((10^dat$Log10_Frequency)*100)
|
|
|
34
|
|
|
35 dat = dat[dat$Frequency >= min_freq,]
|
|
|
36
|
|
|
37 triplets = dat[grepl("VanDongen_cALL_14696", dat$Patient) | grepl("(16278)|(26402)|(26759)", dat$Sample),]
|
|
|
38
|
|
|
39 cat("<tr><td>Normalizing to lowest cell count within locus</td></tr>", file=logfile, append=T)
|
|
|
40
|
|
|
41 dat$locus_V = substring(dat$V_Segment_Major_Gene, 0, 4)
|
|
|
42 dat$locus_J = substring(dat$J_Segment_Major_Gene, 0, 4)
|
|
|
43 min_cell_count = data.frame(data.table(dat)[, list(min_cell_count=min(.SD$Cell_Count)), by=c("Patient", "locus_V", "locus_J")])
|
|
|
44
|
|
|
45 dat$min_cell_paste = paste(dat$Patient, dat$locus_V, dat$locus_J)
|
|
|
46 min_cell_count$min_cell_paste = paste(min_cell_count$Patient, min_cell_count$locus_V, min_cell_count$locus_J)
|
|
|
47
|
|
|
48 min_cell_count = min_cell_count[,c("min_cell_paste", "min_cell_count")]
|
|
|
49 print(paste("rows:", nrow(dat)))
|
|
|
50 dat = merge(dat, min_cell_count, by="min_cell_paste")
|
|
|
51 print(paste("rows:", nrow(dat)))
|
|
|
52 dat$normalized_read_count = round(dat$Clone_Molecule_Count_From_Spikes / dat$Cell_Count * dat$min_cell_count / 2, digits=2) #??????????????????????????????????? wel of geen / 2
|
|
|
53
|
|
|
54 dat = dat[dat$normalized_read_count >= min_cells,]
|
|
|
55
|
|
|
56 dat$paste = paste(dat$Sample, dat$Clone_Sequence)
|
|
|
57
|
|
|
58 patients = split(dat, dat$Patient, drop=T)
|
|
|
59 intervalReads = rev(c(0,10,25,50,100,250,500,750,1000,10000))
|
|
|
60 intervalFreq = rev(c(0,0.01,0.05,0.1,0.5,1,5))
|
|
|
61 V_Segments = c(".*", "IGHV", "IGHD", "IGKV", "IGKV", "IgKINTR", "TRGV", "TRDV", "TRDD" , "TRBV")
|
|
|
62 J_Segments = c(".*", ".*", ".*", "IGKJ", "KDE", ".*", ".*", ".*", ".*", ".*")
|
|
|
63 Titles = c("Total", "IGH-Vh-Jh", "IGH-Dh-Jh", "Vk-Jk", "Vk-Kde" , "Intron-Kde", "TCRG", "TCRD-Vd-Dd", "TCRD-Dd-Dd", "TCRB-Vb-Jb")
|
|
|
64 Titles = factor(Titles, levels=Titles)
|
|
|
65 TitlesOrder = data.frame("Title"=Titles, "TitlesOrder"=1:length(Titles))
|
|
|
66
|
|
|
67 single_patients = data.frame("Patient" = character(0),"Sample" = character(0), "on" = character(0), "Clone_Sequence" = character(0), "Frequency" = numeric(0), "normalized_read_count" = numeric(0), "V_Segment_Major_Gene" = character(0), "J_Segment_Major_Gene" = character(0), "Rearrangement" = character(0))
|
|
|
68
|
|
|
69 patient.merge.list = list() #cache the 'both' table, 2x speedup for more memory...
|
|
|
70 patient.merge.list.second = list()
|
|
|
71 scatter_locus_data_list = list()
|
|
|
72 cat(paste("<table border='0' style='font-family:courier;'>", sep=""), file="multiple_matches.html", append=T)
|
|
|
73 cat(paste("<table border='0' style='font-family:courier;'>", sep=""), file="single_matches.html", append=T)
|
|
|
74 patientCountOnColumn <- function(x, product, interval, on, appendtxt=F){
|
|
|
75 if (!is.data.frame(x) & is.list(x)){
|
|
|
76 x = x[[1]]
|
|
|
77 }
|
|
|
78 #x$Sample = factor(x$Sample, levels=unique(x$Sample))
|
|
|
79 x = data.frame(x,stringsAsFactors=F)
|
|
|
80 onShort = "reads"
|
|
|
81 if(on == "Frequency"){
|
|
|
82 onShort = "freq"
|
|
|
83 }
|
|
|
84 onx = paste(on, ".x", sep="")
|
|
|
85 ony = paste(on, ".y", sep="")
|
|
|
86 splt = split(x, x$Sample, drop=T)
|
|
|
87 type="pair"
|
|
|
88 if(length(splt) == 1){
|
|
|
89 print(paste(paste(x[1,which(colnames(x) == "Patient")]), "has one sample"))
|
|
|
90 splt[[2]] = data.frame("Patient" = character(0), "Receptor" = character(0), "Sample" = character(0), "Cell_Count" = numeric(0), "Clone_Molecule_Count_From_Spikes" = numeric(0), "Log10_Frequency" = numeric(0), "Total_Read_Count" = numeric(0), "dsMol_per_1e6_cells" = numeric(0), "J_Segment_Major_Gene" = character(0), "V_Segment_Major_Gene" = character(0), "Clone_Sequence" = character(0), "CDR3_Sense_Sequence" = character(0), "Related_to_leukemia_clone" = logical(0), "Frequency"= numeric(0), "normalized_read_count" = numeric(0), "paste" = character(0))
|
|
|
91 type="single"
|
|
|
92 }
|
|
|
93 patient1 = splt[[1]]
|
|
|
94 patient2 = splt[[2]]
|
|
|
95
|
|
|
96 threshholdIndex = which(colnames(product) == "interval")
|
|
|
97 V_SegmentIndex = which(colnames(product) == "V_Segments")
|
|
|
98 J_SegmentIndex = which(colnames(product) == "J_Segments")
|
|
|
99 titleIndex = which(colnames(product) == "Titles")
|
|
|
100 sampleIndex = which(colnames(x) == "Sample")
|
|
|
101 patientIndex = which(colnames(x) == "Patient")
|
|
|
102 oneSample = paste(patient1[1,sampleIndex], sep="")
|
|
|
103 twoSample = paste(patient2[1,sampleIndex], sep="")
|
|
|
104 patient = paste(x[1,patientIndex])
|
|
|
105
|
|
|
106 switched = F
|
|
|
107 if(length(grep(".*_Right$", twoSample)) == 1 || length(grep(".*_Dx_BM$", twoSample)) == 1 || length(grep(".*_Dx$", twoSample)) == 1 ){
|
|
|
108 tmp = twoSample
|
|
|
109 twoSample = oneSample
|
|
|
110 oneSample = tmp
|
|
|
111 tmp = patient1
|
|
|
112 patient1 = patient2
|
|
|
113 patient2 = tmp
|
|
|
114 switched = T
|
|
|
115 }
|
|
|
116 if(appendtxt){
|
|
|
117 cat(paste(patient, oneSample, twoSample, type, sep="\t"), file="patients.txt", append=T, sep="", fill=3)
|
|
|
118 }
|
|
|
119 cat(paste("<tr><td>", patient, "</td>", sep=""), file=logfile, append=T)
|
|
|
120
|
|
|
121 if(mergeOn == "Clone_Sequence"){
|
|
|
122 patient1$merge = paste(patient1$Clone_Sequence)
|
|
|
123 patient2$merge = paste(patient2$Clone_Sequence)
|
|
|
124 } else {
|
|
|
125 patient1$merge = paste(patient1$V_Segment_Major_Gene, patient1$J_Segment_Major_Gene, patient1$CDR3_Sense_Sequence)
|
|
|
126 patient2$merge = paste(patient2$V_Segment_Major_Gene, patient2$J_Segment_Major_Gene, patient2$CDR3_Sense_Sequence)
|
|
|
127 }
|
|
|
128
|
|
|
129 scatterplot_data_columns = c("Patient", "Sample", "Frequency", "normalized_read_count", "V_Segment_Major_Gene", "J_Segment_Major_Gene", "merge")
|
|
|
130 #scatterplot_data = rbind(patient1[,scatterplot_data_columns], patient2[,scatterplot_data_columns])
|
|
|
131 scatterplot_data = patient1[NULL,scatterplot_data_columns]
|
|
|
132 #scatterplot_data = scatterplot_data[!duplicated(scatterplot_data$merge),]
|
|
|
133 #scatterplot_data$type = factor(x=oneSample, levels=c(oneSample, twoSample, "In Both"))
|
|
|
134 scatterplot.data.type.factor = c(oneSample, twoSample, paste(c(oneSample, twoSample), "In Both"))
|
|
|
135 #scatterplot_data$type = factor(x=NULL, levels=scatterplot.data.type.factor)
|
|
|
136 scatterplot_data$type = character(0)
|
|
|
137 scatterplot_data$link = numeric(0)
|
|
|
138 scatterplot_data$on = character(0)
|
|
|
139
|
|
|
140 #patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge") #merge alles 'fuzzy'
|
|
|
141 patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge")[NULL,] #blegh
|
|
|
142
|
|
|
143 cs.exact.matches = patient1[patient1$Clone_Sequence %in% patient2$Clone_Sequence,]$Clone_Sequence
|
|
|
144
|
|
|
145 start.time = proc.time()
|
|
|
146 merge.list = c()
|
|
|
147
|
|
|
148 if(patient %in% names(patient.merge.list)){
|
|
|
149 patientMerge = patient.merge.list[[patient]]
|
|
|
150 merge.list[["second"]] = patient.merge.list.second[[patient]]
|
|
|
151 scatterplot_data = scatter_locus_data_list[[patient]]
|
|
|
152 cat(paste("<td>", nrow(patient1), " in ", oneSample, " and ", nrow(patient2), " in ", twoSample, ", ", nrow(patientMerge), " in both (fetched from cache)</td></tr>", sep=""), file=logfile, append=T)
|
|
|
153
|
|
|
154 print(names(patient.merge.list))
|
|
|
155 } else {
|
|
|
156 #fuzzy matching here...
|
|
|
157 #merge.list = patientMerge$merge
|
|
|
158
|
|
|
159 #patient1.fuzzy = patient1[!(patient1$merge %in% merge.list),]
|
|
|
160 #patient2.fuzzy = patient2[!(patient2$merge %in% merge.list),]
|
|
|
161
|
|
|
162 patient1.fuzzy = patient1
|
|
|
163 patient2.fuzzy = patient2
|
|
|
164
|
|
|
165 #patient1.fuzzy$merge = paste(patient1.fuzzy$V_Segment_Major_Gene, patient1.fuzzy$J_Segment_Major_Gene, patient1.fuzzy$CDR3_Sense_Sequence)
|
|
|
166 #patient2.fuzzy$merge = paste(patient2.fuzzy$V_Segment_Major_Gene, patient2.fuzzy$J_Segment_Major_Gene, patient2.fuzzy$CDR3_Sense_Sequence)
|
|
|
167
|
|
|
168 #patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J, patient1.fuzzy$CDR3_Sense_Sequence)
|
|
|
169 #patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J, patient2.fuzzy$CDR3_Sense_Sequence)
|
|
|
170
|
|
|
171 patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J)
|
|
|
172 patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J)
|
|
|
173
|
|
|
174 #merge.freq.table = data.frame(table(c(patient1.fuzzy[!duplicated(patient1.fuzzy$merge),"merge"], patient2.fuzzy[!duplicated(patient2.fuzzy$merge),"merge"]))) #also remove?
|
|
|
175 #merge.freq.table.gt.1 = merge.freq.table[merge.freq.table$Freq > 1,]
|
|
|
176
|
|
|
177 #patient1.fuzzy = patient1.fuzzy[patient1.fuzzy$merge %in% merge.freq.table.gt.1$Var1,]
|
|
|
178 #patient2.fuzzy = patient2.fuzzy[patient2.fuzzy$merge %in% merge.freq.table.gt.1$Var1,]
|
|
|
179
|
|
|
180 patient.fuzzy = rbind(patient1.fuzzy, patient2.fuzzy)
|
|
|
181 patient.fuzzy = patient.fuzzy[order(nchar(patient.fuzzy$Clone_Sequence)),]
|
|
|
182
|
|
|
183 merge.list = list()
|
|
|
184
|
|
|
185 merge.list[["second"]] = vector()
|
|
|
186
|
|
|
187 link.count = 1
|
|
|
188
|
|
|
189 while(nrow(patient.fuzzy) > 1){
|
|
|
190 first.merge = patient.fuzzy[1,"merge"]
|
|
|
191 first.clone.sequence = patient.fuzzy[1,"Clone_Sequence"]
|
|
|
192 first.sample = patient.fuzzy[1,"Sample"]
|
|
|
193 merge.filter = first.merge == patient.fuzzy$merge
|
|
|
194
|
|
|
195 #length.filter = nchar(patient.fuzzy$Clone_Sequence) - nchar(first.clone.sequence) <= 9
|
|
|
196
|
|
|
197 first.sample.filter = first.sample == patient.fuzzy$Sample
|
|
|
198 second.sample.filter = first.sample != patient.fuzzy$Sample
|
|
|
199
|
|
|
200 #first match same sample, sum to a single row, same for other sample
|
|
|
201 #then merge rows like 'normal'
|
|
|
202
|
|
|
203 sequence.filter = grepl(paste("^", first.clone.sequence, sep=""), patient.fuzzy$Clone_Sequence)
|
|
|
204
|
|
|
205
|
|
|
206
|
|
|
207 #match.filter = merge.filter & grepl(first.clone.sequence, patient.fuzzy$Clone_Sequence) & length.filter & sample.filter
|
|
|
208 first.match.filter = merge.filter & sequence.filter & first.sample.filter
|
|
|
209 second.match.filter = merge.filter & sequence.filter & second.sample.filter
|
|
|
210
|
|
|
211 first.rows = patient.fuzzy[first.match.filter,]
|
|
|
212 second.rows = patient.fuzzy[second.match.filter,]
|
|
|
213
|
|
|
214 first.rows.v = table(first.rows$V_Segment_Major_Gene)
|
|
|
215 first.rows.v = names(first.rows.v[which.max(first.rows.v)])
|
|
|
216 first.rows.j = table(first.rows$J_Segment_Major_Gene)
|
|
|
217 first.rows.j = names(first.rows.j[which.max(first.rows.j)])
|
|
|
218
|
|
|
219 first.sum = data.frame(merge = first.clone.sequence,
|
|
|
220 Patient = patient,
|
|
|
221 Receptor = first.rows[1,"Receptor"],
|
|
|
222 Sample = first.rows[1,"Sample"],
|
|
|
223 Cell_Count = first.rows[1,"Cell_Count"],
|
|
|
224 Clone_Molecule_Count_From_Spikes = sum(first.rows$Clone_Molecule_Count_From_Spikes),
|
|
|
225 Log10_Frequency = log10(sum(first.rows$Frequency)),
|
|
|
226 Total_Read_Count = sum(first.rows$Total_Read_Count),
|
|
|
227 dsPerM = sum(first.rows$dsPerM),
|
|
|
228 J_Segment_Major_Gene = first.rows.j,
|
|
|
229 V_Segment_Major_Gene = first.rows.v,
|
|
|
230 Clone_Sequence = first.clone.sequence,
|
|
|
231 CDR3_Sense_Sequence = first.rows[1,"CDR3_Sense_Sequence"],
|
|
|
232 Related_to_leukemia_clone = F,
|
|
|
233 Frequency = sum(first.rows$Frequency),
|
|
|
234 locus_V = first.rows[1,"locus_V"],
|
|
|
235 locus_J = first.rows[1,"locus_J"],
|
|
|
236 min_cell_count = first.rows[1,"min_cell_count"],
|
|
|
237 normalized_read_count = sum(first.rows$normalized_read_count),
|
|
|
238 paste = first.rows[1,"paste"],
|
|
|
239 min_cell_paste = first.rows[1,"min_cell_paste"])
|
|
|
240
|
|
|
241 if(nrow(second.rows) > 0){
|
|
|
242 second.rows.v = table(second.rows$V_Segment_Major_Gene)
|
|
|
243 second.rows.v = names(second.rows.v[which.max(second.rows.v)])
|
|
|
244 second.rows.j = table(second.rows$J_Segment_Major_Gene)
|
|
|
245 second.rows.j = names(second.rows.j[which.max(second.rows.j)])
|
|
|
246
|
|
|
247 second.sum = data.frame(merge = first.clone.sequence,
|
|
|
248 Patient = patient,
|
|
|
249 Receptor = second.rows[1,"Receptor"],
|
|
|
250 Sample = second.rows[1,"Sample"],
|
|
|
251 Cell_Count = second.rows[1,"Cell_Count"],
|
|
|
252 Clone_Molecule_Count_From_Spikes = sum(second.rows$Clone_Molecule_Count_From_Spikes),
|
|
|
253 Log10_Frequency = log10(sum(second.rows$Frequency)),
|
|
|
254 Total_Read_Count = sum(second.rows$Total_Read_Count),
|
|
|
255 dsPerM = sum(second.rows$dsPerM),
|
|
|
256 J_Segment_Major_Gene = second.rows.j,
|
|
|
257 V_Segment_Major_Gene = second.rows.v,
|
|
|
258 Clone_Sequence = first.clone.sequence,
|
|
|
259 CDR3_Sense_Sequence = second.rows[1,"CDR3_Sense_Sequence"],
|
|
|
260 Related_to_leukemia_clone = F,
|
|
|
261 Frequency = sum(second.rows$Frequency),
|
|
|
262 locus_V = second.rows[1,"locus_V"],
|
|
|
263 locus_J = second.rows[1,"locus_J"],
|
|
|
264 min_cell_count = second.rows[1,"min_cell_count"],
|
|
|
265 normalized_read_count = sum(second.rows$normalized_read_count),
|
|
|
266 paste = second.rows[1,"paste"],
|
|
|
267 min_cell_paste = second.rows[1,"min_cell_paste"])
|
|
|
268
|
|
|
269 patientMerge = rbind(patientMerge, merge(first.sum, second.sum, by="merge"))
|
|
|
270 patient.fuzzy = patient.fuzzy[!(first.match.filter | second.match.filter),]
|
|
|
271
|
|
|
272 hidden.clone.sequences = c(first.rows[-1,"Clone_Sequence"], second.rows[second.rows$Clone_Sequence != first.clone.sequence,"Clone_Sequence"])
|
|
|
273 merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences)
|
|
|
274
|
|
|
275 tmp.rows = rbind(first.rows, second.rows)
|
|
|
276 tmp.rows = tmp.rows[order(nchar(tmp.rows$Clone_Sequence)),]
|
|
|
277
|
|
|
278
|
|
|
279 #add to the scatterplot data
|
|
|
280 scatterplot.row = first.sum[,scatterplot_data_columns]
|
|
|
281 scatterplot.row$type = paste(first.sum[,"Sample"], "In Both")
|
|
|
282 scatterplot.row$link = link.count
|
|
|
283 scatterplot.row$on = onShort
|
|
|
284
|
|
|
285 scatterplot_data = rbind(scatterplot_data, scatterplot.row)
|
|
|
286
|
|
|
287 scatterplot.row = second.sum[,scatterplot_data_columns]
|
|
|
288 scatterplot.row$type = paste(second.sum[,"Sample"], "In Both")
|
|
|
289 scatterplot.row$link = link.count
|
|
|
290 scatterplot.row$on = onShort
|
|
|
291
|
|
|
292 scatterplot_data = rbind(scatterplot_data, scatterplot.row)
|
|
|
293
|
|
|
294 #write some information about the match to a log file
|
|
|
295 if (nrow(first.rows) > 1 | nrow(second.rows) > 1) {
|
|
|
296 cat(paste("<tr><td>", patient, " row ", 1:nrow(tmp.rows), "</td><td>", tmp.rows$Sample, ":</td><td>", tmp.rows$Clone_Sequence, "</td><td>", tmp.rows$normalized_read_count, "</td></tr>", sep=""), file="multiple_matches.html", append=T)
|
|
|
297 } else {
|
|
|
298 second.clone.sequence = second.rows[1,"Clone_Sequence"]
|
|
|
299 if(nchar(first.clone.sequence) != nchar(second.clone.sequence)){
|
|
|
300 cat(paste("<tr bgcolor='#DDD'><td>", patient, " row ", 1:nrow(tmp.rows), "</td><td>", tmp.rows$Sample, ":</td><td>", tmp.rows$Clone_Sequence, "</td><td>", tmp.rows$normalized_read_count, "</td></tr>", sep=""), file="single_matches.html", append=T)
|
|
|
301 } else {
|
|
|
302 #cat(paste("<tr><td>", patient, " row ", 1:nrow(tmp.rows), "</td><td>", tmp.rows$Sample, ":</td><td>", tmp.rows$Clone_Sequence, "</td><td>", tmp.rows$normalized_read_count, "</td></tr>", sep=""), file="single_matches.html", append=T)
|
|
|
303 }
|
|
|
304 }
|
|
|
305
|
|
|
306 } else if(nrow(first.rows) > 1) {
|
|
|
307 if(patient1[1,"Sample"] == first.sample){
|
|
|
308 patient1 = patient1[!(patient1$Clone_Sequence %in% first.rows$Clone_Sequence),]
|
|
|
309 patient1 = rbind(patient1, first.sum)
|
|
|
310 } else {
|
|
|
311 patient2 = patient2[!(patient2$Clone_Sequence %in% first.rows$Clone_Sequence),]
|
|
|
312 patient2 = rbind(patient2, first.sum)
|
|
|
313 }
|
|
|
314
|
|
|
315 hidden.clone.sequences = c(first.rows[-1,"Clone_Sequence"])
|
|
|
316 merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences)
|
|
|
317
|
|
|
318 patient.fuzzy = patient.fuzzy[-first.match.filter,]
|
|
|
319
|
|
|
320 #add to the scatterplot data
|
|
|
321 scatterplot.row = first.sum[,scatterplot_data_columns]
|
|
|
322 scatterplot.row$type = first.sum[,"Sample"]
|
|
|
323 scatterplot.row$link = link.count
|
|
|
324 scatterplot.row$on = onShort
|
|
|
325
|
|
|
326 scatterplot_data = rbind(scatterplot_data, scatterplot.row)
|
|
|
327
|
|
|
328 cat(paste("<tr bgcolor='#DDF'><td>", patient, " row ", 1:nrow(first.rows), "</td><td>", first.rows$Sample, ":</td><td>", first.rows$Clone_Sequence, "</td><td>", first.rows$normalized_read_count, "</td></tr>", sep=""), file="single_matches.html", append=T)
|
|
|
329 } else {
|
|
|
330 patient.fuzzy = patient.fuzzy[-1,]
|
|
|
331
|
|
|
332 #add to the scatterplot data
|
|
|
333 scatterplot.row = first.sum[,scatterplot_data_columns]
|
|
|
334 scatterplot.row$type = first.sum[,"Sample"]
|
|
|
335 scatterplot.row$link = link.count
|
|
|
336 scatterplot.row$on = onShort
|
|
|
337
|
|
|
338 scatterplot_data = rbind(scatterplot_data, scatterplot.row)
|
|
|
339 }
|
|
|
340 link.count = link.count + 1
|
|
|
341 }
|
|
|
342 patient.merge.list[[patient]] <<- patientMerge
|
|
|
343 patient.merge.list.second[[patient]] <<- merge.list[["second"]]
|
|
|
344
|
|
|
345 sample.order = data.frame(type = c(oneSample, twoSample, paste(c(oneSample, twoSample), "In Both")),type.order = 1:4)
|
|
|
346 scatterplot_data = merge(scatterplot_data, sample.order, by="type")
|
|
|
347
|
|
|
348 scatter_locus_data_list[[patient]] <<- scatterplot_data
|
|
|
349 cat(paste("<td>", nrow(patient1), " in ", oneSample, " and ", nrow(patient2), " in ", twoSample, ", ", nrow(patientMerge), " in both (finding both took ", (proc.time() - start.time)[[3]], "s)</td></tr>", sep=""), file=logfile, append=T)
|
|
|
350 }
|
|
|
351
|
|
|
352 patient1 = patient1[!(patient1$Clone_Sequence %in% patient.merge.list.second[[patient]]),]
|
|
|
353 patient2 = patient2[!(patient2$Clone_Sequence %in% patient.merge.list.second[[patient]]),]
|
|
|
354
|
|
|
355
|
|
|
356 patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony])
|
|
|
357 #patientMerge$thresholdValue = pmin(patientMerge[,onx], patientMerge[,ony])
|
|
|
358 res1 = vector()
|
|
|
359 res2 = vector()
|
|
|
360 resBoth = vector()
|
|
|
361 read1Count = vector()
|
|
|
362 read2Count = vector()
|
|
|
363 locussum1 = vector()
|
|
|
364 locussum2 = vector()
|
|
|
365
|
|
|
366 #for(iter in 1){
|
|
|
367 for(iter in 1:length(product[,1])){
|
|
|
368 threshhold = product[iter,threshholdIndex]
|
|
|
369 V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="")
|
|
|
370 J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="")
|
|
|
371 #both = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge[,onx] > threshhold & patientMerge[,ony] > threshhold) #both higher than threshold
|
|
|
372 both = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge$thresholdValue > threshhold) #highest of both is higher than threshold
|
|
|
373 one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$merge %in% patientMerge[both,]$merge))
|
|
|
374 two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$merge %in% patientMerge[both,]$merge))
|
|
|
375 read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count))
|
|
|
376 read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count))
|
|
|
377 res1 = append(res1, sum(one))
|
|
|
378 res2 = append(res2, sum(two))
|
|
|
379 resBoth = append(resBoth, sum(both))
|
|
|
380 locussum1 = append(locussum1, sum(patient1[(grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene)),]$normalized_read_count))
|
|
|
381 locussum2 = append(locussum2, sum(patient2[(grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene)),]$normalized_read_count))
|
|
|
382 #threshhold = 0
|
|
|
383 if(threshhold != 0){
|
|
|
384 if(sum(one) > 0){
|
|
|
385 dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
|
|
|
386 colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone Sequence", "Related_to_leukemia_clone")
|
|
|
387 filenameOne = paste(oneSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
388 write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
389 }
|
|
|
390 if(sum(two) > 0){
|
|
|
391 dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
|
|
|
392 colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone Sequence", "Related_to_leukemia_clone")
|
|
|
393 filenameTwo = paste(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
394 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
395 }
|
|
|
396 } else {
|
|
|
397 scatterplot_locus_data = scatterplot_data[grepl(V_Segment, scatterplot_data$V_Segment_Major_Gene) & grepl(J_Segment, scatterplot_data$J_Segment_Major_Gene),]
|
|
|
398 if(nrow(scatterplot_locus_data) > 0){
|
|
|
399 scatterplot_locus_data$Rearrangement = product[iter, titleIndex]
|
|
|
400 }
|
|
|
401
|
|
|
402
|
|
|
403
|
|
|
404 p = NULL
|
|
|
405 print(paste("nrow scatterplot_locus_data", nrow(scatterplot_locus_data)))
|
|
|
406 if(nrow(scatterplot_locus_data) != 0){
|
|
|
407 if(on == "normalized_read_count"){
|
|
|
408 write.table(scatterplot_locus_data, file=paste(oneSample, twoSample, product[iter, titleIndex], "scatterplot_locus_data.txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
409 scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count))))
|
|
|
410 p = ggplot(scatterplot_locus_data, aes(factor(reorder(type, type.order)), normalized_read_count, group=link)) + geom_line() + scale_y_log10(breaks=scales,labels=scales, limits=c(1,1e6)) + scale_x_discrete(breaks=levels(scatterplot_data$type), labels=levels(scatterplot_data$type), drop=FALSE)
|
|
|
411 } else {
|
|
|
412 p = ggplot(scatterplot_locus_data, aes(factor(reorder(type, type.order)), Frequency, group=link)) + geom_line() + scale_y_log10(limits=c(0.0001,100), breaks=c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100), labels=c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100")) + scale_x_discrete(breaks=levels(scatterplot_data$type), labels=levels(scatterplot_data$type), drop=FALSE)
|
|
|
413 }
|
|
|
414 p = p + geom_point(aes(colour=type), position="dodge")
|
|
|
415 p = p + xlab("In one or both samples") + ylab(onShort) + ggtitle(paste(patient1[1,patientIndex], patient1[1,sampleIndex], patient2[1,sampleIndex], onShort, product[iter, titleIndex]))
|
|
|
416 } else {
|
|
|
417 p = ggplot(NULL, aes(x=c("In one", "In Both"),y=0)) + geom_blank(NULL) + xlab("In one or both of the samples") + ylab(onShort) + ggtitle(paste(patient1[1,patientIndex], patient1[1,sampleIndex], patient2[1,sampleIndex], onShort, product[iter, titleIndex]))
|
|
|
418 }
|
|
|
419 png(paste(patient1[1,patientIndex], "_", patient1[1,sampleIndex], "_", patient2[1,sampleIndex], "_", onShort, "_", product[iter, titleIndex],"_scatter.png", sep=""))
|
|
|
420 print(p)
|
|
|
421 dev.off()
|
|
|
422 }
|
|
|
423 if(sum(both) > 0){
|
|
|
424 dfBoth = patientMerge[both,c("V_Segment_Major_Gene.x", "J_Segment_Major_Gene.x", "normalized_read_count.x", "Frequency.x", "Related_to_leukemia_clone.x", "Clone_Sequence.x", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "normalized_read_count.y", "Frequency.y", "Related_to_leukemia_clone.y")]
|
|
|
425 colnames(dfBoth) = c(paste("Proximal segment", oneSample), paste("Distal segment", oneSample), paste("Normalized_Read_Count", oneSample), paste("Frequency", oneSample), paste("Related_to_leukemia_clone", oneSample),"Clone Sequence", paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Normalized_Read_Count", twoSample), paste("Frequency", twoSample), paste("Related_to_leukemia_clone", twoSample))
|
|
|
426 filenameBoth = paste(oneSample, "_", twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
427 write.table(dfBoth, file=paste(filenameBoth, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
428 }
|
|
|
429 }
|
|
|
430 patientResult = data.frame("Locus"=product$Titles, "J_Segment"=product$J_Segments, "V_Segment"=product$V_Segments, "cut_off_value"=paste(">", product$interval, sep=""), "Both"=resBoth, "tmp1"=res1, "read_count1" = round(read1Count), "tmp2"=res2, "read_count2"= round(read2Count), "Sum"=res1 + res2 + resBoth, "percentage" = round((resBoth/(res1 + res2 + resBoth)) * 100, digits=2), "Locus_sum1"=locussum1, "Locus_sum2"=locussum2)
|
|
|
431 if(sum(is.na(patientResult$percentage)) > 0){
|
|
|
432 patientResult[is.na(patientResult$percentage),]$percentage = 0
|
|
|
433 }
|
|
|
434 colnames(patientResult)[6] = oneSample
|
|
|
435 colnames(patientResult)[8] = twoSample
|
|
|
436 colnamesBak = colnames(patientResult)
|
|
|
437 colnames(patientResult) = c("Ig/TCR gene rearrangement type", "Distal Gene segment", "Proximal gene segment", "cut_off_value", paste("Number of sequences ", patient, "_Both", sep=""), paste("Number of sequences", oneSample, sep=""), paste("Normalized Read Count", oneSample), paste("Number of sequences", twoSample, sep=""), paste("Normalized Read Count", twoSample), paste("Sum number of sequences", patient), paste("Percentage of sequences ", patient, "_Both", sep=""), paste("Locus Sum", oneSample), paste("Locus Sum", twoSample))
|
|
|
438 write.table(patientResult, file=paste(patient, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
439 colnames(patientResult) = colnamesBak
|
|
|
440
|
|
|
441 patientResult$Locus = factor(patientResult$Locus, Titles)
|
|
|
442 patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))
|
|
|
443
|
|
|
444 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "Both")])
|
|
|
445 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=Both), stat='identity', position="dodge", fill="#79c36a")
|
|
|
446 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
447 plt = plt + geom_text(aes(ymax=max(Both), x=cut_off_value,y=Both,label=Both), angle=90, hjust=0)
|
|
|
448 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in both")
|
|
|
449 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
|
450 png(paste(patient, "_", onShort, ".png", sep=""), width=1920, height=1080)
|
|
|
451 print(plt)
|
|
|
452 dev.off()
|
|
|
453 #(t,r,b,l)
|
|
|
454 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "percentage")])
|
|
|
455 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=percentage), stat='identity', position="dodge", fill="#79c36a")
|
|
|
456 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
457 plt = plt + geom_text(aes(ymax=max(percentage), x=cut_off_value,y=percentage,label=percentage), angle=90, hjust=0)
|
|
|
458 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("% clones in both left and right")
|
|
|
459 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
|
460 png(paste(patient, "_percent_", onShort, ".png", sep=""), width=1920, height=1080)
|
|
|
461 print(plt)
|
|
|
462 dev.off()
|
|
|
463
|
|
|
464 patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample)] ,id.vars=1:2)
|
|
|
465 patientResult$relativeValue = patientResult$value * 10
|
|
|
466 patientResult[patientResult$relativeValue == 0,]$relativeValue = 1
|
|
|
467 plt = ggplot(patientResult)
|
|
|
468 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge")
|
|
|
469 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
470 plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9)))
|
|
|
471 plt = plt + geom_text(data=patientResult[patientResult$variable == oneSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=-0.2)
|
|
|
472 plt = plt + geom_text(data=patientResult[patientResult$variable == twoSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=0.8)
|
|
|
473 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle(paste("Number of clones in only ", oneSample, " and only ", twoSample, sep=""))
|
|
|
474 png(paste(patient, "_", onShort, "_both.png", sep=""), width=1920, height=1080)
|
|
|
475 print(plt)
|
|
|
476 dev.off()
|
|
|
477 }
|
|
|
478
|
|
|
479 cat("<tr><td>Starting Frequency analysis</td></tr>", file=logfile, append=T)
|
|
|
480
|
|
|
481 interval = intervalFreq
|
|
|
482 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
|
483 product = data.frame("Titles"=rep(Titles, each=length(interval)), "interval"=rep(interval, times=10), "V_Segments"=rep(V_Segments, each=length(interval)), "J_Segments"=rep(J_Segments, each=length(interval)))
|
|
|
484 lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
|
|
|
485
|
|
|
486 cat("<tr><td>Starting Cell Count analysis</td></tr>", file=logfile, append=T)
|
|
|
487
|
|
|
488 interval = intervalReads
|
|
|
489 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
|
490 product = data.frame("Titles"=rep(Titles, each=length(interval)), "interval"=rep(interval, times=10), "V_Segments"=rep(V_Segments, each=length(interval)), "J_Segments"=rep(J_Segments, each=length(interval)))
|
|
|
491 lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="normalized_read_count")
|
|
|
492
|
|
|
493 if(nrow(single_patients) > 0){
|
|
|
494 scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count))))
|
|
|
495 p = ggplot(single_patients, aes(Rearrangement, normalized_read_count)) + scale_y_log10(breaks=scales,labels=as.character(scales)) + expand_limits(y=c(0,1000000))
|
|
|
496 p = p + geom_point(aes(colour=type), position="jitter")
|
|
|
497 p = p + xlab("In one or both samples") + ylab("Reads")
|
|
|
498 p = p + facet_grid(.~Patient) + ggtitle("Scatterplot of the reads of the patients with a single sample")
|
|
|
499 png("singles_reads_scatterplot.png", width=640 * length(unique(single_patients$Patient)) + 100, height=1080)
|
|
|
500 print(p)
|
|
|
501 dev.off()
|
|
|
502
|
|
|
503 #p = ggplot(single_patients, aes(Rearrangement, Frequency)) + scale_y_continuous(limits = c(0, 100)) + expand_limits(y=c(0,100))
|
|
|
504 p = ggplot(single_patients, aes(Rearrangement, Frequency)) + scale_y_log10(limits=c(0.0001,100), breaks=c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100), labels=c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100")) + expand_limits(y=c(0,100))
|
|
|
505 p = p + geom_point(aes(colour=type), position="jitter")
|
|
|
506 p = p + xlab("In one or both samples") + ylab("Frequency")
|
|
|
507 p = p + facet_grid(.~Patient) + ggtitle("Scatterplot of the frequency of the patients with a single sample")
|
|
|
508 png("singles_freq_scatterplot.png", width=640 * length(unique(single_patients$Patient)) + 100, height=1080)
|
|
|
509 print(p)
|
|
|
510 dev.off()
|
|
|
511 } else {
|
|
|
512 empty <- data.frame()
|
|
|
513 p = ggplot(empty) + geom_point() + xlim(0, 10) + ylim(0, 100) + xlab("In one or both samples") + ylab("Frequency") + ggtitle("Scatterplot of the frequency of the patients with a single sample")
|
|
|
514
|
|
|
515 png("singles_reads_scatterplot.png", width=400, height=300)
|
|
|
516 print(p)
|
|
|
517 dev.off()
|
|
|
518
|
|
|
519 png("singles_freq_scatterplot.png", width=400, height=300)
|
|
|
520 print(p)
|
|
|
521 dev.off()
|
|
|
522 }
|
|
|
523
|
|
|
524 patient.merge.list = list() #cache the 'both' table, 2x speedup for more memory...
|
|
|
525 patient.merge.list.second = list()
|
|
|
526
|
|
|
527 tripletAnalysis <- function(patient1, label1, patient2, label2, patient3, label3, product, interval, on, appendTriplets= FALSE){
|
|
|
528 onShort = "reads"
|
|
|
529 if(on == "Frequency"){
|
|
|
530 onShort = "freq"
|
|
|
531 }
|
|
|
532 onx = paste(on, ".x", sep="")
|
|
|
533 ony = paste(on, ".y", sep="")
|
|
|
534 onz = paste(on, ".z", sep="")
|
|
|
535 type="triplet"
|
|
|
536
|
|
|
537 threshholdIndex = which(colnames(product) == "interval")
|
|
|
538 V_SegmentIndex = which(colnames(product) == "V_Segments")
|
|
|
539 J_SegmentIndex = which(colnames(product) == "J_Segments")
|
|
|
540 titleIndex = which(colnames(product) == "Titles")
|
|
|
541 sampleIndex = which(colnames(patient1) == "Sample")
|
|
|
542 patientIndex = which(colnames(patient1) == "Patient")
|
|
|
543 oneSample = paste(patient1[1,sampleIndex], sep="")
|
|
|
544 twoSample = paste(patient2[1,sampleIndex], sep="")
|
|
|
545 threeSample = paste(patient3[1,sampleIndex], sep="")
|
|
|
546
|
|
|
547 if(mergeOn == "Clone_Sequence"){
|
|
|
548 patient1$merge = paste(patient1$Clone_Sequence)
|
|
|
549 patient2$merge = paste(patient2$Clone_Sequence)
|
|
|
550 patient3$merge = paste(patient3$Clone_Sequence)
|
|
|
551
|
|
|
552 } else {
|
|
|
553 patient1$merge = paste(patient1$V_Segment_Major_Gene, patient1$J_Segment_Major_Gene, patient1$CDR3_Sense_Sequence)
|
|
|
554 patient2$merge = paste(patient2$V_Segment_Major_Gene, patient2$J_Segment_Major_Gene, patient2$CDR3_Sense_Sequence)
|
|
|
555 patient3$merge = paste(patient3$V_Segment_Major_Gene, patient3$J_Segment_Major_Gene, patient3$CDR3_Sense_Sequence)
|
|
|
556 }
|
|
|
557
|
|
|
558 #patientMerge = merge(patient1, patient2, by="merge")[NULL,]
|
|
|
559 patient1.fuzzy = patient1
|
|
|
560 patient2.fuzzy = patient2
|
|
|
561 patient3.fuzzy = patient3
|
|
|
562
|
|
|
563 cat(paste("<tr><td>", label1, "</td>", sep=""), file=logfile, append=T)
|
|
|
564
|
|
|
565 patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J)
|
|
|
566 patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J)
|
|
|
567 patient3.fuzzy$merge = paste(patient3.fuzzy$locus_V, patient3.fuzzy$locus_J)
|
|
|
568
|
|
|
569 patient.fuzzy = rbind(patient1.fuzzy ,patient2.fuzzy, patient3.fuzzy)
|
|
|
570 patient.fuzzy = patient.fuzzy[order(nchar(patient.fuzzy$Clone_Sequence)),]
|
|
|
571
|
|
|
572 other.sample.list = list()
|
|
|
573 other.sample.list[[oneSample]] = c(twoSample, threeSample)
|
|
|
574 other.sample.list[[twoSample]] = c(oneSample, threeSample)
|
|
|
575 other.sample.list[[threeSample]] = c(oneSample, twoSample)
|
|
|
576
|
|
|
577 patientMerge = merge(patient1, patient2, by="merge")
|
|
|
578 patientMerge = merge(patientMerge, patient3, by="merge")
|
|
|
579 colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge)))] = paste(colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge), perl=T))], ".z", sep="")
|
|
|
580 #patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz])
|
|
|
581 patientMerge = patientMerge[NULL,]
|
|
|
582
|
|
|
583 duo.merge.list = list()
|
|
|
584
|
|
|
585 patientMerge12 = merge(patient1, patient2, by="merge")
|
|
|
586 #patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony])
|
|
|
587 patientMerge12 = patientMerge12[NULL,]
|
|
|
588 duo.merge.list[[paste(oneSample, twoSample)]] = patientMerge12
|
|
|
589 duo.merge.list[[paste(twoSample, oneSample)]] = patientMerge12
|
|
|
590
|
|
|
591 patientMerge13 = merge(patient1, patient3, by="merge")
|
|
|
592 #patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony])
|
|
|
593 patientMerge13 = patientMerge13[NULL,]
|
|
|
594 duo.merge.list[[paste(oneSample, threeSample)]] = patientMerge13
|
|
|
595 duo.merge.list[[paste(threeSample, oneSample)]] = patientMerge13
|
|
|
596
|
|
|
597 patientMerge23 = merge(patient2, patient3, by="merge")
|
|
|
598 #patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony])
|
|
|
599 patientMerge23 = patientMerge23[NULL,]
|
|
|
600 duo.merge.list[[paste(twoSample, threeSample)]] = patientMerge23
|
|
|
601 duo.merge.list[[paste(threeSample, twoSample)]] = patientMerge23
|
|
|
602
|
|
|
603 merge.list = list()
|
|
|
604 merge.list[["second"]] = vector()
|
|
|
605
|
|
|
606 start.time = proc.time()
|
|
|
607 if(paste(label1, "123") %in% names(patient.merge.list)){
|
|
|
608 patientMerge = patient.merge.list[[paste(label1, "123")]]
|
|
|
609 patientMerge12 = patient.merge.list[[paste(label1, "12")]]
|
|
|
610 patientMerge13 = patient.merge.list[[paste(label1, "13")]]
|
|
|
611 patientMerge23 = patient.merge.list[[paste(label1, "23")]]
|
|
|
612
|
|
|
613 merge.list[["second"]] = patient.merge.list.second[[label1]]
|
|
|
614
|
|
|
615 cat(paste("<td>", nrow(patient1), " in ", label1, " and ", nrow(patient2), " in ", label2, nrow(patient3), " in ", label3, ", ", nrow(patientMerge), " in both (fetched from cache)</td></tr>", sep=""), file=logfile, append=T)
|
|
|
616 } else {
|
|
|
617 while(nrow(patient.fuzzy) > 0){
|
|
|
618 first.merge = patient.fuzzy[1,"merge"]
|
|
|
619 first.clone.sequence = patient.fuzzy[1,"Clone_Sequence"]
|
|
|
620 first.sample = patient.fuzzy[1,"Sample"]
|
|
|
621
|
|
|
622 merge.filter = first.merge == patient.fuzzy$merge
|
|
|
623
|
|
|
624 second.sample = other.sample.list[[first.sample]][1]
|
|
|
625 third.sample = other.sample.list[[first.sample]][2]
|
|
|
626
|
|
|
627 sample.filter.1 = first.sample == patient.fuzzy$Sample
|
|
|
628 sample.filter.2 = second.sample == patient.fuzzy$Sample
|
|
|
629 sample.filter.3 = third.sample == patient.fuzzy$Sample
|
|
|
630
|
|
|
631 sequence.filter = grepl(paste("^", first.clone.sequence, sep=""), patient.fuzzy$Clone_Sequence)
|
|
|
632
|
|
|
633 match.filter.1 = sample.filter.1 & sequence.filter & merge.filter
|
|
|
634 match.filter.2 = sample.filter.2 & sequence.filter & merge.filter
|
|
|
635 match.filter.3 = sample.filter.3 & sequence.filter & merge.filter
|
|
|
636
|
|
|
637 matches.in.1 = any(match.filter.1)
|
|
|
638 matches.in.2 = any(match.filter.2)
|
|
|
639 matches.in.3 = any(match.filter.3)
|
|
|
640
|
|
|
641
|
|
|
642
|
|
|
643 rows.1 = patient.fuzzy[match.filter.1,]
|
|
|
644
|
|
|
645 sum.1 = data.frame(merge = first.clone.sequence,
|
|
|
646 Patient = label1,
|
|
|
647 Receptor = rows.1[1,"Receptor"],
|
|
|
648 Sample = rows.1[1,"Sample"],
|
|
|
649 Cell_Count = rows.1[1,"Cell_Count"],
|
|
|
650 Clone_Molecule_Count_From_Spikes = sum(rows.1$Clone_Molecule_Count_From_Spikes),
|
|
|
651 Log10_Frequency = log10(sum(rows.1$Frequency)),
|
|
|
652 Total_Read_Count = sum(rows.1$Total_Read_Count),
|
|
|
653 dsPerM = sum(rows.1$dsPerM),
|
|
|
654 J_Segment_Major_Gene = rows.1[1,"J_Segment_Major_Gene"],
|
|
|
655 V_Segment_Major_Gene = rows.1[1,"V_Segment_Major_Gene"],
|
|
|
656 Clone_Sequence = first.clone.sequence,
|
|
|
657 CDR3_Sense_Sequence = rows.1[1,"CDR3_Sense_Sequence"],
|
|
|
658 Related_to_leukemia_clone = F,
|
|
|
659 Frequency = sum(rows.1$Frequency),
|
|
|
660 locus_V = rows.1[1,"locus_V"],
|
|
|
661 locus_J = rows.1[1,"locus_J"],
|
|
|
662 uniqueID = rows.1[1,"uniqueID"],
|
|
|
663 normalized_read_count = sum(rows.1$normalized_read_count))
|
|
|
664 sum.2 = sum.1[NULL,]
|
|
|
665 rows.2 = patient.fuzzy[match.filter.2,]
|
|
|
666 if(matches.in.2){
|
|
|
667 sum.2 = data.frame(merge = first.clone.sequence,
|
|
|
668 Patient = label1,
|
|
|
669 Receptor = rows.2[1,"Receptor"],
|
|
|
670 Sample = rows.2[1,"Sample"],
|
|
|
671 Cell_Count = rows.2[1,"Cell_Count"],
|
|
|
672 Clone_Molecule_Count_From_Spikes = sum(rows.2$Clone_Molecule_Count_From_Spikes),
|
|
|
673 Log10_Frequency = log10(sum(rows.2$Frequency)),
|
|
|
674 Total_Read_Count = sum(rows.2$Total_Read_Count),
|
|
|
675 dsPerM = sum(rows.2$dsPerM),
|
|
|
676 J_Segment_Major_Gene = rows.2[1,"J_Segment_Major_Gene"],
|
|
|
677 V_Segment_Major_Gene = rows.2[1,"V_Segment_Major_Gene"],
|
|
|
678 Clone_Sequence = first.clone.sequence,
|
|
|
679 CDR3_Sense_Sequence = rows.2[1,"CDR3_Sense_Sequence"],
|
|
|
680 Related_to_leukemia_clone = F,
|
|
|
681 Frequency = sum(rows.2$Frequency),
|
|
|
682 locus_V = rows.2[1,"locus_V"],
|
|
|
683 locus_J = rows.2[1,"locus_J"],
|
|
|
684 uniqueID = rows.2[1,"uniqueID"],
|
|
|
685 normalized_read_count = sum(rows.2$normalized_read_count))
|
|
|
686 }
|
|
|
687
|
|
|
688 sum.3 = sum.1[NULL,]
|
|
|
689 rows.3 = patient.fuzzy[match.filter.3,]
|
|
|
690 if(matches.in.3){
|
|
|
691 sum.3 = data.frame(merge = first.clone.sequence,
|
|
|
692 Patient = label1,
|
|
|
693 Receptor = rows.3[1,"Receptor"],
|
|
|
694 Sample = rows.3[1,"Sample"],
|
|
|
695 Cell_Count = rows.3[1,"Cell_Count"],
|
|
|
696 Clone_Molecule_Count_From_Spikes = sum(rows.3$Clone_Molecule_Count_From_Spikes),
|
|
|
697 Log10_Frequency = log10(sum(rows.3$Frequency)),
|
|
|
698 Total_Read_Count = sum(rows.3$Total_Read_Count),
|
|
|
699 dsPerM = sum(rows.3$dsPerM),
|
|
|
700 J_Segment_Major_Gene = rows.3[1,"J_Segment_Major_Gene"],
|
|
|
701 V_Segment_Major_Gene = rows.3[1,"V_Segment_Major_Gene"],
|
|
|
702 Clone_Sequence = first.clone.sequence,
|
|
|
703 CDR3_Sense_Sequence = rows.3[1,"CDR3_Sense_Sequence"],
|
|
|
704 Related_to_leukemia_clone = F,
|
|
|
705 Frequency = sum(rows.3$Frequency),
|
|
|
706 locus_V = rows.3[1,"locus_V"],
|
|
|
707 locus_J = rows.3[1,"locus_J"],
|
|
|
708 uniqueID = rows.3[1,"uniqueID"],
|
|
|
709 normalized_read_count = sum(rows.3$normalized_read_count))
|
|
|
710 }
|
|
|
711
|
|
|
712 if(matches.in.2 & matches.in.3){
|
|
|
713 merge.123 = merge(sum.1, sum.2, by="merge")
|
|
|
714 merge.123 = merge(merge.123, sum.3, by="merge")
|
|
|
715 colnames(merge.123)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(merge.123)))] = paste(colnames(merge.123)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(merge.123), perl=T))], ".z", sep="")
|
|
|
716 #merge.123$thresholdValue = pmax(merge.123[,onx], merge.123[,ony], merge.123[,onz])
|
|
|
717
|
|
|
718 patientMerge = rbind(patientMerge, merge.123)
|
|
|
719 patient.fuzzy = patient.fuzzy[!(match.filter.1 | match.filter.2 | match.filter.3),]
|
|
|
720
|
|
|
721 hidden.clone.sequences = c(rows.1[-1,"Clone_Sequence"], rows.2[rows.2$Clone_Sequence != first.clone.sequence,"Clone_Sequence"], rows.3[rows.3$Clone_Sequence != first.clone.sequence,"Clone_Sequence"])
|
|
|
722 merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences)
|
|
|
723
|
|
|
724 } else if (matches.in.2) {
|
|
|
725 #other.sample1 = other.sample.list[[first.sample]][1]
|
|
|
726 #other.sample2 = other.sample.list[[first.sample]][2]
|
|
|
727
|
|
|
728 second.sample = sum.2[,"Sample"]
|
|
|
729
|
|
|
730 current.merge.list = duo.merge.list[[paste(first.sample, second.sample)]]
|
|
|
731
|
|
|
732 merge.12 = merge(sum.1, sum.2, by="merge")
|
|
|
733
|
|
|
734 current.merge.list = rbind(current.merge.list, merge.12)
|
|
|
735 duo.merge.list[[paste(first.sample, second.sample)]] = current.merge.list
|
|
|
736
|
|
|
737 patient.fuzzy = patient.fuzzy[!(match.filter.1 | match.filter.2),]
|
|
|
738
|
|
|
739 hidden.clone.sequences = c(rows.1[-1,"Clone_Sequence"], rows.2[rows.2$Clone_Sequence != first.clone.sequence,"Clone_Sequence"])
|
|
|
740 merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences)
|
|
|
741
|
|
|
742 } else if (matches.in.3) {
|
|
|
743
|
|
|
744 #other.sample1 = other.sample.list[[first.sample]][1]
|
|
|
745 #other.sample2 = other.sample.list[[first.sample]][2]
|
|
|
746
|
|
|
747 second.sample = sum.3[,"Sample"]
|
|
|
748
|
|
|
749 current.merge.list = duo.merge.list[[paste(first.sample, second.sample)]]
|
|
|
750
|
|
|
751 merge.13 = merge(sum.1, sum.3, by="merge")
|
|
|
752
|
|
|
753 current.merge.list = rbind(current.merge.list, merge.13)
|
|
|
754 duo.merge.list[[paste(first.sample, second.sample)]] = current.merge.list
|
|
|
755
|
|
|
756 patient.fuzzy = patient.fuzzy[!(match.filter.1 | match.filter.3),]
|
|
|
757
|
|
|
758 hidden.clone.sequences = c(rows.1[-1,"Clone_Sequence"], rows.3[rows.3$Clone_Sequence != first.clone.sequence,"Clone_Sequence"])
|
|
|
759 merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences)
|
|
|
760
|
|
|
761 } else if(nrow(rows.1) > 1){
|
|
|
762 patient1 = patient1[!(patient1$Clone_Sequence %in% rows.1$Clone_Sequence),]
|
|
|
763 print(names(patient1)[names(patient1) %in% sum.1])
|
|
|
764 print(names(patient1)[!(names(patient1) %in% sum.1)])
|
|
|
765 print(names(patient1))
|
|
|
766 print(names(sum.1))
|
|
|
767 print(summary(sum.1))
|
|
|
768 print(summary(patient1))
|
|
|
769 print(dim(sum.1))
|
|
|
770 print(dim(patient1))
|
|
|
771 print(head(sum.1[,names(patient1)]))
|
|
|
772 patient1 = rbind(patient1, sum.1[,names(patient1)])
|
|
|
773 patient.fuzzy = patient.fuzzy[-match.filter.1,]
|
|
|
774 } else {
|
|
|
775 patient.fuzzy = patient.fuzzy[-1,]
|
|
|
776 }
|
|
|
777
|
|
|
778 tmp.rows = rbind(rows.1, rows.2, rows.3)
|
|
|
779 tmp.rows = tmp.rows[order(nchar(tmp.rows$Clone_Sequence)),]
|
|
|
780
|
|
|
781 if (sum(match.filter.1) > 1 | sum(match.filter.2) > 1 | sum(match.filter.1) > 1) {
|
|
|
782 cat(paste("<tr><td>", label1, " row ", 1:nrow(tmp.rows), "</td><td>", tmp.rows$Sample, ":</td><td>", tmp.rows$Clone_Sequence, "</td><td>", tmp.rows$normalized_read_count, "</td></tr>", sep=""), file="multiple_matches.html", append=T)
|
|
|
783 } else {
|
|
|
784 }
|
|
|
785
|
|
|
786 }
|
|
|
787 patient.merge.list[[paste(label1, "123")]] = patientMerge
|
|
|
788
|
|
|
789 patientMerge12 = duo.merge.list[[paste(oneSample, twoSample)]]
|
|
|
790 patientMerge13 = duo.merge.list[[paste(oneSample, threeSample)]]
|
|
|
791 patientMerge23 = duo.merge.list[[paste(twoSample, threeSample)]]
|
|
|
792
|
|
|
793 patient.merge.list[[paste(label1, "12")]] = patientMerge12
|
|
|
794 patient.merge.list[[paste(label1, "13")]] = patientMerge13
|
|
|
795 patient.merge.list[[paste(label1, "23")]] = patientMerge23
|
|
|
796
|
|
|
797 patient.merge.list.second[[label1]] = merge.list[["second"]]
|
|
|
798 }
|
|
|
799 cat(paste("<td>", nrow(patient1), " in ", label1, " and ", nrow(patient2), " in ", label2, nrow(patient3), " in ", label3, ", ", nrow(patientMerge), " in both (finding both took ", (proc.time() - start.time)[[3]], "s)</td></tr>", sep=""), file=logfile, append=T)
|
|
|
800 patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz])
|
|
|
801 patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony])
|
|
|
802 patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony])
|
|
|
803 patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony])
|
|
|
804
|
|
|
805 #patientMerge$thresholdValue = pmin(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz])
|
|
|
806 #patientMerge12$thresholdValue = pmin(patientMerge12[,onx], patientMerge12[,ony])
|
|
|
807 #patientMerge13$thresholdValue = pmin(patientMerge13[,onx], patientMerge13[,ony])
|
|
|
808 #patientMerge23$thresholdValue = pmin(patientMerge23[,onx], patientMerge23[,ony])
|
|
|
809
|
|
|
810 patient1 = patient1[!(patient1$Clone_Sequence %in% merge.list[["second"]]),]
|
|
|
811 patient2 = patient2[!(patient2$Clone_Sequence %in% merge.list[["second"]]),]
|
|
|
812 patient3 = patient3[!(patient3$Clone_Sequence %in% merge.list[["second"]]),]
|
|
|
813
|
|
|
814 if(F){
|
|
|
815 patientMerge = merge(patient1, patient2, by="merge")
|
|
|
816 patientMerge = merge(patientMerge, patient3, by="merge")
|
|
|
817 colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge)))] = paste(colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge), perl=T))], ".z", sep="")
|
|
|
818 patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz])
|
|
|
819 patientMerge12 = merge(patient1, patient2, by="merge")
|
|
|
820 patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony])
|
|
|
821 patientMerge13 = merge(patient1, patient3, by="merge")
|
|
|
822 patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony])
|
|
|
823 patientMerge23 = merge(patient2, patient3, by="merge")
|
|
|
824 patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony])
|
|
|
825 }
|
|
|
826
|
|
|
827 scatterplot_data_columns = c("Clone_Sequence", "Frequency", "normalized_read_count", "V_Segment_Major_Gene", "J_Segment_Major_Gene", "merge")
|
|
|
828 scatterplot_data = rbind(patient1[,scatterplot_data_columns], patient2[,scatterplot_data_columns], patient3[,scatterplot_data_columns])
|
|
|
829 scatterplot_data = scatterplot_data[!duplicated(scatterplot_data$merge),]
|
|
|
830 scatterplot_data$type = factor(x="In one", levels=c("In one", "In two", "In three", "In multiple"))
|
|
|
831
|
|
|
832 res1 = vector()
|
|
|
833 res2 = vector()
|
|
|
834 res3 = vector()
|
|
|
835 res12 = vector()
|
|
|
836 res13 = vector()
|
|
|
837 res23 = vector()
|
|
|
838 resAll = vector()
|
|
|
839 read1Count = vector()
|
|
|
840 read2Count = vector()
|
|
|
841 read3Count = vector()
|
|
|
842
|
|
|
843 if(appendTriplets){
|
|
|
844 cat(paste(label1, label2, label3, sep="\t"), file="triplets.txt", append=T, sep="", fill=3)
|
|
|
845 }
|
|
|
846 for(iter in 1:length(product[,1])){
|
|
|
847 threshhold = product[iter,threshholdIndex]
|
|
|
848 V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="")
|
|
|
849 J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="")
|
|
|
850 #all = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge[,onx] > threshhold & patientMerge[,ony] > threshhold & patientMerge[,onz] > threshhold)
|
|
|
851 all = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge$thresholdValue > threshhold)
|
|
|
852
|
|
|
853 one_two = (grepl(V_Segment, patientMerge12$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge12$J_Segment_Major_Gene.x) & patientMerge12$thresholdValue > threshhold & !(patientMerge12$merge %in% patientMerge[all,]$merge))
|
|
|
854 one_three = (grepl(V_Segment, patientMerge13$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge13$J_Segment_Major_Gene.x) & patientMerge13$thresholdValue > threshhold & !(patientMerge13$merge %in% patientMerge[all,]$merge))
|
|
|
855 two_three = (grepl(V_Segment, patientMerge23$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge23$J_Segment_Major_Gene.x) & patientMerge23$thresholdValue > threshhold & !(patientMerge23$merge %in% patientMerge[all,]$merge))
|
|
|
856
|
|
|
857 one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$merge %in% patientMerge[all,]$merge) & !(patient1$merge %in% patientMerge12[one_two,]$merge) & !(patient1$merge %in% patientMerge13[one_three,]$merge))
|
|
|
858 two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$merge %in% patientMerge[all,]$merge) & !(patient2$merge %in% patientMerge12[one_two,]$merge) & !(patient2$merge %in% patientMerge23[two_three,]$merge))
|
|
|
859 three = (grepl(V_Segment, patient3$V_Segment_Major_Gene) & grepl(J_Segment, patient3$J_Segment_Major_Gene) & patient3[,on] > threshhold & !(patient3$merge %in% patientMerge[all,]$merge) & !(patient3$merge %in% patientMerge13[one_three,]$merge) & !(patient3$merge %in% patientMerge23[two_three,]$merge))
|
|
|
860
|
|
|
861 read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.x))
|
|
|
862 read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.y))
|
|
|
863 read3Count = append(read3Count, sum(patient3[three,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.z))
|
|
|
864 res1 = append(res1, sum(one))
|
|
|
865 res2 = append(res2, sum(two))
|
|
|
866 res3 = append(res3, sum(three))
|
|
|
867 resAll = append(resAll, sum(all))
|
|
|
868 res12 = append(res12, sum(one_two))
|
|
|
869 res13 = append(res13, sum(one_three))
|
|
|
870 res23 = append(res23, sum(two_three))
|
|
|
871 #threshhold = 0
|
|
|
872 if(threshhold != 0){
|
|
|
873 if(sum(one) > 0){
|
|
|
874 dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
|
|
|
875 colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")
|
|
|
876 filenameOne = paste(label1, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
877 write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
878 }
|
|
|
879 if(sum(two) > 0){
|
|
|
880 dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
|
|
|
881 colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")
|
|
|
882 filenameTwo = paste(label2, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
883 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
884 }
|
|
|
885 if(sum(three) > 0){
|
|
|
886 dfThree = patient3[three,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
|
|
|
887 colnames(dfThree) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")
|
|
|
888 filenameThree = paste(label3, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
889 write.table(dfThree, file=paste(filenameThree, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
890 }
|
|
|
891 if(sum(one_two) > 0){
|
|
|
892 dfOne_two = patientMerge12[one_two,c("V_Segment_Major_Gene.x", "J_Segment_Major_Gene.x", "normalized_read_count.x", "Frequency.x", "Related_to_leukemia_clone.x", "Clone_Sequence.x", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "normalized_read_count.y", "Frequency.y", "Related_to_leukemia_clone.y")]
|
|
|
893 colnames(dfOne_two) = c(paste("Proximal segment", oneSample), paste("Distal segment", oneSample), paste("Normalized_Read_Count", oneSample), paste("Frequency", oneSample), paste("Related_to_leukemia_clone", oneSample),"Clone_Sequence", paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Normalized_Read_Count", twoSample), paste("Frequency", twoSample), paste("Related_to_leukemia_clone", twoSample))
|
|
|
894 filenameOne_two = paste(label1, "_", label2, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="")
|
|
|
895 write.table(dfOne_two, file=paste(filenameOne_two, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
896 }
|
|
|
897 if(sum(one_three) > 0){
|
|
|
898 dfOne_three = patientMerge13[one_three,c("V_Segment_Major_Gene.x", "J_Segment_Major_Gene.x", "normalized_read_count.x", "Frequency.x", "Related_to_leukemia_clone.x", "Clone_Sequence.x", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "normalized_read_count.y", "Frequency.y", "Related_to_leukemia_clone.y")]
|
|
|
899 colnames(dfOne_three) = c(paste("Proximal segment", oneSample), paste("Distal segment", oneSample), paste("Normalized_Read_Count", oneSample), paste("Frequency", oneSample), paste("Related_to_leukemia_clone", oneSample),"Clone_Sequence", paste("Proximal segment", threeSample), paste("Distal segment", threeSample), paste("Normalized_Read_Count", threeSample), paste("Frequency", threeSample), paste("Related_to_leukemia_clone", threeSample))
|
|
|
900 filenameOne_three = paste(label1, "_", label3, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="")
|
|
|
901 write.table(dfOne_three, file=paste(filenameOne_three, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
902 }
|
|
|
903 if(sum(two_three) > 0){
|
|
|
904 dfTwo_three = patientMerge23[two_three,c("V_Segment_Major_Gene.x", "J_Segment_Major_Gene.x", "normalized_read_count.x", "Frequency.x", "Related_to_leukemia_clone.x", "Clone_Sequence.x", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "normalized_read_count.y", "Frequency.y", "Related_to_leukemia_clone.y")]
|
|
|
905 colnames(dfTwo_three) = c(paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Normalized_Read_Count", twoSample), paste("Frequency", twoSample), paste("Related_to_leukemia_clone", twoSample),"Clone_Sequence", paste("Proximal segment", threeSample), paste("Distal segment", threeSample), paste("Normalized_Read_Count", threeSample), paste("Frequency", threeSample), paste("Related_to_leukemia_clone", threeSample))
|
|
|
906 filenameTwo_three = paste(label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="")
|
|
|
907 write.table(dfTwo_three, file=paste(filenameTwo_three, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
908 }
|
|
|
909 } else { #scatterplot data
|
|
|
910 scatterplot_locus_data = scatterplot_data[grepl(V_Segment, scatterplot_data$V_Segment_Major_Gene) & grepl(J_Segment, scatterplot_data$J_Segment_Major_Gene),]
|
|
|
911 scatterplot_locus_data = scatterplot_locus_data[!(scatterplot_locus_data$merge %in% merge.list[["second"]]),]
|
|
|
912 in_two = (scatterplot_locus_data$merge %in% patientMerge12[one_two,]$merge) | (scatterplot_locus_data$merge %in% patientMerge13[one_three,]$merge) | (scatterplot_locus_data$merge %in% patientMerge23[two_three,]$merge)
|
|
|
913 if(sum(in_two) > 0){
|
|
|
914 scatterplot_locus_data[in_two,]$type = "In two"
|
|
|
915 }
|
|
|
916 in_three = (scatterplot_locus_data$merge %in% patientMerge[all,]$merge)
|
|
|
917 if(sum(in_three)> 0){
|
|
|
918 scatterplot_locus_data[in_three,]$type = "In three"
|
|
|
919 }
|
|
|
920 not_in_one = scatterplot_locus_data$type != "In one"
|
|
|
921 if(sum(not_in_one) > 0){
|
|
|
922 #scatterplot_locus_data[not_in_one,]$type = "In multiple"
|
|
|
923 }
|
|
|
924 p = NULL
|
|
|
925 if(nrow(scatterplot_locus_data) != 0){
|
|
|
926 if(on == "normalized_read_count"){
|
|
|
927 scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count))))
|
|
|
928 p = ggplot(scatterplot_locus_data, aes(type, normalized_read_count)) + scale_y_log10(breaks=scales,labels=scales, limits=c(1, 1e6))
|
|
|
929 } else {
|
|
|
930 p = ggplot(scatterplot_locus_data, aes(type, Frequency)) + scale_y_log10(limits=c(0.0001,100), breaks=c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100), labels=c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100")) + expand_limits(y=c(0,100))
|
|
|
931 #p = ggplot(scatterplot_locus_data, aes(type, Frequency)) + scale_y_continuous(limits = c(0, 100)) + expand_limits(y=c(0,100))
|
|
|
932 }
|
|
|
933 p = p + geom_point(aes(colour=type), position="jitter")
|
|
|
934 p = p + xlab("In one or in multiple samples") + ylab(onShort) + ggtitle(paste(label1, label2, label3, onShort, product[iter, titleIndex]))
|
|
|
935 } else {
|
|
|
936 p = ggplot(NULL, aes(x=c("In one", "In multiple"),y=0)) + geom_blank(NULL) + xlab("In two or in three of the samples") + ylab(onShort) + ggtitle(paste(label1, label2, label3, onShort, product[iter, titleIndex]))
|
|
|
937 }
|
|
|
938 png(paste(label1, "_", label2, "_", label3, "_", onShort, "_", product[iter, titleIndex],"_scatter.png", sep=""))
|
|
|
939 print(p)
|
|
|
940 dev.off()
|
|
|
941 }
|
|
|
942 if(sum(all) > 0){
|
|
|
943 dfAll = patientMerge[all,c("V_Segment_Major_Gene.x", "J_Segment_Major_Gene.x", "normalized_read_count.x", "Frequency.x", "Related_to_leukemia_clone.x", "Clone_Sequence.x", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "normalized_read_count.y", "Frequency.y", "Related_to_leukemia_clone.y", "V_Segment_Major_Gene.z", "J_Segment_Major_Gene.z", "normalized_read_count.z", "Frequency.z", "Related_to_leukemia_clone.z")]
|
|
|
944 colnames(dfAll) = c(paste("Proximal segment", oneSample), paste("Distal segment", oneSample), paste("Normalized_Read_Count", oneSample), paste("Frequency", oneSample), paste("Related_to_leukemia_clone", oneSample),"Clone_Sequence", paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Normalized_Read_Count", twoSample), paste("Frequency", twoSample), paste("Related_to_leukemia_clone", twoSample), paste("Proximal segment", threeSample), paste("Distal segment", threeSample), paste("Normalized_Read_Count", threeSample), paste("Frequency", threeSample), paste("Related_to_leukemia_clone", threeSample))
|
|
|
945 filenameAll = paste(label1, "_", label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
946 write.table(dfAll, file=paste(filenameAll, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
947 }
|
|
|
948 }
|
|
|
949 #patientResult = data.frame("Locus"=product$Titles, "J_Segment"=product$J_Segments, "V_Segment"=product$V_Segments, "cut_off_value"=paste(">", product$interval, sep=""), "All"=resAll, "tmp1"=res1, "read_count1" = round(read1Count), "tmp2"=res2, "read_count2"= round(read2Count), "tmp3"=res3, "read_count3"=round(read3Count))
|
|
|
950 patientResult = data.frame("Locus"=product$Titles, "J_Segment"=product$J_Segments, "V_Segment"=product$V_Segments, "cut_off_value"=paste(">", product$interval, sep=""), "All"=resAll, "tmp1"=res1, "tmp2"=res2, "tmp3"=res3, "tmp12"=res12, "tmp13"=res13, "tmp23"=res23)
|
|
|
951 colnames(patientResult)[6] = oneSample
|
|
|
952 colnames(patientResult)[7] = twoSample
|
|
|
953 colnames(patientResult)[8] = threeSample
|
|
|
954 colnames(patientResult)[9] = paste(oneSample, twoSample, sep="_")
|
|
|
955 colnames(patientResult)[10] = paste(oneSample, twoSample, sep="_")
|
|
|
956 colnames(patientResult)[11] = paste(oneSample, twoSample, sep="_")
|
|
|
957
|
|
|
958 colnamesBak = colnames(patientResult)
|
|
|
959 colnames(patientResult) = c("Ig/TCR gene rearrangement type", "Distal Gene segment", "Proximal gene segment", "cut_off_value", "Number of sequences All", paste("Number of sequences", oneSample), paste("Number of sequences", twoSample), paste("Number of sequences", threeSample), paste("Number of sequences", oneSample, twoSample), paste("Number of sequences", oneSample, threeSample), paste("Number of sequences", twoSample, threeSample))
|
|
|
960 write.table(patientResult, file=paste(label1, "_", label2, "_", label3, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
961 colnames(patientResult) = colnamesBak
|
|
|
962
|
|
|
963 patientResult$Locus = factor(patientResult$Locus, Titles)
|
|
|
964 patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))
|
|
|
965
|
|
|
966 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "All")])
|
|
|
967 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=All), stat='identity', position="dodge", fill="#79c36a")
|
|
|
968 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
969 plt = plt + geom_text(aes(ymax=max(All), x=cut_off_value,y=All,label=All), angle=90, hjust=0)
|
|
|
970 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in All")
|
|
|
971 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
|
972 png(paste(label1, "_", label2, "_", label3, "_", onShort, "_total_all.png", sep=""), width=1920, height=1080)
|
|
|
973 print(plt)
|
|
|
974 dev.off()
|
|
|
975
|
|
|
976 fontSize = 4
|
|
|
977
|
|
|
978 bak = patientResult
|
|
|
979 patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample, threeSample)] ,id.vars=1:2)
|
|
|
980 patientResult$relativeValue = patientResult$value * 10
|
|
|
981 patientResult[patientResult$relativeValue == 0,]$relativeValue = 1
|
|
|
982 plt = ggplot(patientResult)
|
|
|
983 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge")
|
|
|
984 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
985 plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9)))
|
|
|
986 plt = plt + geom_text(data=patientResult[patientResult$variable == oneSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=-0.7, size=fontSize)
|
|
|
987 plt = plt + geom_text(data=patientResult[patientResult$variable == twoSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=0.4, size=fontSize)
|
|
|
988 plt = plt + geom_text(data=patientResult[patientResult$variable == threeSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=1.5, size=fontSize)
|
|
|
989 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in only one sample")
|
|
|
990 png(paste(label1, "_", label2, "_", label3, "_", onShort, "_indiv_all.png", sep=""), width=1920, height=1080)
|
|
|
991 print(plt)
|
|
|
992 dev.off()
|
|
|
993 }
|
|
|
994
|
|
|
995 if(nrow(triplets) != 0){
|
|
|
996
|
|
|
997 cat("<tr><td>Starting triplet analysis</td></tr>", file=logfile, append=T)
|
|
|
998
|
|
|
999 triplets$uniqueID = "ID"
|
|
|
1000
|
|
|
1001 triplets[grepl("16278_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
|
|
|
1002 triplets[grepl("26402_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
|
|
|
1003 triplets[grepl("26759_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
|
|
|
1004
|
|
|
1005 triplets[grepl("16278_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
|
|
|
1006 triplets[grepl("26402_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
|
|
|
1007 triplets[grepl("26759_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
|
|
|
1008
|
|
|
1009 triplets[grepl("14696", triplets$Patient),]$uniqueID = "14696"
|
|
|
1010
|
|
|
1011 cat("<tr><td>Normalizing to lowest cell count within locus</td></tr>", file=logfile, append=T)
|
|
|
1012
|
|
|
1013 triplets$locus_V = substring(triplets$V_Segment_Major_Gene, 0, 4)
|
|
|
1014 triplets$locus_J = substring(triplets$J_Segment_Major_Gene, 0, 4)
|
|
|
1015 min_cell_count = data.frame(data.table(triplets)[, list(min_cell_count=min(.SD$Cell_Count)), by=c("uniqueID", "locus_V", "locus_J")])
|
|
|
1016
|
|
|
1017 triplets$min_cell_paste = paste(triplets$uniqueID, triplets$locus_V, triplets$locus_J)
|
|
|
1018 min_cell_count$min_cell_paste = paste(min_cell_count$uniqueID, min_cell_count$locus_V, min_cell_count$locus_J)
|
|
|
1019
|
|
|
1020 min_cell_count = min_cell_count[,c("min_cell_paste", "min_cell_count")]
|
|
|
1021
|
|
|
1022 triplets = merge(triplets, min_cell_count, by="min_cell_paste")
|
|
|
1023
|
|
|
1024 triplets$normalized_read_count = round(triplets$Clone_Molecule_Count_From_Spikes / triplets$Cell_Count * triplets$min_cell_count / 2, digits=2) #??????????????????????????????????? wel of geen / 2
|
|
|
1025
|
|
|
1026 triplets = triplets[triplets$normalized_read_count >= min_cells,]
|
|
|
1027
|
|
|
1028 column_drops = c("min_cell_count", "min_cell_paste")
|
|
|
1029
|
|
|
1030 triplets = triplets[,!(colnames(triplets) %in% column_drops)]
|
|
|
1031
|
|
|
1032 cat("<tr><td>Starting Cell Count analysis</td></tr>", file=logfile, append=T)
|
|
|
1033
|
|
|
1034 interval = intervalReads
|
|
|
1035 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
|
1036 product = data.frame("Titles"=rep(Titles, each=length(interval)), "interval"=rep(interval, times=10), "V_Segments"=rep(V_Segments, each=length(interval)), "J_Segments"=rep(J_Segments, each=length(interval)))
|
|
|
1037
|
|
|
1038 one = triplets[triplets$Sample == "14696_reg_BM",]
|
|
|
1039 two = triplets[triplets$Sample == "24536_reg_BM",]
|
|
|
1040 three = triplets[triplets$Sample == "24062_reg_BM",]
|
|
|
1041 tripletAnalysis(one, "14696_1_Trio", two, "14696_2_Trio", three, "14696_3_Trio", product=product, interval=interval, on="normalized_read_count", T)
|
|
|
1042
|
|
|
1043 one = triplets[triplets$Sample == "16278_Left",]
|
|
|
1044 two = triplets[triplets$Sample == "26402_Left",]
|
|
|
1045 three = triplets[triplets$Sample == "26759_Left",]
|
|
|
1046 tripletAnalysis(one, "16278_Left_Trio", two, "26402_Left_Trio", three, "26759_Left_Trio", product=product, interval=interval, on="normalized_read_count", T)
|
|
|
1047
|
|
|
1048 one = triplets[triplets$Sample == "16278_Right",]
|
|
|
1049 two = triplets[triplets$Sample == "26402_Right",]
|
|
|
1050 three = triplets[triplets$Sample == "26759_Right",]
|
|
|
1051 tripletAnalysis(one, "16278_Right_Trio", two, "26402_Right_Trio", three, "26759_Right_Trio", product=product, interval=interval, on="normalized_read_count", T)
|
|
|
1052
|
|
|
1053 cat("<tr><td>Starting Frequency analysis</td></tr>", file=logfile, append=T)
|
|
|
1054
|
|
|
1055 interval = intervalFreq
|
|
|
1056 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
|
1057 product = data.frame("Titles"=rep(Titles, each=length(interval)), "interval"=rep(interval, times=10), "V_Segments"=rep(V_Segments, each=length(interval)), "J_Segments"=rep(J_Segments, each=length(interval)))
|
|
|
1058
|
|
|
1059 one = triplets[triplets$Sample == "14696_reg_BM",]
|
|
|
1060 two = triplets[triplets$Sample == "24536_reg_BM",]
|
|
|
1061 three = triplets[triplets$Sample == "24062_reg_BM",]
|
|
|
1062 tripletAnalysis(one, "14696_1_Trio", two, "14696_2_Trio", three, "14696_3_Trio", product=product, interval=interval, on="Frequency", F)
|
|
|
1063
|
|
|
1064 one = triplets[triplets$Sample == "16278_Left",]
|
|
|
1065 two = triplets[triplets$Sample == "26402_Left",]
|
|
|
1066 three = triplets[triplets$Sample == "26759_Left",]
|
|
|
1067 tripletAnalysis(one, "16278_Left_Trio", two, "26402_Left_Trio", three, "26759_Left_Trio", product=product, interval=interval, on="Frequency", F)
|
|
|
1068
|
|
|
1069 one = triplets[triplets$Sample == "16278_Right",]
|
|
|
1070 two = triplets[triplets$Sample == "26402_Right",]
|
|
|
1071 three = triplets[triplets$Sample == "26759_Right",]
|
|
|
1072 tripletAnalysis(one, "16278_Right_Trio", two, "26402_Right_Trio", three, "26759_Right_Trio", product=product, interval=interval, on="Frequency", F)
|
|
|
1073 } else {
|
|
|
1074 cat("", file="triplets.txt")
|
|
|
1075 }
|
|
|
1076 cat("</table></html>", file=logfile, append=T)
|
