# HG changeset patch
# User davidvanzessen
# Date 1472635907 14400
# Node ID ed6885c85660e10bb1bc3e9b573f5a9c5b43bcb0
Uploaded
diff -r 000000000000 -r ed6885c85660 ALL.xml
--- /dev/null Thu Jan 01 00:00:00 1970 +0000
+++ b/ALL.xml Wed Aug 31 05:31:47 2016 -0400
@@ -0,0 +1,46 @@
+
Starting analysis | |||||||||||||||||||||||||||||||||||||||||||
Reading input | |||||||||||||||||||||||||||||||||||||||||||
Selecting first V/J Genes | |||||||||||||||||||||||||||||||||||||||||||
Calculating Frequency | |||||||||||||||||||||||||||||||||||||||||||
Normalizing to lowest cell count within locus |
", patient, " | ", sep=""), file=logfile, append=T) + + if(mergeOn == "Clone_Sequence"){ + patient1$merge = paste(patient1$Clone_Sequence) + patient2$merge = paste(patient2$Clone_Sequence) + } 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) + } + + scatterplot_data_columns = c("Patient", "Sample", "Frequency", "normalized_read_count", "V_Segment_Major_Gene", "J_Segment_Major_Gene", "merge") + #scatterplot_data = rbind(patient1[,scatterplot_data_columns], patient2[,scatterplot_data_columns]) + scatterplot_data = patient1[NULL,scatterplot_data_columns] + #scatterplot_data = scatterplot_data[!duplicated(scatterplot_data$merge),] + #scatterplot_data$type = factor(x=oneSample, levels=c(oneSample, twoSample, "In Both")) + scatterplot.data.type.factor = c(oneSample, twoSample, paste(c(oneSample, twoSample), "In Both")) + #scatterplot_data$type = factor(x=NULL, levels=scatterplot.data.type.factor) + scatterplot_data$type = character(0) + scatterplot_data$link = numeric(0) + scatterplot_data$on = character(0) + + #patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge") #merge alles 'fuzzy' + patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge")[NULL,] #blegh + + cs.exact.matches = patient1[patient1$Clone_Sequence %in% patient2$Clone_Sequence,]$Clone_Sequence + + start.time = proc.time() + merge.list = c() + + if(patient %in% names(patient.merge.list)){ + patientMerge = patient.merge.list[[patient]] + merge.list[["second"]] = patient.merge.list.second[[patient]] + scatterplot_data = scatter_locus_data_list[[patient]] + cat(paste("", nrow(patient1), " in ", oneSample, " and ", nrow(patient2), " in ", twoSample, ", ", nrow(patientMerge), " in both (fetched from cache) | ||
", patient, " row ", 1:nrow(tmp.rows), " | ", tmp.rows$Sample, ": | ", tmp.rows$Clone_Sequence, " | ", tmp.rows$normalized_read_count, " |
", patient, " row ", 1:nrow(tmp.rows), " | ", tmp.rows$Sample, ": | ", tmp.rows$Clone_Sequence, " | ", tmp.rows$normalized_read_count, " |
", patient, " row ", 1:nrow(tmp.rows), " | ", tmp.rows$Sample, ": | ", tmp.rows$Clone_Sequence, " | ", tmp.rows$normalized_read_count, " |
", patient, " row ", 1:nrow(first.rows), " | ", first.rows$Sample, ": | ", first.rows$Clone_Sequence, " | ", first.rows$normalized_read_count, " | ", nrow(patient1), " in ", oneSample, " and ", nrow(patient2), " in ", twoSample, ", ", nrow(patientMerge), " in both (finding both took ", (proc.time() - start.time)[[3]], "s) | ", sep=""), file=logfile, append=T) + } + + patient1 = patient1[!(patient1$Clone_Sequence %in% patient.merge.list.second[[patient]]),] + patient2 = patient2[!(patient2$Clone_Sequence %in% patient.merge.list.second[[patient]]),] + + + patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony]) + #patientMerge$thresholdValue = pmin(patientMerge[,onx], patientMerge[,ony]) + res1 = vector() + res2 = vector() + resBoth = vector() + read1Count = vector() + read2Count = vector() + locussum1 = vector() + locussum2 = vector() + + #for(iter in 1){ + 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="") + #both = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge[,onx] > threshhold & patientMerge[,ony] > threshhold) #both higher than threshold + both = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge$thresholdValue > threshhold) #highest of both is higher than threshold + one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$merge %in% patientMerge[both,]$merge)) + two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$merge %in% patientMerge[both,]$merge)) + read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count)) + read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count)) + res1 = append(res1, sum(one)) + res2 = append(res2, sum(two)) + resBoth = append(resBoth, sum(both)) + locussum1 = append(locussum1, sum(patient1[(grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene)),]$normalized_read_count)) + locussum2 = append(locussum2, sum(patient2[(grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene)),]$normalized_read_count)) + #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(oneSample, "_", 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(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="") + write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) + } + } else { + scatterplot_locus_data = scatterplot_data[grepl(V_Segment, scatterplot_data$V_Segment_Major_Gene) & grepl(J_Segment, scatterplot_data$J_Segment_Major_Gene),] + if(nrow(scatterplot_locus_data) > 0){ + scatterplot_locus_data$Rearrangement = product[iter, titleIndex] + } + + + + p = NULL + print(paste("nrow scatterplot_locus_data", nrow(scatterplot_locus_data))) + if(nrow(scatterplot_locus_data) != 0){ + if(on == "normalized_read_count"){ + write.table(scatterplot_locus_data, file=paste(oneSample, twoSample, product[iter, titleIndex], "scatterplot_locus_data.txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) + scales = 10^(0:6) #(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count)))) + p = ggplot(scatterplot_locus_data, aes(factor(reorder(type, type.order)), normalized_read_count, group=link)) + geom_line() + scale_y_log10(breaks=scales,labels=scales, limits=c(1,1e6)) + scale_x_discrete(breaks=levels(scatterplot_data$type), labels=levels(scatterplot_data$type), drop=FALSE) + } else { + p = ggplot(scatterplot_locus_data, aes(factor(reorder(type, type.order)), Frequency, group=link)) + geom_line() + scale_y_log10(limits=c(0.0001,100), breaks=c(0.0001, 0.001, 0.01, 0.1, 1, 10, 100), labels=c("0.0001", "0.001", "0.01", "0.1", "1", "10", "100")) + scale_x_discrete(breaks=levels(scatterplot_data$type), labels=levels(scatterplot_data$type), drop=FALSE) + } + p = p + geom_point(aes(colour=type), position="dodge") + p = p + xlab("In one or both samples") + ylab(onShort) + ggtitle(paste(patient1[1,patientIndex], patient1[1,sampleIndex], patient2[1,sampleIndex], onShort, product[iter, titleIndex])) + } else { + p = ggplot(NULL, aes(x=c("In one", "In Both"),y=0)) + geom_blank(NULL) + xlab("In one or both of the samples") + ylab(onShort) + ggtitle(paste(patient1[1,patientIndex], patient1[1,sampleIndex], patient2[1,sampleIndex], onShort, product[iter, titleIndex])) + } + png(paste(patient1[1,patientIndex], "_", patient1[1,sampleIndex], "_", patient2[1,sampleIndex], "_", onShort, "_", product[iter, titleIndex],"_scatter.png", sep="")) + print(p) + dev.off() + } + if(sum(both) > 0){ + dfBoth = patientMerge[both,c("V_Segment_Major_Gene.x", "J_Segment_Major_Gene.x", "normalized_read_count.x", "Frequency.x", "Related_to_leukemia_clone.x", "Clone_Sequence.x", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "normalized_read_count.y", "Frequency.y", "Related_to_leukemia_clone.y")] + colnames(dfBoth) = c(paste("Proximal segment", oneSample), paste("Distal segment", oneSample), paste("Normalized_Read_Count", oneSample), paste("Frequency", oneSample), paste("Related_to_leukemia_clone", oneSample),"Clone Sequence", paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Normalized_Read_Count", twoSample), paste("Frequency", twoSample), paste("Related_to_leukemia_clone", twoSample)) + filenameBoth = paste(oneSample, "_", twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="") + write.table(dfBoth, file=paste(filenameBoth, ".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=""), "Both"=resBoth, "tmp1"=res1, "read_count1" = round(read1Count), "tmp2"=res2, "read_count2"= round(read2Count), "Sum"=res1 + res2 + resBoth, "percentage" = round((resBoth/(res1 + res2 + resBoth)) * 100, digits=2), "Locus_sum1"=locussum1, "Locus_sum2"=locussum2) + if(sum(is.na(patientResult$percentage)) > 0){ + patientResult[is.na(patientResult$percentage),]$percentage = 0 + } + colnames(patientResult)[6] = oneSample + colnames(patientResult)[8] = twoSample + colnamesBak = colnames(patientResult) + colnames(patientResult) = c("Ig/TCR gene rearrangement type", "Distal Gene segment", "Proximal gene segment", "cut_off_value", paste("Number of sequences ", patient, "_Both", sep=""), paste("Number of sequences", oneSample, sep=""), paste("Normalized Read Count", oneSample), paste("Number of sequences", twoSample, sep=""), paste("Normalized Read Count", twoSample), paste("Sum number of sequences", patient), paste("Percentage of sequences ", patient, "_Both", sep=""), paste("Locus Sum", oneSample), paste("Locus Sum", twoSample)) + write.table(patientResult, file=paste(patient, "_", 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", "Both")]) + plt = plt + geom_bar( aes( x=factor(cut_off_value), y=Both), 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(Both), x=cut_off_value,y=Both,label=Both), angle=90, hjust=0) + plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in both") + plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines")) + png(paste(patient, "_", onShort, ".png", sep=""), width=1920, height=1080) + print(plt) + dev.off() + #(t,r,b,l) + plt = ggplot(patientResult[,c("Locus", "cut_off_value", "percentage")]) + plt = plt + geom_bar( aes( x=factor(cut_off_value), y=percentage), 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(percentage), x=cut_off_value,y=percentage,label=percentage), angle=90, hjust=0) + plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("% clones in both left and right") + plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines")) + png(paste(patient, "_percent_", onShort, ".png", sep=""), width=1920, height=1080) + print(plt) + dev.off() + + patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample)] ,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.2) + plt = plt + geom_text(data=patientResult[patientResult$variable == twoSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=0.8) + plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle(paste("Number of clones in only ", oneSample, " and only ", twoSample, sep="")) + png(paste(patient, "_", onShort, "_both.png", sep=""), width=1920, height=1080) + print(plt) + dev.off() +} + +cat("
Starting Frequency analysis | |||
Starting Cell Count analysis | |||
", 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) + + 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) + + 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() + + 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 + + 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() + + 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]] + + cat(paste("", nrow(patient1), " in ", label1, " and ", nrow(patient2), " in ", label2, nrow(patient3), " in ", label3, ", ", nrow(patientMerge), " in both (fetched from cache) | ||
", label1, " row ", 1:nrow(tmp.rows), " | ", tmp.rows$Sample, ": | ", tmp.rows$Clone_Sequence, " | ", tmp.rows$normalized_read_count, " | ", 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]) + + 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"]]),] + + 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() +} + +if(nrow(triplets) != 0){ + + cat("
Starting triplet analysis | |||
Normalizing to lowest cell count within locus | |||
Starting Cell Count analysis | |||
Starting Frequency analysis |
t |
" >> "$html"
+ echo "
|
" >> "$html"
+ echo "
|
Singles (Frequency scatterplot, Reads scatterplot): | |||||||||||||||||||||||||||||||||||||||||||
$patient | |||||||||||||||||||||||||||||||||||||||||||
Pairs (Left & Right): | |||||||||||||||||||||||||||||||||||||||||||
$patient | |||||||||||||||||||||||||||||||||||||||||||
Pairs (BM & PB): | |||||||||||||||||||||||||||||||||||||||||||
$patient | |||||||||||||||||||||||||||||||||||||||||||
Pairs (Dx & R): | |||||||||||||||||||||||||||||||||||||||||||
$patient | |||||||||||||||||||||||||||||||||||||||||||
Triplets: | |||||||||||||||||||||||||||||||||||||||||||
$patient |
" >> "$html"
+ echo "
|
" >> "$html"
+ echo "
|