Mercurial > repos > davidvanzessen > prisca
changeset 1:75853bceec00 draft
Uploaded
| author | davidvanzessen |
|---|---|
| date | Tue, 17 Jan 2017 07:24:44 -0500 |
| parents | ed6885c85660 |
| children | 7ffd0fba8cf4 |
| files | ALL.xml RScript.r wrapper.sh |
| diffstat | 3 files changed, 494 insertions(+), 500 deletions(-) [+] |
line wrap: on
line diff
--- a/ALL.xml Wed Aug 31 05:31:47 2016 -0400 +++ b/ALL.xml Tue Jan 17 07:24:44 2017 -0500 @@ -8,8 +8,8 @@ <param name="min_freq" type="text" label="Minimum Frequency, between 0 and 100 in percentage" value='0'/> <param name="min_cells" type="text" label="Minimum cell count" value='0'/> <param name="merge_on" type="select" label="Merge On"> - <option value="Clone_Sequence">Clone_Sequence</option> - <option value="V_J_CDR3">V+J+CDR3</option> + <option value="Clone_Sequence">Clone_Sequence</option> + <option value="V_J_CDR3">V+J+CDR3</option> </param> </inputs> <outputs> @@ -40,6 +40,7 @@ | CDR3_Sense_Sequence | The CDR3 sequence region. | +----------------------------------+----------------------------------------------+ +It will handle patients with one, two or three samples. And generate a detailed HTML report on the sequences found in indiviual samples and in both samples. </help>
--- a/RScript.r Wed Aug 31 05:31:47 2016 -0400 +++ b/RScript.r Tue Jan 17 07:24:44 2017 -0500 @@ -34,7 +34,19 @@ dat = dat[dat$Frequency >= min_freq,] -triplets = dat[grepl("VanDongen_cALL_14696", dat$Patient) | grepl("(16278)|(26402)|(26759)", dat$Sample),] +patient.sample.counts = data.frame(data.table(dat)[, list(count=.N), by=c("Patient", "Sample")]) +patient.sample.counts = data.frame(data.table(patient.sample.counts)[, list(count=.N), by=c("Patient")]) + +print("Found the following patients with number of samples:") +print(patient.sample.counts) + +patient.sample.counts.pairs = patient.sample.counts[patient.sample.counts$count %in% 1:2,] +patient.sample.counts.triplets = patient.sample.counts[patient.sample.counts$count == 3,] + + + +triplets = dat[dat$Patient %in% patient.sample.counts.triplets$Patient,] +dat = dat[dat$Patient %in% patient.sample.counts.pairs$Patient,] cat("<tr><td>Normalizing to lowest cell count within locus</td></tr>", file=logfile, append=T) @@ -475,21 +487,21 @@ print(plt) dev.off() } - -cat("<tr><td>Starting Frequency analysis</td></tr>", file=logfile, append=T) +if(length(patients) > 0){ + cat("<tr><td>Starting Frequency analysis</td></tr>", file=logfile, append=T) -interval = intervalFreq -intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval)) -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))) -lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T) + interval = intervalFreq + intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval)) + 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))) + lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T) -cat("<tr><td>Starting Cell Count analysis</td></tr>", file=logfile, append=T) + cat("<tr><td>Starting Cell Count analysis</td></tr>", file=logfile, append=T) -interval = intervalReads -intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval)) -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))) -lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="normalized_read_count") - + interval = intervalReads + intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval)) + 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))) + lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="normalized_read_count") +} if(nrow(single_patients) > 0){ scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count)))) p = ggplot(single_patients, aes(Rearrangement, normalized_read_count)) + scale_y_log10(breaks=scales,labels=as.character(scales)) + expand_limits(y=c(0,1000000)) @@ -525,551 +537,532 @@ patient.merge.list.second = list() tripletAnalysis <- function(patient1, label1, patient2, label2, patient3, label3, product, interval, on, appendTriplets= FALSE){ - onShort = "reads" - if(on == "Frequency"){ - onShort = "freq" - } - onx = paste(on, ".x", sep="") - ony = paste(on, ".y", sep="") - onz = paste(on, ".z", sep="") - type="triplet" - - threshholdIndex = which(colnames(product) == "interval") - V_SegmentIndex = which(colnames(product) == "V_Segments") - J_SegmentIndex = which(colnames(product) == "J_Segments") - titleIndex = which(colnames(product) == "Titles") - sampleIndex = which(colnames(patient1) == "Sample") - patientIndex = which(colnames(patient1) == "Patient") - oneSample = paste(patient1[1,sampleIndex], sep="") - twoSample = paste(patient2[1,sampleIndex], sep="") - threeSample = paste(patient3[1,sampleIndex], sep="") + onShort = "reads" + if(on == "Frequency"){ + onShort = "freq" + } + onx = paste(on, ".x", sep="") + ony = paste(on, ".y", sep="") + onz = paste(on, ".z", sep="") + type="triplet" - if(mergeOn == "Clone_Sequence"){ - patient1$merge = paste(patient1$Clone_Sequence) + threshholdIndex = which(colnames(product) == "interval") + V_SegmentIndex = which(colnames(product) == "V_Segments") + J_SegmentIndex = which(colnames(product) == "J_Segments") + titleIndex = which(colnames(product) == "Titles") + sampleIndex = which(colnames(patient1) == "Sample") + patientIndex = which(colnames(patient1) == "Patient") + oneSample = paste(patient1[1,sampleIndex], sep="") + twoSample = paste(patient2[1,sampleIndex], sep="") + threeSample = paste(patient3[1,sampleIndex], sep="") + + if(mergeOn == "Clone_Sequence"){ + patient1$merge = paste(patient1$Clone_Sequence) patient2$merge = paste(patient2$Clone_Sequence) patient3$merge = paste(patient3$Clone_Sequence) - } else { + } else { patient1$merge = paste(patient1$V_Segment_Major_Gene, patient1$J_Segment_Major_Gene, patient1$CDR3_Sense_Sequence) patient2$merge = paste(patient2$V_Segment_Major_Gene, patient2$J_Segment_Major_Gene, patient2$CDR3_Sense_Sequence) patient3$merge = paste(patient3$V_Segment_Major_Gene, patient3$J_Segment_Major_Gene, patient3$CDR3_Sense_Sequence) - } + } - #patientMerge = merge(patient1, patient2, by="merge")[NULL,] - patient1.fuzzy = patient1 - patient2.fuzzy = patient2 - patient3.fuzzy = patient3 + #patientMerge = merge(patient1, patient2, by="merge")[NULL,] + patient1.fuzzy = patient1 + patient2.fuzzy = patient2 + patient3.fuzzy = patient3 - cat(paste("<tr><td>", label1, "</td>", sep=""), file=logfile, append=T) + cat(paste("<tr><td>", label1, "</td>", sep=""), file=logfile, append=T) - patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J) - patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J) - patient3.fuzzy$merge = paste(patient3.fuzzy$locus_V, patient3.fuzzy$locus_J) + patient1.fuzzy$merge = paste(patient1.fuzzy$locus_V, patient1.fuzzy$locus_J) + patient2.fuzzy$merge = paste(patient2.fuzzy$locus_V, patient2.fuzzy$locus_J) + patient3.fuzzy$merge = paste(patient3.fuzzy$locus_V, patient3.fuzzy$locus_J) - patient.fuzzy = rbind(patient1.fuzzy ,patient2.fuzzy, patient3.fuzzy) - patient.fuzzy = patient.fuzzy[order(nchar(patient.fuzzy$Clone_Sequence)),] + patient.fuzzy = rbind(patient1.fuzzy ,patient2.fuzzy, patient3.fuzzy) + patient.fuzzy = patient.fuzzy[order(nchar(patient.fuzzy$Clone_Sequence)),] - other.sample.list = list() - other.sample.list[[oneSample]] = c(twoSample, threeSample) - other.sample.list[[twoSample]] = c(oneSample, threeSample) - other.sample.list[[threeSample]] = c(oneSample, twoSample) + other.sample.list = list() + other.sample.list[[oneSample]] = c(twoSample, threeSample) + other.sample.list[[twoSample]] = c(oneSample, threeSample) + other.sample.list[[threeSample]] = c(oneSample, twoSample) - patientMerge = merge(patient1, patient2, by="merge") - patientMerge = merge(patientMerge, patient3, by="merge") - colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge)))] = paste(colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge), perl=T))], ".z", sep="") - #patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz]) - patientMerge = patientMerge[NULL,] + patientMerge = merge(patient1, patient2, by="merge") + patientMerge = merge(patientMerge, patient3, by="merge") + colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge)))] = paste(colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge), perl=T))], ".z", sep="") + #patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz]) + patientMerge = patientMerge[NULL,] - duo.merge.list = list() + duo.merge.list = list() - patientMerge12 = merge(patient1, patient2, by="merge") - #patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony]) - patientMerge12 = patientMerge12[NULL,] - duo.merge.list[[paste(oneSample, twoSample)]] = patientMerge12 - duo.merge.list[[paste(twoSample, oneSample)]] = patientMerge12 + patientMerge12 = merge(patient1, patient2, by="merge") + #patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony]) + patientMerge12 = patientMerge12[NULL,] + duo.merge.list[[paste(oneSample, twoSample)]] = patientMerge12 + duo.merge.list[[paste(twoSample, oneSample)]] = patientMerge12 - patientMerge13 = merge(patient1, patient3, by="merge") - #patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony]) - patientMerge13 = patientMerge13[NULL,] - duo.merge.list[[paste(oneSample, threeSample)]] = patientMerge13 - duo.merge.list[[paste(threeSample, oneSample)]] = patientMerge13 + patientMerge13 = merge(patient1, patient3, by="merge") + #patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony]) + patientMerge13 = patientMerge13[NULL,] + duo.merge.list[[paste(oneSample, threeSample)]] = patientMerge13 + duo.merge.list[[paste(threeSample, oneSample)]] = patientMerge13 - patientMerge23 = merge(patient2, patient3, by="merge") - #patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony]) - patientMerge23 = patientMerge23[NULL,] - duo.merge.list[[paste(twoSample, threeSample)]] = patientMerge23 - duo.merge.list[[paste(threeSample, twoSample)]] = patientMerge23 - - merge.list = list() - merge.list[["second"]] = vector() + patientMerge23 = merge(patient2, patient3, by="merge") + #patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony]) + patientMerge23 = patientMerge23[NULL,] + duo.merge.list[[paste(twoSample, threeSample)]] = patientMerge23 + duo.merge.list[[paste(threeSample, twoSample)]] = patientMerge23 - start.time = proc.time() - if(paste(label1, "123") %in% names(patient.merge.list)){ - patientMerge = patient.merge.list[[paste(label1, "123")]] - patientMerge12 = patient.merge.list[[paste(label1, "12")]] - patientMerge13 = patient.merge.list[[paste(label1, "13")]] - patientMerge23 = patient.merge.list[[paste(label1, "23")]] + merge.list = list() + merge.list[["second"]] = vector() + + #print(paste(nrow(patient1), nrow(patient2), nrow(patient3), label1, label2, label3)) + + start.time = proc.time() + if(paste(label1, "123") %in% names(patient.merge.list)){ + patientMerge = patient.merge.list[[paste(label1, "123")]] + patientMerge12 = patient.merge.list[[paste(label1, "12")]] + patientMerge13 = patient.merge.list[[paste(label1, "13")]] + patientMerge23 = patient.merge.list[[paste(label1, "23")]] - merge.list[["second"]] = patient.merge.list.second[[label1]] + #merge.list[["second"]] = patient.merge.list.second[[label1]] - 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) - } else { - while(nrow(patient.fuzzy) > 0){ - first.merge = patient.fuzzy[1,"merge"] - first.clone.sequence = patient.fuzzy[1,"Clone_Sequence"] - first.sample = patient.fuzzy[1,"Sample"] - - merge.filter = first.merge == patient.fuzzy$merge - - second.sample = other.sample.list[[first.sample]][1] - third.sample = other.sample.list[[first.sample]][2] + 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) + } else { + while(nrow(patient.fuzzy) > 0){ + first.merge = patient.fuzzy[1,"merge"] + first.clone.sequence = patient.fuzzy[1,"Clone_Sequence"] + first.sample = paste(patient.fuzzy[1,"Sample"], sep="") + + merge.filter = first.merge == patient.fuzzy$merge + + second.sample = other.sample.list[[first.sample]][1] + third.sample = other.sample.list[[first.sample]][2] - sample.filter.1 = first.sample == patient.fuzzy$Sample - sample.filter.2 = second.sample == patient.fuzzy$Sample - sample.filter.3 = third.sample == patient.fuzzy$Sample + sample.filter.1 = first.sample == patient.fuzzy$Sample + sample.filter.2 = second.sample == patient.fuzzy$Sample + sample.filter.3 = third.sample == patient.fuzzy$Sample - sequence.filter = grepl(paste("^", first.clone.sequence, sep=""), patient.fuzzy$Clone_Sequence) + sequence.filter = grepl(paste("^", first.clone.sequence, sep=""), patient.fuzzy$Clone_Sequence) - match.filter.1 = sample.filter.1 & sequence.filter & merge.filter - match.filter.2 = sample.filter.2 & sequence.filter & merge.filter - match.filter.3 = sample.filter.3 & sequence.filter & merge.filter + match.filter.1 = sample.filter.1 & sequence.filter & merge.filter + match.filter.2 = sample.filter.2 & sequence.filter & merge.filter + match.filter.3 = sample.filter.3 & sequence.filter & merge.filter - matches.in.1 = any(match.filter.1) - matches.in.2 = any(match.filter.2) - matches.in.3 = any(match.filter.3) + matches.in.1 = any(match.filter.1) + matches.in.2 = any(match.filter.2) + matches.in.3 = any(match.filter.3) - - - rows.1 = patient.fuzzy[match.filter.1,] + rows.1 = patient.fuzzy[match.filter.1,] - sum.1 = data.frame(merge = first.clone.sequence, - Patient = label1, - Receptor = rows.1[1,"Receptor"], - Sample = rows.1[1,"Sample"], - Cell_Count = rows.1[1,"Cell_Count"], - Clone_Molecule_Count_From_Spikes = sum(rows.1$Clone_Molecule_Count_From_Spikes), - Log10_Frequency = log10(sum(rows.1$Frequency)), - Total_Read_Count = sum(rows.1$Total_Read_Count), - dsPerM = sum(rows.1$dsPerM), - J_Segment_Major_Gene = rows.1[1,"J_Segment_Major_Gene"], - V_Segment_Major_Gene = rows.1[1,"V_Segment_Major_Gene"], - Clone_Sequence = first.clone.sequence, - CDR3_Sense_Sequence = rows.1[1,"CDR3_Sense_Sequence"], - Related_to_leukemia_clone = F, - Frequency = sum(rows.1$Frequency), - locus_V = rows.1[1,"locus_V"], - locus_J = rows.1[1,"locus_J"], - uniqueID = rows.1[1,"uniqueID"], - normalized_read_count = sum(rows.1$normalized_read_count)) - sum.2 = sum.1[NULL,] - rows.2 = patient.fuzzy[match.filter.2,] - if(matches.in.2){ - sum.2 = data.frame(merge = first.clone.sequence, - Patient = label1, - Receptor = rows.2[1,"Receptor"], - Sample = rows.2[1,"Sample"], - Cell_Count = rows.2[1,"Cell_Count"], - Clone_Molecule_Count_From_Spikes = sum(rows.2$Clone_Molecule_Count_From_Spikes), - Log10_Frequency = log10(sum(rows.2$Frequency)), - Total_Read_Count = sum(rows.2$Total_Read_Count), - dsPerM = sum(rows.2$dsPerM), - J_Segment_Major_Gene = rows.2[1,"J_Segment_Major_Gene"], - V_Segment_Major_Gene = rows.2[1,"V_Segment_Major_Gene"], - Clone_Sequence = first.clone.sequence, - CDR3_Sense_Sequence = rows.2[1,"CDR3_Sense_Sequence"], - Related_to_leukemia_clone = F, - Frequency = sum(rows.2$Frequency), - locus_V = rows.2[1,"locus_V"], - locus_J = rows.2[1,"locus_J"], - uniqueID = rows.2[1,"uniqueID"], - normalized_read_count = sum(rows.2$normalized_read_count)) - } + sum.1 = data.frame(merge = first.clone.sequence, + Patient = label1, + Receptor = rows.1[1,"Receptor"], + Sample = rows.1[1,"Sample"], + Cell_Count = rows.1[1,"Cell_Count"], + Clone_Molecule_Count_From_Spikes = sum(rows.1$Clone_Molecule_Count_From_Spikes), + Log10_Frequency = log10(sum(rows.1$Frequency)), + Total_Read_Count = sum(rows.1$Total_Read_Count), + dsPerM = sum(rows.1$dsPerM), + J_Segment_Major_Gene = rows.1[1,"J_Segment_Major_Gene"], + V_Segment_Major_Gene = rows.1[1,"V_Segment_Major_Gene"], + Clone_Sequence = first.clone.sequence, + CDR3_Sense_Sequence = rows.1[1,"CDR3_Sense_Sequence"], + Related_to_leukemia_clone = F, + Frequency = sum(rows.1$Frequency), + locus_V = rows.1[1,"locus_V"], + locus_J = rows.1[1,"locus_J"], + uniqueID = rows.1[1,"uniqueID"], + normalized_read_count = sum(rows.1$normalized_read_count)) + sum.2 = sum.1[NULL,] + rows.2 = patient.fuzzy[match.filter.2,] + if(matches.in.2){ + sum.2 = data.frame(merge = first.clone.sequence, + Patient = label1, + Receptor = rows.2[1,"Receptor"], + Sample = rows.2[1,"Sample"], + Cell_Count = rows.2[1,"Cell_Count"], + Clone_Molecule_Count_From_Spikes = sum(rows.2$Clone_Molecule_Count_From_Spikes), + Log10_Frequency = log10(sum(rows.2$Frequency)), + Total_Read_Count = sum(rows.2$Total_Read_Count), + dsPerM = sum(rows.2$dsPerM), + J_Segment_Major_Gene = rows.2[1,"J_Segment_Major_Gene"], + V_Segment_Major_Gene = rows.2[1,"V_Segment_Major_Gene"], + Clone_Sequence = first.clone.sequence, + CDR3_Sense_Sequence = rows.2[1,"CDR3_Sense_Sequence"], + Related_to_leukemia_clone = F, + Frequency = sum(rows.2$Frequency), + locus_V = rows.2[1,"locus_V"], + locus_J = rows.2[1,"locus_J"], + uniqueID = rows.2[1,"uniqueID"], + normalized_read_count = sum(rows.2$normalized_read_count)) + } - sum.3 = sum.1[NULL,] - rows.3 = patient.fuzzy[match.filter.3,] - if(matches.in.3){ - sum.3 = data.frame(merge = first.clone.sequence, - Patient = label1, - Receptor = rows.3[1,"Receptor"], - Sample = rows.3[1,"Sample"], - Cell_Count = rows.3[1,"Cell_Count"], - Clone_Molecule_Count_From_Spikes = sum(rows.3$Clone_Molecule_Count_From_Spikes), - Log10_Frequency = log10(sum(rows.3$Frequency)), - Total_Read_Count = sum(rows.3$Total_Read_Count), - dsPerM = sum(rows.3$dsPerM), - J_Segment_Major_Gene = rows.3[1,"J_Segment_Major_Gene"], - V_Segment_Major_Gene = rows.3[1,"V_Segment_Major_Gene"], - Clone_Sequence = first.clone.sequence, - CDR3_Sense_Sequence = rows.3[1,"CDR3_Sense_Sequence"], - Related_to_leukemia_clone = F, - Frequency = sum(rows.3$Frequency), - locus_V = rows.3[1,"locus_V"], - locus_J = rows.3[1,"locus_J"], - uniqueID = rows.3[1,"uniqueID"], - normalized_read_count = sum(rows.3$normalized_read_count)) - } + sum.3 = sum.1[NULL,] + rows.3 = patient.fuzzy[match.filter.3,] + if(matches.in.3){ + sum.3 = data.frame(merge = first.clone.sequence, + Patient = label1, + Receptor = rows.3[1,"Receptor"], + Sample = rows.3[1,"Sample"], + Cell_Count = rows.3[1,"Cell_Count"], + Clone_Molecule_Count_From_Spikes = sum(rows.3$Clone_Molecule_Count_From_Spikes), + Log10_Frequency = log10(sum(rows.3$Frequency)), + Total_Read_Count = sum(rows.3$Total_Read_Count), + dsPerM = sum(rows.3$dsPerM), + J_Segment_Major_Gene = rows.3[1,"J_Segment_Major_Gene"], + V_Segment_Major_Gene = rows.3[1,"V_Segment_Major_Gene"], + Clone_Sequence = first.clone.sequence, + CDR3_Sense_Sequence = rows.3[1,"CDR3_Sense_Sequence"], + Related_to_leukemia_clone = F, + Frequency = sum(rows.3$Frequency), + locus_V = rows.3[1,"locus_V"], + locus_J = rows.3[1,"locus_J"], + uniqueID = rows.3[1,"uniqueID"], + normalized_read_count = sum(rows.3$normalized_read_count)) + } - if(matches.in.2 & matches.in.3){ - merge.123 = merge(sum.1, sum.2, by="merge") - merge.123 = merge(merge.123, sum.3, by="merge") - 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="") - #merge.123$thresholdValue = pmax(merge.123[,onx], merge.123[,ony], merge.123[,onz]) + if(matches.in.2 & matches.in.3){ + merge.123 = merge(sum.1, sum.2, by="merge") + merge.123 = merge(merge.123, sum.3, by="merge") + 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="") + #merge.123$thresholdValue = pmax(merge.123[,onx], merge.123[,ony], merge.123[,onz]) - patientMerge = rbind(patientMerge, merge.123) - patient.fuzzy = patient.fuzzy[!(match.filter.1 | match.filter.2 | match.filter.3),] + patientMerge = rbind(patientMerge, merge.123) + patient.fuzzy = patient.fuzzy[!(match.filter.1 | match.filter.2 | match.filter.3),] - 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"]) - merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences) + 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"]) + merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences) - } else if (matches.in.2) { - #other.sample1 = other.sample.list[[first.sample]][1] - #other.sample2 = other.sample.list[[first.sample]][2] + } else if (matches.in.2) { + #other.sample1 = other.sample.list[[first.sample]][1] + #other.sample2 = other.sample.list[[first.sample]][2] - second.sample = sum.2[,"Sample"] + second.sample = sum.2[,"Sample"] - current.merge.list = duo.merge.list[[paste(first.sample, second.sample)]] + current.merge.list = duo.merge.list[[paste(first.sample, second.sample)]] - merge.12 = merge(sum.1, sum.2, by="merge") + merge.12 = merge(sum.1, sum.2, by="merge") - current.merge.list = rbind(current.merge.list, merge.12) - duo.merge.list[[paste(first.sample, second.sample)]] = current.merge.list + current.merge.list = rbind(current.merge.list, merge.12) + duo.merge.list[[paste(first.sample, second.sample)]] = current.merge.list - patient.fuzzy = patient.fuzzy[!(match.filter.1 | match.filter.2),] + patient.fuzzy = patient.fuzzy[!(match.filter.1 | match.filter.2),] - hidden.clone.sequences = c(rows.1[-1,"Clone_Sequence"], rows.2[rows.2$Clone_Sequence != first.clone.sequence,"Clone_Sequence"]) - merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences) + hidden.clone.sequences = c(rows.1[-1,"Clone_Sequence"], rows.2[rows.2$Clone_Sequence != first.clone.sequence,"Clone_Sequence"]) + merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences) - } else if (matches.in.3) { + } else if (matches.in.3) { - #other.sample1 = other.sample.list[[first.sample]][1] - #other.sample2 = other.sample.list[[first.sample]][2] + #other.sample1 = other.sample.list[[first.sample]][1] + #other.sample2 = other.sample.list[[first.sample]][2] - second.sample = sum.3[,"Sample"] + second.sample = sum.3[,"Sample"] - current.merge.list = duo.merge.list[[paste(first.sample, second.sample)]] + current.merge.list = duo.merge.list[[paste(first.sample, second.sample)]] - merge.13 = merge(sum.1, sum.3, by="merge") + merge.13 = merge(sum.1, sum.3, by="merge") - current.merge.list = rbind(current.merge.list, merge.13) - duo.merge.list[[paste(first.sample, second.sample)]] = current.merge.list + current.merge.list = rbind(current.merge.list, merge.13) + duo.merge.list[[paste(first.sample, second.sample)]] = current.merge.list - patient.fuzzy = patient.fuzzy[!(match.filter.1 | match.filter.3),] + patient.fuzzy = patient.fuzzy[!(match.filter.1 | match.filter.3),] - hidden.clone.sequences = c(rows.1[-1,"Clone_Sequence"], rows.3[rows.3$Clone_Sequence != first.clone.sequence,"Clone_Sequence"]) - merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences) + hidden.clone.sequences = c(rows.1[-1,"Clone_Sequence"], rows.3[rows.3$Clone_Sequence != first.clone.sequence,"Clone_Sequence"]) + merge.list[["second"]] = append(merge.list[["second"]], hidden.clone.sequences) - } else if(nrow(rows.1) > 1){ - patient1 = patient1[!(patient1$Clone_Sequence %in% rows.1$Clone_Sequence),] - print(names(patient1)[names(patient1) %in% sum.1]) - print(names(patient1)[!(names(patient1) %in% sum.1)]) - print(names(patient1)) - print(names(sum.1)) - print(summary(sum.1)) - print(summary(patient1)) - print(dim(sum.1)) - print(dim(patient1)) - print(head(sum.1[,names(patient1)])) - patient1 = rbind(patient1, sum.1[,names(patient1)]) - patient.fuzzy = patient.fuzzy[-match.filter.1,] - } else { - patient.fuzzy = patient.fuzzy[-1,] - } + } else if(nrow(rows.1) > 1){ + patient1 = patient1[!(patient1$Clone_Sequence %in% rows.1$Clone_Sequence),] + print(names(patient1)[names(patient1) %in% sum.1]) + print(names(patient1)[!(names(patient1) %in% sum.1)]) + print(names(patient1)) + print(names(sum.1)) + print(summary(sum.1)) + print(summary(patient1)) + print(dim(sum.1)) + print(dim(patient1)) + print(head(sum.1[,names(patient1)])) + patient1 = rbind(patient1, sum.1[,names(patient1)]) + patient.fuzzy = patient.fuzzy[-match.filter.1,] + } else { + patient.fuzzy = patient.fuzzy[-1,] + } + + tmp.rows = rbind(rows.1, rows.2, rows.3) + tmp.rows = tmp.rows[order(nchar(tmp.rows$Clone_Sequence)),] - tmp.rows = rbind(rows.1, rows.2, rows.3) - tmp.rows = tmp.rows[order(nchar(tmp.rows$Clone_Sequence)),] + if (sum(match.filter.1) > 1 | sum(match.filter.2) > 1 | sum(match.filter.1) > 1) { + 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) + } else { + } + + } + patient.merge.list[[paste(label1, "123")]] = patientMerge + + patientMerge12 = duo.merge.list[[paste(oneSample, twoSample)]] + patientMerge13 = duo.merge.list[[paste(oneSample, threeSample)]] + patientMerge23 = duo.merge.list[[paste(twoSample, threeSample)]] - if (sum(match.filter.1) > 1 | sum(match.filter.2) > 1 | sum(match.filter.1) > 1) { - 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) - } else { - } + patient.merge.list[[paste(label1, "12")]] = patientMerge12 + patient.merge.list[[paste(label1, "13")]] = patientMerge13 + patient.merge.list[[paste(label1, "23")]] = patientMerge23 + + #patient.merge.list.second[[label1]] = merge.list[["second"]] + } + 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) + patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz]) + patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony]) + patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony]) + patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony]) + + #patientMerge$thresholdValue = pmin(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz]) + #patientMerge12$thresholdValue = pmin(patientMerge12[,onx], patientMerge12[,ony]) + #patientMerge13$thresholdValue = pmin(patientMerge13[,onx], patientMerge13[,ony]) + #patientMerge23$thresholdValue = pmin(patientMerge23[,onx], patientMerge23[,ony]) - } - patient.merge.list[[paste(label1, "123")]] = patientMerge + patient1 = patient1[!(patient1$Clone_Sequence %in% merge.list[["second"]]),] + patient2 = patient2[!(patient2$Clone_Sequence %in% merge.list[["second"]]),] + patient3 = patient3[!(patient3$Clone_Sequence %in% merge.list[["second"]]),] - patientMerge12 = duo.merge.list[[paste(oneSample, twoSample)]] - patientMerge13 = duo.merge.list[[paste(oneSample, threeSample)]] - patientMerge23 = duo.merge.list[[paste(twoSample, threeSample)]] + if(F){ + patientMerge = merge(patient1, patient2, by="merge") + patientMerge = merge(patientMerge, patient3, by="merge") + colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge)))] = paste(colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge), perl=T))], ".z", sep="") + patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz]) + patientMerge12 = merge(patient1, patient2, by="merge") + patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony]) + patientMerge13 = merge(patient1, patient3, by="merge") + patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony]) + patientMerge23 = merge(patient2, patient3, by="merge") + patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony]) + } - patient.merge.list[[paste(label1, "12")]] = patientMerge12 - patient.merge.list[[paste(label1, "13")]] = patientMerge13 - patient.merge.list[[paste(label1, "23")]] = patientMerge23 + scatterplot_data_columns = c("Clone_Sequence", "Frequency", "normalized_read_count", "V_Segment_Major_Gene", "J_Segment_Major_Gene", "merge") + scatterplot_data = rbind(patient1[,scatterplot_data_columns], patient2[,scatterplot_data_columns], patient3[,scatterplot_data_columns]) + scatterplot_data = scatterplot_data[!duplicated(scatterplot_data$merge),] + + scatterplot_data$type = factor(x="In one", levels=c("In one", "In two", "In three", "In multiple")) - patient.merge.list.second[[label1]] = merge.list[["second"]] - } - 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) - patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz]) - patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony]) - patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony]) - patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony]) + res1 = vector() + res2 = vector() + res3 = vector() + res12 = vector() + res13 = vector() + res23 = vector() + resAll = vector() + read1Count = vector() + read2Count = vector() + read3Count = vector() - #patientMerge$thresholdValue = pmin(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz]) - #patientMerge12$thresholdValue = pmin(patientMerge12[,onx], patientMerge12[,ony]) - #patientMerge13$thresholdValue = pmin(patientMerge13[,onx], patientMerge13[,ony]) - #patientMerge23$thresholdValue = pmin(patientMerge23[,onx], patientMerge23[,ony]) + if(appendTriplets){ + cat(paste(label1, label2, label3, sep="\t"), file="triplets.txt", append=T, sep="", fill=3) + } + for(iter in 1:length(product[,1])){ + threshhold = product[iter,threshholdIndex] + V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="") + J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="") + #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) + all = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge$thresholdValue > threshhold) - patient1 = patient1[!(patient1$Clone_Sequence %in% merge.list[["second"]]),] - patient2 = patient2[!(patient2$Clone_Sequence %in% merge.list[["second"]]),] - patient3 = patient3[!(patient3$Clone_Sequence %in% merge.list[["second"]]),] + 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)) + 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)) + 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)) + + 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)) + 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)) + 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)) - if(F){ - patientMerge = merge(patient1, patient2, by="merge") - patientMerge = merge(patientMerge, patient3, by="merge") - colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge)))] = paste(colnames(patientMerge)[which(!grepl("(\\.x$)|(\\.y$)|(merge)", names(patientMerge), perl=T))], ".z", sep="") - patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz]) - patientMerge12 = merge(patient1, patient2, by="merge") - patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony]) - patientMerge13 = merge(patient1, patient3, by="merge") - patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony]) - patientMerge23 = merge(patient2, patient3, by="merge") - patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony]) - } - - scatterplot_data_columns = c("Clone_Sequence", "Frequency", "normalized_read_count", "V_Segment_Major_Gene", "J_Segment_Major_Gene", "merge") - scatterplot_data = rbind(patient1[,scatterplot_data_columns], patient2[,scatterplot_data_columns], patient3[,scatterplot_data_columns]) - scatterplot_data = scatterplot_data[!duplicated(scatterplot_data$merge),] - scatterplot_data$type = factor(x="In one", levels=c("In one", "In two", "In three", "In multiple")) - - res1 = vector() - res2 = vector() - res3 = vector() - res12 = vector() - res13 = vector() - res23 = vector() - resAll = vector() - read1Count = vector() - read2Count = vector() - read3Count = vector() - - if(appendTriplets){ - cat(paste(label1, label2, label3, sep="\t"), file="triplets.txt", append=T, sep="", fill=3) - } - for(iter in 1:length(product[,1])){ - threshhold = product[iter,threshholdIndex] - V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="") - J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="") - #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) - all = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge$thresholdValue > threshhold) - - 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)) - 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)) - 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)) - - 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)) - 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)) - 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)) - - read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.x)) - read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.y)) - read3Count = append(read3Count, sum(patient3[three,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.z)) - res1 = append(res1, sum(one)) - res2 = append(res2, sum(two)) - res3 = append(res3, sum(three)) - resAll = append(resAll, sum(all)) - res12 = append(res12, sum(one_two)) - res13 = append(res13, sum(one_three)) - res23 = append(res23, sum(two_three)) - #threshhold = 0 - if(threshhold != 0){ - if(sum(one) > 0){ - dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")] - colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone") - filenameOne = paste(label1, "_", product[iter, titleIndex], "_", threshhold, sep="") - write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) - } - if(sum(two) > 0){ - dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")] - colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone") - filenameTwo = paste(label2, "_", product[iter, titleIndex], "_", threshhold, sep="") - write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) - } - if(sum(three) > 0){ - dfThree = patient3[three,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")] - colnames(dfThree) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone") - filenameThree = paste(label3, "_", product[iter, titleIndex], "_", threshhold, sep="") - write.table(dfThree, file=paste(filenameThree, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) - } - if(sum(one_two) > 0){ - 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")] - 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)) - filenameOne_two = paste(label1, "_", label2, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="") - write.table(dfOne_two, file=paste(filenameOne_two, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) - } - if(sum(one_three) > 0){ - 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")] - 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)) - filenameOne_three = paste(label1, "_", label3, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="") - write.table(dfOne_three, file=paste(filenameOne_three, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) - } - if(sum(two_three) > 0){ - 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")] - 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)) - filenameTwo_three = paste(label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="") - write.table(dfTwo_three, file=paste(filenameTwo_three, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) - } - } else { #scatterplot data - scatterplot_locus_data = scatterplot_data[grepl(V_Segment, scatterplot_data$V_Segment_Major_Gene) & grepl(J_Segment, scatterplot_data$J_Segment_Major_Gene),] - scatterplot_locus_data = scatterplot_locus_data[!(scatterplot_locus_data$merge %in% merge.list[["second"]]),] - 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) - if(sum(in_two) > 0){ - scatterplot_locus_data[in_two,]$type = "In two" - } - in_three = (scatterplot_locus_data$merge %in% patientMerge[all,]$merge) - if(sum(in_three)> 0){ - scatterplot_locus_data[in_three,]$type = "In three" - } - not_in_one = scatterplot_locus_data$type != "In one" - if(sum(not_in_one) > 0){ - #scatterplot_locus_data[not_in_one,]$type = "In multiple" - } - p = NULL - if(nrow(scatterplot_locus_data) != 0){ - if(on == "normalized_read_count"){ - scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count)))) - p = ggplot(scatterplot_locus_data, aes(type, normalized_read_count)) + scale_y_log10(breaks=scales,labels=scales, limits=c(1, 1e6)) - } else { - 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)) - #p = ggplot(scatterplot_locus_data, aes(type, Frequency)) + scale_y_continuous(limits = c(0, 100)) + expand_limits(y=c(0,100)) - } - p = p + geom_point(aes(colour=type), position="jitter") - p = p + xlab("In one or in multiple samples") + ylab(onShort) + ggtitle(paste(label1, label2, label3, onShort, product[iter, titleIndex])) - } else { - 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])) - } - png(paste(label1, "_", label2, "_", label3, "_", onShort, "_", product[iter, titleIndex],"_scatter.png", sep="")) - print(p) - dev.off() - } - if(sum(all) > 0){ - 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")] - 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)) - filenameAll = paste(label1, "_", label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, sep="") - write.table(dfAll, file=paste(filenameAll, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) - } - } - #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)) - 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) - colnames(patientResult)[6] = oneSample - colnames(patientResult)[7] = twoSample - colnames(patientResult)[8] = threeSample - colnames(patientResult)[9] = paste(oneSample, twoSample, sep="_") - colnames(patientResult)[10] = paste(oneSample, twoSample, sep="_") - colnames(patientResult)[11] = paste(oneSample, twoSample, sep="_") - - colnamesBak = colnames(patientResult) - 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)) - write.table(patientResult, file=paste(label1, "_", label2, "_", label3, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) - colnames(patientResult) = colnamesBak - - patientResult$Locus = factor(patientResult$Locus, Titles) - patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep="")) - - plt = ggplot(patientResult[,c("Locus", "cut_off_value", "All")]) - plt = plt + geom_bar( aes( x=factor(cut_off_value), y=All), stat='identity', position="dodge", fill="#79c36a") - plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) - plt = plt + geom_text(aes(ymax=max(All), x=cut_off_value,y=All,label=All), angle=90, hjust=0) - plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in All") - plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines")) - png(paste(label1, "_", label2, "_", label3, "_", onShort, "_total_all.png", sep=""), width=1920, height=1080) - print(plt) - dev.off() - - fontSize = 4 - - bak = patientResult - patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample, threeSample)] ,id.vars=1:2) - patientResult$relativeValue = patientResult$value * 10 - patientResult[patientResult$relativeValue == 0,]$relativeValue = 1 - plt = ggplot(patientResult) - plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge") - plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) - plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9))) - 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) - 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) - 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) - plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in only one sample") - png(paste(label1, "_", label2, "_", label3, "_", onShort, "_indiv_all.png", sep=""), width=1920, height=1080) - print(plt) - dev.off() + read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.x)) + read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.y)) + read3Count = append(read3Count, sum(patient3[three,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.z)) + res1 = append(res1, sum(one)) + res2 = append(res2, sum(two)) + res3 = append(res3, sum(three)) + resAll = append(resAll, sum(all)) + res12 = append(res12, sum(one_two)) + res13 = append(res13, sum(one_three)) + res23 = append(res23, sum(two_three)) + #threshhold = 0 + if(threshhold != 0){ + if(sum(one) > 0){ + dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")] + colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone") + filenameOne = paste(label1, "_", product[iter, titleIndex], "_", threshhold, sep="") + write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) + } + if(sum(two) > 0){ + dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")] + colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone") + filenameTwo = paste(label2, "_", product[iter, titleIndex], "_", threshhold, sep="") + write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) + } + if(sum(three) > 0){ + dfThree = patient3[three,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")] + colnames(dfThree) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone") + filenameThree = paste(label3, "_", product[iter, titleIndex], "_", threshhold, sep="") + write.table(dfThree, file=paste(filenameThree, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) + } + if(sum(one_two) > 0){ + 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")] + 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)) + filenameOne_two = paste(label1, "_", label2, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="") + write.table(dfOne_two, file=paste(filenameOne_two, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) + } + if(sum(one_three) > 0){ + 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")] + 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)) + filenameOne_three = paste(label1, "_", label3, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="") + write.table(dfOne_three, file=paste(filenameOne_three, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) + } + if(sum(two_three) > 0){ + 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")] + 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)) + filenameTwo_three = paste(label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="") + write.table(dfTwo_three, file=paste(filenameTwo_three, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) + } + } else { #scatterplot data + scatterplot_locus_data = scatterplot_data[grepl(V_Segment, scatterplot_data$V_Segment_Major_Gene) & grepl(J_Segment, scatterplot_data$J_Segment_Major_Gene),] + scatterplot_locus_data = scatterplot_locus_data[!(scatterplot_locus_data$merge %in% merge.list[["second"]]),] + 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) + if(sum(in_two) > 0){ + scatterplot_locus_data[in_two,]$type = "In two" + } + in_three = (scatterplot_locus_data$merge %in% patientMerge[all,]$merge) + if(sum(in_three)> 0){ + scatterplot_locus_data[in_three,]$type = "In three" + } + not_in_one = scatterplot_locus_data$type != "In one" + if(sum(not_in_one) > 0){ + #scatterplot_locus_data[not_in_one,]$type = "In multiple" + } + p = NULL + if(nrow(scatterplot_locus_data) != 0){ + if(on == "normalized_read_count"){ + scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count)))) + p = ggplot(scatterplot_locus_data, aes(type, normalized_read_count)) + scale_y_log10(breaks=scales,labels=scales, limits=c(1, 1e6)) + } else { + 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)) + #p = ggplot(scatterplot_locus_data, aes(type, Frequency)) + scale_y_continuous(limits = c(0, 100)) + expand_limits(y=c(0,100)) + } + p = p + geom_point(aes(colour=type), position="jitter") + p = p + xlab("In one or in multiple samples") + ylab(onShort) + ggtitle(paste(label1, label2, label3, onShort, product[iter, titleIndex])) + } else { + 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])) + } + png(paste(label1, "_", label2, "_", label3, "_", onShort, "_", product[iter, titleIndex],"_scatter.png", sep="")) + print(p) + dev.off() + } + if(sum(all) > 0){ + 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")] + 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)) + filenameAll = paste(label1, "_", label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, sep="") + write.table(dfAll, file=paste(filenameAll, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) + } + } + #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)) + 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) + colnames(patientResult)[6] = oneSample + colnames(patientResult)[7] = twoSample + colnames(patientResult)[8] = threeSample + colnames(patientResult)[9] = paste(oneSample, twoSample, sep="_") + colnames(patientResult)[10] = paste(oneSample, twoSample, sep="_") + colnames(patientResult)[11] = paste(oneSample, twoSample, sep="_") + + colnamesBak = colnames(patientResult) + 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)) + write.table(patientResult, file=paste(label1, "_", label2, "_", label3, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) + colnames(patientResult) = colnamesBak + + patientResult$Locus = factor(patientResult$Locus, Titles) + patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep="")) + + plt = ggplot(patientResult[,c("Locus", "cut_off_value", "All")]) + plt = plt + geom_bar( aes( x=factor(cut_off_value), y=All), stat='identity', position="dodge", fill="#79c36a") + plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + plt = plt + geom_text(aes(ymax=max(All), x=cut_off_value,y=All,label=All), angle=90, hjust=0) + plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in All") + plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines")) + png(paste(label1, "_", label2, "_", label3, "_", onShort, "_total_all.png", sep=""), width=1920, height=1080) + print(plt) + dev.off() + + fontSize = 4 + + bak = patientResult + patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample, threeSample)] ,id.vars=1:2) + patientResult$relativeValue = patientResult$value * 10 + patientResult[patientResult$relativeValue == 0,]$relativeValue = 1 + plt = ggplot(patientResult) + plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge") + plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1)) + plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9))) + 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) + 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) + 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) + plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in only one sample") + png(paste(label1, "_", label2, "_", label3, "_", onShort, "_indiv_all.png", sep=""), width=1920, height=1080) + print(plt) + dev.off() } if(nrow(triplets) != 0){ - cat("<tr><td>Starting triplet analysis</td></tr>", file=logfile, append=T) + cat("<tr><td>Starting triplet analysis</td></tr>", file=logfile, append=T) + + triplets$uniqueID = paste(triplets$Patient, triplets$Sample, sep="_") + + cat("<tr><td>Normalizing to lowest cell count within locus</td></tr>", file=logfile, append=T) - triplets$uniqueID = "ID" - - triplets[grepl("16278_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left" - triplets[grepl("26402_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left" - triplets[grepl("26759_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left" - - triplets[grepl("16278_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right" - triplets[grepl("26402_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right" - triplets[grepl("26759_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right" - - triplets[grepl("14696", triplets$Patient),]$uniqueID = "14696" + triplets$locus_V = substring(triplets$V_Segment_Major_Gene, 0, 4) + triplets$locus_J = substring(triplets$J_Segment_Major_Gene, 0, 4) + min_cell_count = data.frame(data.table(triplets)[, list(min_cell_count=min(.SD$Cell_Count)), by=c("uniqueID", "locus_V", "locus_J")]) - cat("<tr><td>Normalizing to lowest cell count within locus</td></tr>", file=logfile, append=T) + triplets$min_cell_paste = paste(triplets$uniqueID, triplets$locus_V, triplets$locus_J) + min_cell_count$min_cell_paste = paste(min_cell_count$uniqueID, min_cell_count$locus_V, min_cell_count$locus_J) - triplets$locus_V = substring(triplets$V_Segment_Major_Gene, 0, 4) - triplets$locus_J = substring(triplets$J_Segment_Major_Gene, 0, 4) - min_cell_count = data.frame(data.table(triplets)[, list(min_cell_count=min(.SD$Cell_Count)), by=c("uniqueID", "locus_V", "locus_J")]) - - triplets$min_cell_paste = paste(triplets$uniqueID, triplets$locus_V, triplets$locus_J) - min_cell_count$min_cell_paste = paste(min_cell_count$uniqueID, min_cell_count$locus_V, min_cell_count$locus_J) - - min_cell_count = min_cell_count[,c("min_cell_paste", "min_cell_count")] - - triplets = merge(triplets, min_cell_count, by="min_cell_paste") - - 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 - - triplets = triplets[triplets$normalized_read_count >= min_cells,] - - column_drops = c("min_cell_count", "min_cell_paste") - - triplets = triplets[,!(colnames(triplets) %in% column_drops)] + min_cell_count = min_cell_count[,c("min_cell_paste", "min_cell_count")] + + triplets = merge(triplets, min_cell_count, by="min_cell_paste") + + triplets$normalized_read_count = round(triplets$Clone_Molecule_Count_From_Spikes / triplets$Cell_Count * triplets$min_cell_count / 2, digits=2) - cat("<tr><td>Starting Cell Count analysis</td></tr>", file=logfile, append=T) + triplets = triplets[triplets$normalized_read_count >= min_cells,] + + column_drops = c("min_cell_count", "min_cell_paste") + + triplets = triplets[,!(colnames(triplets) %in% column_drops)] + + cat("<tr><td>Starting Cell Count analysis</td></tr>", file=logfile, append=T) - interval = intervalReads - intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval)) - 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))) - - one = triplets[triplets$Sample == "14696_reg_BM",] - two = triplets[triplets$Sample == "24536_reg_BM",] - three = triplets[triplets$Sample == "24062_reg_BM",] - tripletAnalysis(one, "14696_1_Trio", two, "14696_2_Trio", three, "14696_3_Trio", product=product, interval=interval, on="normalized_read_count", T) - - one = triplets[triplets$Sample == "16278_Left",] - two = triplets[triplets$Sample == "26402_Left",] - three = triplets[triplets$Sample == "26759_Left",] - tripletAnalysis(one, "16278_Left_Trio", two, "26402_Left_Trio", three, "26759_Left_Trio", product=product, interval=interval, on="normalized_read_count", T) - - one = triplets[triplets$Sample == "16278_Right",] - two = triplets[triplets$Sample == "26402_Right",] - three = triplets[triplets$Sample == "26759_Right",] - tripletAnalysis(one, "16278_Right_Trio", two, "26402_Right_Trio", three, "26759_Right_Trio", product=product, interval=interval, on="normalized_read_count", T) - - cat("<tr><td>Starting Frequency analysis</td></tr>", file=logfile, append=T) + interval = intervalReads + intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval)) + 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))) - interval = intervalFreq - intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval)) - 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))) - - one = triplets[triplets$Sample == "14696_reg_BM",] - two = triplets[triplets$Sample == "24536_reg_BM",] - three = triplets[triplets$Sample == "24062_reg_BM",] - tripletAnalysis(one, "14696_1_Trio", two, "14696_2_Trio", three, "14696_3_Trio", product=product, interval=interval, on="Frequency", F) - - one = triplets[triplets$Sample == "16278_Left",] - two = triplets[triplets$Sample == "26402_Left",] - three = triplets[triplets$Sample == "26759_Left",] - tripletAnalysis(one, "16278_Left_Trio", two, "26402_Left_Trio", three, "26759_Left_Trio", product=product, interval=interval, on="Frequency", F) - - one = triplets[triplets$Sample == "16278_Right",] - two = triplets[triplets$Sample == "26402_Right",] - three = triplets[triplets$Sample == "26759_Right",] - tripletAnalysis(one, "16278_Right_Trio", two, "26402_Right_Trio", three, "26759_Right_Trio", product=product, interval=interval, on="Frequency", F) + triplets = split(triplets, triplets$Patient, drop=T) + print(nrow(triplets)) + for(triplet in triplets){ + samples = unique(triplet$Sample) + one = triplet[triplet$Sample == samples[1],] + two = triplet[triplet$Sample == samples[2],] + three = triplet[triplet$Sample == samples[3],] + + print(paste(nrow(triplet), nrow(one), nrow(two), nrow(three))) + tripletAnalysis(one, one[1,"uniqueID"], two, two[1,"uniqueID"], three, three[1,"uniqueID"], product=product, interval=interval, on="normalized_read_count", T) + } + + cat("<tr><td>Starting Frequency analysis</td></tr>", file=logfile, append=T) + + interval = intervalFreq + intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval)) + 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))) + + for(triplet in triplets){ + samples = unique(triplet$Sample) + one = triplet[triplet$Sample == samples[1],] + two = triplet[triplet$Sample == samples[2],] + three = triplet[triplet$Sample == samples[3],] + tripletAnalysis(one, one[1,"uniqueID"], two, two[1,"uniqueID"], three, three[1,"uniqueID"], product=product, interval=interval, on="Frequency", F) + } } else { cat("", file="triplets.txt") }
--- a/wrapper.sh Wed Aug 31 05:31:47 2016 -0400 +++ b/wrapper.sh Tue Jan 17 07:24:44 2017 -0500 @@ -46,7 +46,7 @@ oldLocus="" sample1="$(echo ${sample1} | tr -d '\r' | tr -d '\n')" sample2="$(echo ${sample2} | tr -d '\r' | tr -d '\n')" - tail -n+2 ${patient}_freq.txt | sed "s/>//" > tmp.txt + tail -n+2 "${patient}_freq.txt" | sed "s/>//" > tmp.txt echo "<div class='tabber'>" >> "$html" echo "<div class='tabbertab' title='Data frequency'>" >> "$html" echo "<table><tr><td style='vertical-align:top;'>" >> "$html" @@ -99,7 +99,7 @@ echo "<a href='${patient}_percent_freq.png'><img src='${patient}_percent_freq.png' width='1280' height='720' /></a></div>" >> "$html" echo "${scatterplot_tab}</tr></table></div>" >> "$html" - tail -n+2 ${patient}_reads.txt | sed "s/>//" > tmp.txt + tail -n+2 "${patient}_reads.txt" | sed "s/>//" > tmp.txt echo "<div class='tabbertab' title='Data reads'>" >> "$html" echo "<table><tr><td style='vertical-align:top;'>" >> "$html" echo "<table border = 1 class='result_table summary_table' id='summary_table_${patient}_reads'>" >> "$html" @@ -189,9 +189,9 @@ echo "$patient" html="${patient}.html" echo "<tr><td><a href='${patient}.html'>$patient</a></td></tr>" >> "index.html" - echo "$header" > $html + echo "$header" > "$html" oldLocus="" - tail -n+2 ${patient}_freq.txt | sed "s/>//" > tmp.txt + tail -n+2 "${patient}_freq.txt" | sed "s/>//" > tmp.txt echo "<div class='tabber'>" >> "$html" echo "<div class='tabbertab' title='Data frequency'>" >> "$html" echo "<table><tr><td style='vertical-align:top;'>" >> "$html" @@ -261,7 +261,7 @@ echo "<a href='${patient}_freq_indiv_all.png'><img src='${patient}_freq_indiv_all.png' width='1280' height='720' /></a><br /></div>" >> "$html" echo "${scatterplot_tab}</tr></table></div>" >> "$html" - tail -n+2 ${patient}_reads.txt | sed "s/>//" > tmp.txt + tail -n+2 "${patient}_reads.txt" | sed "s/>//" > tmp.txt echo "<div class='tabbertab' title='Data reads'>" >> "$html" echo "<table><tr><td style='vertical-align:top;'>" >> "$html" echo "<table border = 1 class='result_table summary_table' id='summary_table_${patient}_reads'>" >> "$html"