Mercurial > repos > davidvanzessen > clonal_sequences_in_paired_samples

comparison RScript.r @ 49:7658e9f3d416 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Thu, 08 Oct 2015 05:46:16 -0400
parents 1b5b862b055b
children 7dd7cefcf72d
comparison
equal deleted inserted replaced
48:1b5b862b055b 49:7658e9f3d416
52 52
53 dat = dat[dat$normalized_read_count >= min_cells,] 53 dat = dat[dat$normalized_read_count >= min_cells,]
54 54
55 dat$paste = paste(dat$Sample, dat$Clone_Sequence) 55 dat$paste = paste(dat$Sample, dat$Clone_Sequence)
56 56
57 #remove duplicate V+J+CDR3, add together numerical values
58 if(mergeOn != "Clone_Sequence"){
59 cat("<tr><td>Adding duplicate V+J+CDR3 sequences</td></tr>", file=logfile, append=T)
60 dat= data.frame(data.table(dat)[, list(Receptor=unique(.SD$Receptor),
61 Cell_Count=unique(.SD$Cell_Count),
62 Clone_Molecule_Count_From_Spikes=sum(.SD$Clone_Molecule_Count_From_Spikes),
63 Total_Read_Count=sum(.SD$Total_Read_Count),
64 dsPerM=ifelse("dsPerM" %in% names(dat), sum(.SD$dsPerM), 0),
65 Related_to_leukemia_clone=all(.SD$Related_to_leukemia_clone),
66 Frequency=sum(.SD$Frequency),
67 locus_V=unique(.SD$locus_V),
68 locus_J=unique(.SD$locus_J),
69 min_cell_count=unique(.SD$min_cell_count),
70 normalized_read_count=sum(.SD$normalized_read_count),
71 Log10_Frequency=sum(.SD$Log10_Frequency),
72 Clone_Sequence=.SD$Clone_Sequence[1],
73 min_cell_paste=.SD$min_cell_paste[1],
74 paste=unique(.SD$paste)), by=c("Patient", "Sample", "V_Segment_Major_Gene", "J_Segment_Major_Gene", "CDR3_Sense_Sequence")])
75 }
76
77 patients = split(dat, dat$Patient, drop=T) 57 patients = split(dat, dat$Patient, drop=T)
78 intervalReads = rev(c(0,10,25,50,100,250,500,750,1000,10000)) 58 intervalReads = rev(c(0,10,25,50,100,250,500,750,1000,10000))
79 intervalFreq = rev(c(0,0.01,0.05,0.1,0.5,1,5)) 59 intervalFreq = rev(c(0,0.01,0.05,0.1,0.5,1,5))
80 V_Segments = c(".*", "IGHV", "IGHD", "IGKV", "IGKV", "IgKINTR", "TRGV", "TRDV", "TRDD" , "TRBV") 60 V_Segments = c(".*", "IGHV", "IGHD", "IGKV", "IGKV", "IgKINTR", "TRGV", "TRDV", "TRDD" , "TRBV")
81 J_Segments = c(".*", ".*", ".*", "IGKJ", "KDE", ".*", ".*", ".*", ".*", ".*") 61 J_Segments = c(".*", ".*", ".*", "IGKJ", "KDE", ".*", ".*", ".*", ".*", ".*")
144 scatterplot_data = rbind(patient1[,scatterplot_data_columns], patient2[,scatterplot_data_columns]) 124 scatterplot_data = rbind(patient1[,scatterplot_data_columns], patient2[,scatterplot_data_columns])
145 scatterplot_data = scatterplot_data[!duplicated(scatterplot_data$merge),] 125 scatterplot_data = scatterplot_data[!duplicated(scatterplot_data$merge),]
146 scatterplot_data$type = factor(x=oneSample, levels=c(oneSample, twoSample, "In Both")) 126 scatterplot_data$type = factor(x=oneSample, levels=c(oneSample, twoSample, "In Both"))
147 scatterplot_data$on = onShort 127 scatterplot_data$on = onShort
148 128
149 patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge") 129 #patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge") #merge alles 'fuzzy'
150 130 patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge")[NULL,] #blegh
131
132 cs.exact.matches = patient1[patient1$Clone_Sequence %in% patient2$Clone_Sequence,]$Clone_Sequence
133
151 134
152 #fuzzy matching here... 135 #fuzzy matching here...
153 if(mergeOn == "Clone_Sequence"){ 136 if(mergeOn == "Clone_Sequence"){
154 merge.list = patientMerge$merge 137 #merge.list = patientMerge$merge
155 138
156 patient1.fuzzy = patient1[!(patient1$merge %in% merge.list),] 139 #patient1.fuzzy = patient1[!(patient1$merge %in% merge.list),]
157 patient2.fuzzy = patient2[!(patient2$merge %in% merge.list),] 140 #patient2.fuzzy = patient2[!(patient2$merge %in% merge.list),]
158 141
142 patient1.fuzzy = patient1
143 patient2.fuzzy = patient2
144
159 #patient1.fuzzy$merge = paste(patient1.fuzzy$V_Segment_Major_Gene, patient1.fuzzy$J_Segment_Major_Gene, patient1.fuzzy$CDR3_Sense_Sequence) 145 #patient1.fuzzy$merge = paste(patient1.fuzzy$V_Segment_Major_Gene, patient1.fuzzy$J_Segment_Major_Gene, patient1.fuzzy$CDR3_Sense_Sequence)
160 #patient2.fuzzy$merge = paste(patient2.fuzzy$V_Segment_Major_Gene, patient2.fuzzy$J_Segment_Major_Gene, patient2.fuzzy$CDR3_Sense_Sequence) 146 #patient2.fuzzy$merge = paste(patient2.fuzzy$V_Segment_Major_Gene, patient2.fuzzy$J_Segment_Major_Gene, patient2.fuzzy$CDR3_Sense_Sequence)
161 147
162 #patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J, patient1.fuzzy$CDR3_Sense_Sequence) 148 #patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J, patient1.fuzzy$CDR3_Sense_Sequence)
163 #patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J, patient2.fuzzy$CDR3_Sense_Sequence) 149 #patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J, patient2.fuzzy$CDR3_Sense_Sequence)
164 150
165 patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J) 151 patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J)
166 patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J) 152 patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J)
167 153
168 merge.freq.table = data.frame(table(c(patient1.fuzzy[!duplicated(patient1.fuzzy$merge),"merge"], patient2.fuzzy[!duplicated(patient2.fuzzy$merge),"merge"]))) 154 #merge.freq.table = data.frame(table(c(patient1.fuzzy[!duplicated(patient1.fuzzy$merge),"merge"], patient2.fuzzy[!duplicated(patient2.fuzzy$merge),"merge"]))) #also remove?
169 merge.freq.table.gt.1 = merge.freq.table[merge.freq.table$Freq > 1,] 155 #merge.freq.table.gt.1 = merge.freq.table[merge.freq.table$Freq > 1,]
170 156
171 patient1.fuzzy = patient1.fuzzy[patient1.fuzzy$merge %in% merge.freq.table.gt.1$Var1,] 157 #patient1.fuzzy = patient1.fuzzy[patient1.fuzzy$merge %in% merge.freq.table.gt.1$Var1,]
172 patient2.fuzzy = patient2.fuzzy[patient2.fuzzy$merge %in% merge.freq.table.gt.1$Var1,] 158 #patient2.fuzzy = patient2.fuzzy[patient2.fuzzy$merge %in% merge.freq.table.gt.1$Var1,]
173 159
174 patient.fuzzy = rbind(patient1.fuzzy, patient2.fuzzy) 160 patient.fuzzy = rbind(patient1.fuzzy, patient2.fuzzy)
175 patient.fuzzy = patient.fuzzy[order(nchar(patient.fuzzy$Clone_Sequence)),] 161 patient.fuzzy = patient.fuzzy[order(nchar(patient.fuzzy$Clone_Sequence)),]
176 162
163 merge.list = list()
164
165 merge.list[["second"]] = vector()
166
167
177 while(nrow(patient.fuzzy) > 1){ 168 while(nrow(patient.fuzzy) > 1){
178 first.merge = patient.fuzzy[1,"merge"] 169 first.merge = patient.fuzzy[1,"merge"]
179 first.clone.sequence = patient.fuzzy[1,"Clone_Sequence"] 170 first.clone.sequence = patient.fuzzy[1,"Clone_Sequence"]
180 171 first.sample = patient.fuzzy[1,"Sample"]
181 merge.filter = first.merge == patient.fuzzy$merge 172 merge.filter = first.merge == patient.fuzzy$merge
182 173
183 length.filter = nchar(patient.fuzzy$Clone_Sequence) - nchar(first.clone.sequence) <= 9 174 #length.filter = nchar(patient.fuzzy$Clone_Sequence) - nchar(first.clone.sequence) <= 9
184 175
185 sample.filter = patient.fuzzy[1,"Sample"] != patient.fuzzy$Sample 176 first.sample.filter = first.sample == patient.fuzzy$Sample
177 second.sample.filter = first.sample != patient.fuzzy$Sample
178
179 #first match same sample, sum to a single row, same for other sample
180 #then merge rows like 'normal'
186 181
187 sequence.filter = grepl(paste("^", first.clone.sequence, sep=""), patient.fuzzy$Clone_Sequence) 182 sequence.filter = grepl(paste("^", first.clone.sequence, sep=""), patient.fuzzy$Clone_Sequence)
188 183
184
185
189 #match.filter = merge.filter & grepl(first.clone.sequence, patient.fuzzy$Clone_Sequence) & length.filter & sample.filter 186 #match.filter = merge.filter & grepl(first.clone.sequence, patient.fuzzy$Clone_Sequence) & length.filter & sample.filter
190 match.filter = merge.filter & sequence.filter & sample.filter 187 first.match.filter = merge.filter & sequence.filter & first.sample.filter
191 188 second.match.filter = merge.filter & sequence.filter & second.sample.filter
192 if(sum(match.filter) == 1){ 189
193 second.match = which(match.filter)[1] 190 first.rows = patient.fuzzy[first.match.filter,]
194 second.clone.sequence = patient.fuzzy[second.match,"Clone_Sequence"] 191 second.rows = patient.fuzzy[second.match.filter,]
195 first.sample = patient.fuzzy[1,"Sample"] 192
196 second.sample = patient.fuzzy[second.match,"Sample"] 193 first.sum = data.frame(merge = first.clone.sequence,
197 194 Patient = patient,
198 first.match.row = patient.fuzzy[1,] 195 Receptor = first.rows[1,"Receptor"],
199 second.match.row = patient.fuzzy[second.match,] 196 Sample = first.rows[1,"Sample"],
200 print(paste(first.merge, first.match.row$normalized_read_count, second.match.row$normalized_read_count, first.clone.sequence, second.clone.sequence)) 197 Cell_Count = first.rows[1,"Cell_Count"],
201 patientMerge.new.row = data.frame(merge=first.clone.sequence, 198 Clone_Molecule_Count_From_Spikes = sum(first.rows$Clone_Molecule_Count_From_Spikes),
202 min_cell_paste.x=first.match.row[1,"min_cell_paste"], 199 Log10_Frequency = log10(sum(first.rows$Frequency)),
203 Patient.x=first.match.row[1,"Patient"], 200 Total_Read_Count = sum(first.rows$Total_Read_Count),
204 Receptor.x=first.match.row[1,"Receptor"], 201 dsPerM = sum(first.rows$dsPerM),
205 Sample.x=first.match.row[1,"Sample"], 202 J_Segment_Major_Gene = sort(table(first.rows$J_Segment_Major_Gene),decreasing=TRUE)[1],
206 Cell_Count.x=first.match.row[1,"Cell_Count"], 203 V_Segment_Major_Gene = sort(table(first.rows$V_Segment_Major_Gene),decreasing=TRUE)[1],
207 Clone_Molecule_Count_From_Spikes.x=first.match.row[1,"Clone_Molecule_Count_From_Spikes"], 204 Clone_Sequence = first.clone.sequence,
208 Log10_Frequency.x=first.match.row[1,"Log10_Frequency"], 205 CDR3_Sense_Sequence = first.rows[1,"CDR3_Sense_Sequence"],
209 Total_Read_Count.x=first.match.row[1,"Total_Read_Count"], 206 Related_to_leukemia_clone = F,
210 dsPerM.x=first.match.row[1,"dsPerM"], 207 Frequency = sum(first.rows$Frequency),
211 J_Segment_Major_Gene.x=first.match.row[1,"J_Segment_Major_Gene"], 208 locus_V = first.rows[1,"locus_V"],
212 V_Segment_Major_Gene.x=first.match.row[1,"V_Segment_Major_Gene"], 209 locus_J = first.rows[1,"locus_J"],
213 Clone_Sequence.x=first.match.row[1,"Clone_Sequence"], 210 min_cell_count = first.rows[1,"min_cell_count"],
214 CDR3_Sense_Sequence.x=first.match.row[1,"CDR3_Sense_Sequence"], 211 normalized_read_count = sum(first.rows$normalized_read_count),
215 Related_to_leukemia_clone.x=first.match.row[1,"Related_to_leukemia_clone"], 212 paste = first.rows[1,"paste"],
216 Frequency.x=first.match.row[1,"Frequency"], 213 min_cell_paste = first.rows[1,"min_cell_paste"])
217 locus_V.x=first.match.row[1,"locus_V"], 214
218 locus_J.x=first.match.row[1,"locus_J"], 215 if(nrow(second.rows) > 0){
219 min_cell_count.x=first.match.row[1,"min_cell_count"], 216 second.sum = data.frame(merge = first.clone.sequence,
220 normalized_read_count.x=first.match.row[1,"normalized_read_count"], 217 Patient = patient,
221 paste.x=first.match.row[1,"paste"], 218 Receptor = second.rows[1,"Receptor"],
222 min_cell_paste.y=second.match.row[1,"min_cell_paste"], 219 Sample = second.rows[1,"Sample"],
223 Patient.y=second.match.row[1,"Patient"], 220 Cell_Count = second.rows[1,"Cell_Count"],
224 Receptor.y=second.match.row[1,"Receptor"], 221 Clone_Molecule_Count_From_Spikes = sum(second.rows$Clone_Molecule_Count_From_Spikes),
225 Sample.y=second.match.row[1,"Sample"], 222 Log10_Frequency = log10(sum(second.rows$Frequency)),
226 Cell_Count.y=second.match.row[1,"Cell_Count"], 223 Total_Read_Count = sum(second.rows$Total_Read_Count),
227 Clone_Molecule_Count_From_Spikes.y=second.match.row[1,"Clone_Molecule_Count_From_Spikes"], 224 dsPerM = sum(second.rows$dsPerM),
228 Log10_Frequency.y=second.match.row[1,"Log10_Frequency"], 225 J_Segment_Major_Gene = sort(table(second.rows$J_Segment_Major_Gene),decreasing=TRUE)[1],
229 Total_Read_Count.y=second.match.row[1,"Total_Read_Count"], 226 V_Segment_Major_Gene = sort(table(second.rows$V_Segment_Major_Gene),decreasing=TRUE)[1],
230 dsPerM.y=second.match.row[1,"dsPerM"], 227 Clone_Sequence = first.clone.sequence,
231 J_Segment_Major_Gene.y=second.match.row[1,"J_Segment_Major_Gene"], 228 CDR3_Sense_Sequence = second.rows[1,"CDR3_Sense_Sequence"],
232 V_Segment_Major_Gene.y=second.match.row[1,"V_Segment_Major_Gene"], 229 Related_to_leukemia_clone = F,
233 Clone_Sequence.y=second.match.row[1,"Clone_Sequence"], 230 Frequency = sum(second.rows$Frequency),
234 CDR3_Sense_Sequence.y=second.match.row[1,"CDR3_Sense_Sequence"], 231 locus_V = second.rows[1,"locus_V"],
235 Related_to_leukemia_clone.y=second.match.row[1,"Related_to_leukemia_clone"], 232 locus_J = second.rows[1,"locus_J"],
236 Frequency.y=second.match.row[1,"Frequency"], 233 min_cell_count = second.rows[1,"min_cell_count"],
237 locus_V.y=second.match.row[1,"locus_V"], 234 normalized_read_count = sum(second.rows$normalized_read_count),
238 locus_J.y=second.match.row[1,"locus_J"], 235 paste = second.rows[1,"paste"],
239 min_cell_count.y=second.match.row[1,"min_cell_count"], 236 min_cell_paste = second.rows[1,"min_cell_paste"])
240 normalized_read_count.y=second.match.row[1,"normalized_read_count"], 237
241 paste.y=first.match.row[1,"paste"]) 238 patientMerge = rbind(patientMerge, merge(first.sum, second.sum, by="merge"))
242 239 patient.fuzzy = patient.fuzzy[!(first.match.filter | second.match.filter),]
243 240
244 patientMerge = rbind(patientMerge, patientMerge.new.row) 241
245 patient.fuzzy = patient.fuzzy[-match.filter,] 242 if(sum(first.match.filter) == 1 & sum(second.match.filter) == 1){
246 243 second.clone.sequence = patient.fuzzy[second.match.filter, "Clone_Sequence"]
247 patient1 = patient1[!(patient1$Clone_Sequence %in% c(first.clone.sequence, second.clone.sequence)),] 244 if(nchar(first.clone.sequence) == nchar(second.clone.sequence)){
248 patient2 = patient2[!(patient2$Clone_Sequence %in% c(first.clone.sequence, second.clone.sequence)),] 245 merge.list[["second"]] = append(merge.list[["second"]], second.clone.sequence)
249 246 }
250 scatterplot_data = scatterplot_data[scatterplot_data$merge != second.clone.sequence,] 247 }
251 248
252 } else if (sum(match.filter) > 1){ 249 if(nrow(first.rows) > 1 | nrow(second.rows) > 1){
253 cat(paste("<tr><td>", "Multiple matches (", sum(match.filter), ") found for", first.merge, "in", patient, "</td></tr>", sep=" "), file=logfile, append=T) 250
254 patient.fuzzy = patient.fuzzy[-1,] 251 }
252
255 } else { 253 } else {
256 patient.fuzzy = patient.fuzzy[-1,] 254 patient.fuzzy = patient.fuzzy[-1,]
257 } 255 }
258
259
260 } 256 }
261 257
262 } 258 }
263 259
264 260
301 filenameTwo = paste(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="") 297 filenameTwo = paste(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
302 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) 298 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
303 } 299 }
304 } else { 300 } else {
305 scatterplot_locus_data = scatterplot_data[grepl(V_Segment, scatterplot_data$V_Segment_Major_Gene) & grepl(J_Segment, scatterplot_data$J_Segment_Major_Gene),] 301 scatterplot_locus_data = scatterplot_data[grepl(V_Segment, scatterplot_data$V_Segment_Major_Gene) & grepl(J_Segment, scatterplot_data$J_Segment_Major_Gene),]
302 #scatterplot_locus_data = scatterplot_locus_data[!(scatterplot_locus_data$merge %in% merge.list[[twoSample]]),]
303 scatterplot_locus_data = scatterplot_locus_data[!(scatterplot_locus_data$merge %in% merge.list[["second"]]),]
306 if(nrow(scatterplot_locus_data) > 0){ 304 if(nrow(scatterplot_locus_data) > 0){
307 scatterplot_locus_data$Rearrangement = product[iter, titleIndex] 305 scatterplot_locus_data$Rearrangement = product[iter, titleIndex]
308 } 306 }
309 in_one = (scatterplot_locus_data$merge %in% patient1$merge) 307 in_one = (scatterplot_locus_data$merge %in% patient1$merge)
310 in_two = (scatterplot_locus_data$merge %in% patient2$merge) 308 in_two = (scatterplot_locus_data$merge %in% patient2$merge)
311 not_in_one = !in_one
312 if(any(in_two)){ 309 if(any(in_two)){
313 scatterplot_locus_data[not_in_one,]$type = twoSample 310 scatterplot_locus_data[in_two,]$type = twoSample
314 } 311 }
315 in_both = (scatterplot_locus_data$merge %in% patientMerge[both,]$merge) 312 in_both = (scatterplot_locus_data$merge %in% patientMerge$merge)
313 #merge.list.filter = (scatterplot_locus_data$merge %in% merge.list[[oneSample]])
314 #exact.matches.filter = (scatterplot_locus_data$merge %in% cs.exact.matches)
316 if(any(in_both)){ 315 if(any(in_both)){
317 scatterplot_locus_data[in_both,]$type = "In Both" 316 scatterplot_locus_data[in_both,]$type = "In Both"
318 } 317 }
319 if(type == "single"){ 318 if(type == "single"){
320 single_patients <<- rbind(single_patients, scatterplot_locus_data) 319 single_patients <<- rbind(single_patients, scatterplot_locus_data)
321 } 320 }
322 p = NULL 321 p = NULL
323 if(nrow(scatterplot_locus_data) != 0){ 322 if(nrow(scatterplot_locus_data) != 0){
324 if(on == "normalized_read_count"){ 323 if(on == "normalized_read_count"){
325 scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count)))) 324 scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count))))
326 p = ggplot(scatterplot_locus_data, aes(type, normalized_read_count)) + scale_y_log10(breaks=scales,labels=scales) + expand_limits(y=10^6) 325 p = ggplot(scatterplot_locus_data, aes(type, normalized_read_count)) + scale_y_log10(breaks=scales,labels=scales) + expand_limits(y=10^6) + scale_x_discrete(breaks=levels(scatterplot_data$type), labels=levels(scatterplot_data$type), drop=FALSE)
327 } else { 326 } else {
328 p = ggplot(scatterplot_locus_data, aes(type, Frequency)) + scale_y_continuous(limits = c(0, 100)) + expand_limits(y=c(0,100)) 327 p = ggplot(scatterplot_locus_data, aes(type, Frequency)) + scale_y_continuous(limits = c(0, 100)) + expand_limits(y=c(0,100)) + scale_x_discrete(breaks=levels(scatterplot_data$type), labels=levels(scatterplot_data$type), drop=FALSE)
329 } 328 }
330 p = p + geom_point(aes(colour=type), position="jitter") 329 p = p + geom_point(aes(colour=type), position="jitter")
331 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])) 330 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]))
332 } else { 331 } else {
333 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])) 332 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]))