@@ -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.
diff -r ed6885c85660 -r 75853bceec00 RScript.r
--- 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("| Normalizing to lowest cell count within locus |
", file=logfile, append=T)
@@ -475,21 +487,21 @@
print(plt)
dev.off()
}
-
-cat("| Starting Frequency analysis |
", file=logfile, append=T)
+if(length(patients) > 0){
+ cat("| Starting Frequency analysis |
", 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("| Starting Cell Count analysis |
", file=logfile, append=T)
+ cat("| Starting Cell Count analysis |
", 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("| ", label1, " | ", sep=""), file=logfile, append=T)
+ cat(paste("
| ", label1, " | ", 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("", nrow(patient1), " in ", label1, " and ", nrow(patient2), " in ", label2, nrow(patient3), " in ", label3, ", ", nrow(patientMerge), " in both (fetched from cache) |
", 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("", nrow(patient1), " in ", label1, " and ", nrow(patient2), " in ", label2, nrow(patient3), " in ", label3, ", ", nrow(patientMerge), " in both (fetched from cache) | ", 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("| ", label1, " row ", 1:nrow(tmp.rows), " | ", tmp.rows$Sample, ": | ", tmp.rows$Clone_Sequence, " | ", tmp.rows$normalized_read_count, " |
", 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("| ", label1, " row ", 1:nrow(tmp.rows), " | ", tmp.rows$Sample, ": | ", tmp.rows$Clone_Sequence, " | ", tmp.rows$normalized_read_count, " |
", 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("", 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) | ", 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("", 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) | ", 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("| Starting triplet analysis |
", file=logfile, append=T)
+ cat("| Starting triplet analysis |
", file=logfile, append=T)
+
+ triplets$uniqueID = paste(triplets$Patient, triplets$Sample, sep="_")
+
+ cat("| Normalizing to lowest cell count within locus |
", 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("| Normalizing to lowest cell count within locus |
", 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("| Starting Cell Count analysis |
", 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("| Starting Cell Count analysis |
", 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("| Starting Frequency analysis |
", 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("| Starting Frequency analysis |
", 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")
}
diff -r ed6885c85660 -r 75853bceec00 wrapper.sh
--- 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 "" >> "$html"
echo "
" >> "$html"
echo "
" >> "$html"
@@ -99,7 +99,7 @@
echo " " >> "$html"
echo "${scatterplot_tab} |
" >> "$html"
- tail -n+2 ${patient}_reads.txt | sed "s/>//" > tmp.txt
+ tail -n+2 "${patient}_reads.txt" | sed "s/>//" > tmp.txt
echo "
" >> "$html"
echo "
" >> "$html"
echo "" >> "$html"
@@ -189,9 +189,9 @@
echo "$patient"
html="${patient}.html"
echo "| $patient | " >> "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 "" >> "$html"
echo " " >> "$html"
echo " " >> "$html"
@@ -261,7 +261,7 @@
echo "
" >> "$html"
echo "${scatterplot_tab} |
" >> "$html"
- tail -n+2 ${patient}_reads.txt | sed "s/>//" > tmp.txt
+ tail -n+2 "${patient}_reads.txt" | sed "s/>//" > tmp.txt
echo " |