|
0
|
1 args <- commandArgs(trailingOnly = TRUE)
|
|
|
2
|
|
|
3 inFile = args[1]
|
|
|
4 outDir = args[2]
|
|
3
|
5 logfile = args[3]
|
|
|
6 min_freq = as.numeric(args[4])
|
|
|
7 min_cells = as.numeric(args[5])
|
|
|
8
|
|
|
9 cat("<html><table><tr><td>Starting analysis</td></tr>", file=logfile, append=F)
|
|
0
|
10
|
|
4
|
11 library(ggplot2)
|
|
|
12 library(reshape2)
|
|
|
13 library(data.table)
|
|
|
14 library(grid)
|
|
|
15 library(parallel)
|
|
0
|
16 #require(xtable)
|
|
3
|
17 cat("<tr><td>Reading input</td></tr>", file=logfile, append=T)
|
|
9
|
18 dat = read.table(inFile, header=T, sep="\t", dec=",", fill=T, stringsAsFactors=F)
|
|
|
19 dat = dat[!is.na(dat$Patient),]
|
|
12
|
20 dat = dat[!duplicated(dat$Clone_Sequence), ]
|
|
9
|
21
|
|
0
|
22 setwd(outDir)
|
|
3
|
23 cat("<tr><td>Selecting first V/J Genes</td></tr>", file=logfile, append=T)
|
|
2
|
24 dat$V_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$V_Segment_Major_Gene), "; "), "[[", 1)))
|
|
|
25 dat$J_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$J_Segment_Major_Gene), "; "), "[[", 1)))
|
|
|
26
|
|
12
|
27 dat$Frequency = ((10^dat$Log10_Frequency)*100)
|
|
|
28
|
|
9
|
29 cat("<tr><td>Deduplication</td></tr>", file=logfile, append=T)
|
|
11
|
30 #dat = data.frame(data.table(dat)[, list(Patient=unique(.SD$Patient), Clone_Molecule_Count_From_Spikes=sum(.SD$Clone_Molecule_Count_From_Spikes), Log10_Frequency=sum(.SD$Log10_Frequency), Total_Read_Count=sum(.SD$Total_Read_Count), Related_to_leukemia_clone=any(.SD$Related_to_leukemia_clone)), by=c("Sample", "Cell_Count", "J_Segment_Major_Gene", "V_Segment_Major_Gene", "CDR3_Sense_Sequence")])
|
|
|
31
|
|
12
|
32 most.common = function(x, ret="V"){
|
|
|
33 past = paste(x$V_Segment_Major_Gene, x$J_Segment_Major_Gene, sep=";")
|
|
|
34 ux = unique(past)
|
|
11
|
35 if(length(ux) > 1){
|
|
12
|
36 xtdf = data.frame(table(past))
|
|
11
|
37 #print(xtdf)
|
|
12
|
38 res = unlist(strsplit(as.character(xtdf$past[which.max(xtdf$Freq)]), ";"))
|
|
|
39 #print(res)
|
|
|
40 if(ret == "V"){
|
|
|
41 return(res[1])
|
|
|
42 } else {
|
|
|
43 return(res[2])
|
|
|
44 }
|
|
11
|
45 }
|
|
12
|
46
|
|
|
47 if(ret == "V"){
|
|
|
48 return(unique(x$V_Segment_Major_Gene))
|
|
|
49 } else {
|
|
|
50 return(unique(x$J_Segment_Major_Gene))
|
|
|
51 }
|
|
11
|
52 }
|
|
|
53
|
|
12
|
54 dat = data.frame(data.table(dat)[, list(Patient=unique(.SD$Patient), V_Segment_Major_Gene= as.character(most.common(.SD, ret="V")), J_Segment_Major_Gene= as.character(most.common(.SD, ret="J")), Clone_Molecule_Count_From_Spikes=sum(.SD$Clone_Molecule_Count_From_Spikes), Log10_Frequency=sum(.SD$Log10_Frequency), Frequency=sum(.SD$Frequency), Total_Read_Count=sum(.SD$Total_Read_Count), Related_to_leukemia_clone=any(.SD$Related_to_leukemia_clone)), by=c("Sample", "Cell_Count", "CDR3_Sense_Sequence")])
|
|
|
55 dat = data.frame(data.table(dat)[, list(Patient=unique(.SD$Patient), Clone_Molecule_Count_From_Spikes=sum(.SD$Clone_Molecule_Count_From_Spikes), Log10_Frequency=sum(.SD$Log10_Frequency), Frequency=sum(.SD$Frequency), Total_Read_Count=sum(.SD$Total_Read_Count), Related_to_leukemia_clone=any(.SD$Related_to_leukemia_clone)), by=c("Sample", "Cell_Count", "J_Segment_Major_Gene", "V_Segment_Major_Gene", "CDR3_Sense_Sequence")])
|
|
9
|
56
|
|
3
|
57 cat("<tr><td>Calculating Frequency</td></tr>", file=logfile, append=T)
|
|
12
|
58
|
|
2
|
59
|
|
3
|
60 dat = dat[dat$Frequency >= min_freq,]
|
|
|
61
|
|
12
|
62 #cat("<tr><td>Normalizing cell count to 1.000.000</td></tr>", file=logfile, append=T)
|
|
|
63 #dat$normalized_read_count = round(dat$Clone_Molecule_Count_From_Spikes / dat$Cell_Count * 1000000 / 2)
|
|
|
64 dat$normalized_read_count = dat$Clone_Molecule_Count_From_Spikes
|
|
3
|
65 dat = dat[dat$normalized_read_count >= min_cells,]
|
|
2
|
66 dat$paste = paste(dat$Sample, dat$V_Segment_Major_Gene, dat$J_Segment_Major_Gene, dat$CDR3_Sense_Sequence)
|
|
9
|
67 triplets = dat[grepl("VanDongen_cALL_14696", dat$Patient) | grepl("(16278)|(26402)|(26759)", dat$Sample),]
|
|
|
68
|
|
0
|
69 patients = split(dat, dat$Patient, drop=T)
|
|
9
|
70 intervalReads = rev(c(0,10,25,50,100,250,500,750,1000,10000))
|
|
6
|
71 intervalFreq = rev(c(0,0.01,0.05,0.1,0.5,1,5))
|
|
0
|
72 V_Segments = c(".*", "IGHV", "IGHD", "IGKV", "IGKV", "IgKINTR", "TRGV", "TRDV", "TRDD" , "TRBV")
|
|
|
73 J_Segments = c(".*", ".*", ".*", "IGKJ", "KDE", ".*", ".*", ".*", ".*", ".*")
|
|
|
74 Titles = c("Total", "IGH-Vh-Jh", "IGH-Dh-Jh", "Vk-Jk", "Vk-Kde" , "Intron-Kde", "TCRG", "TCRD-Vd-Dd", "TCRD-Dd-Dd", "TCRB-Vb-Jb")
|
|
|
75 Titles = factor(Titles, levels=Titles)
|
|
|
76 TitlesOrder = data.frame("Title"=Titles, "TitlesOrder"=1:length(Titles))
|
|
|
77
|
|
|
78 patientCountOnColumn <- function(x, product, interval, on, appendtxt=F){
|
|
|
79 x$Sample = factor(x$Sample, levels=unique(x$Sample))
|
|
|
80 onShort = "reads"
|
|
|
81 if(on == "Frequency"){
|
|
|
82 onShort = "freq"
|
|
|
83 }
|
|
|
84 splt = split(x, x$Sample, drop=T)
|
|
4
|
85 type="pair"
|
|
2
|
86 if(length(splt) == 1){
|
|
3
|
87 print(paste(paste(x[1,which(colnames(x) == "Patient")]), "has one sample"))
|
|
4
|
88 splt[[2]] = data.frame("Patient" = character(0), "Receptor" = character(0), "Sample" = character(0), "Cell_Count" = numeric(0), "Clone_Molecule_Count_From_Spikes" = numeric(0), "Log10_Frequency" = numeric(0), "Total_Read_Count" = numeric(0), "dsMol_per_1e6_cells" = numeric(0), "J_Segment_Major_Gene" = character(0), "V_Segment_Major_Gene" = character(0), "Clone_Sequence" = character(0), "CDR3_Sense_Sequence" = character(0), "Related_to_leukemia_clone" = logical(0), "Frequency"= numeric(0), "normalized_read_count" = numeric(0), "paste" = character(0))
|
|
|
89 type="single"
|
|
2
|
90 }
|
|
0
|
91 patient1 = splt[[1]]
|
|
|
92 patient2 = splt[[2]]
|
|
|
93
|
|
|
94 threshholdIndex = which(colnames(product) == "interval")
|
|
|
95 V_SegmentIndex = which(colnames(product) == "V_Segments")
|
|
|
96 J_SegmentIndex = which(colnames(product) == "J_Segments")
|
|
|
97 titleIndex = which(colnames(product) == "Titles")
|
|
|
98 sampleIndex = which(colnames(x) == "Sample")
|
|
|
99 patientIndex = which(colnames(x) == "Patient")
|
|
|
100 oneSample = paste(patient1[1,sampleIndex], sep="")
|
|
|
101 twoSample = paste(patient2[1,sampleIndex], sep="")
|
|
|
102 patient = paste(x[1,patientIndex])
|
|
3
|
103
|
|
0
|
104 switched = F
|
|
|
105 if(length(grep(".*_Right$", twoSample)) == 1 || length(grep(".*_Dx_BM$", twoSample)) == 1 || length(grep(".*_Dx$", twoSample)) == 1 ){
|
|
|
106 tmp = twoSample
|
|
|
107 twoSample = oneSample
|
|
|
108 oneSample = tmp
|
|
|
109 tmp = patient1
|
|
|
110 patient1 = patient2
|
|
|
111 patient2 = tmp
|
|
|
112 switched = T
|
|
|
113 }
|
|
|
114 if(appendtxt){
|
|
4
|
115 cat(paste(patient, oneSample, twoSample, type, sep="\t"), file="patients.txt", append=T, sep="", fill=3)
|
|
0
|
116 }
|
|
3
|
117 cat(paste("<tr><td>", patient, "</td></tr>", sep=""), file=logfile, append=T)
|
|
9
|
118
|
|
12
|
119 #patient1$merge = paste(patient1$V_Segment_Major_Gene, patient1$J_Segment_Major_Gene, patient1$CDR3_Sense_Sequence)
|
|
|
120 #patient2$merge = paste(patient2$V_Segment_Major_Gene, patient2$J_Segment_Major_Gene, patient2$CDR3_Sense_Sequence)
|
|
|
121 patient1$merge = paste(patient1$CDR3_Sense_Sequence)
|
|
|
122 patient2$merge = paste(patient2$CDR3_Sense_Sequence)
|
|
9
|
123
|
|
12
|
124 #patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge")
|
|
9
|
125 patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge")
|
|
0
|
126 res1 = vector()
|
|
|
127 res2 = vector()
|
|
|
128 resBoth = vector()
|
|
|
129 read1Count = vector()
|
|
|
130 read2Count = vector()
|
|
2
|
131 locussum1 = vector()
|
|
|
132 locussum2 = vector()
|
|
9
|
133
|
|
|
134 print(patient)
|
|
0
|
135 #for(iter in 1){
|
|
|
136 for(iter in 1:length(product[,1])){
|
|
|
137 threshhold = product[iter,threshholdIndex]
|
|
|
138 V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="")
|
|
|
139 J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="")
|
|
|
140 both = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge[,paste(on, ".x", sep="")] > threshhold & patientMerge[,paste(on, ".y", sep="")] > threshhold)
|
|
10
|
141 one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$CDR3_Sense_Sequence %in% patientMerge[both,]$CDR3_Sense_Sequence.x))
|
|
|
142 two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$CDR3_Sense_Sequence %in% patientMerge[both,]$CDR3_Sense_Sequence.x))
|
|
2
|
143 read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count) + sum(patientMerge[both,]$normalized_read_count.x))
|
|
|
144 read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count) + sum(patientMerge[both,]$normalized_read_count.y))
|
|
0
|
145 res1 = append(res1, sum(one))
|
|
2
|
146 res2 = append(res2, sum(two))
|
|
0
|
147 resBoth = append(resBoth, sum(both))
|
|
2
|
148 locussum1 = append(locussum1, sum(patient1[(grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene)),]$normalized_read_count))
|
|
|
149 locussum2 = append(locussum2, sum(patient2[(grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene)),]$normalized_read_count))
|
|
0
|
150 #threshhold = 0
|
|
|
151 if(threshhold != 0){
|
|
|
152 if(sum(one) > 0){
|
|
9
|
153 dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "CDR3_Sense_Sequence", "Related_to_leukemia_clone")]
|
|
|
154 colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "CDR3 Sequence", "Related_to_leukemia_clone")
|
|
0
|
155 filenameOne = paste(oneSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
156 write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
157 }
|
|
|
158 if(sum(two) > 0){
|
|
9
|
159 dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "CDR3_Sense_Sequence", "Related_to_leukemia_clone")]
|
|
|
160 colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "CDR3 Sequence", "Related_to_leukemia_clone")
|
|
0
|
161 filenameTwo = paste(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
162 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
163 }
|
|
|
164 }
|
|
|
165 if(sum(both) > 0){
|
|
9
|
166 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", "CDR3_Sense_Sequence.x", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "normalized_read_count.y", "Frequency.y", "Related_to_leukemia_clone.y")]
|
|
|
167 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),"CDR3 Sequence", paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Normalized_Read_Count", twoSample), paste("Frequency", twoSample), paste("Related_to_leukemia_clone", twoSample))
|
|
0
|
168 filenameBoth = paste(oneSample, "_", twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
169 write.table(dfBoth, file=paste(filenameBoth, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
170 }
|
|
|
171 }
|
|
2
|
172 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)
|
|
0
|
173 if(sum(is.na(patientResult$percentage)) > 0){
|
|
|
174 patientResult[is.na(patientResult$percentage),]$percentage = 0
|
|
|
175 }
|
|
|
176 colnames(patientResult)[6] = oneSample
|
|
|
177 colnames(patientResult)[8] = twoSample
|
|
|
178 colnamesBak = colnames(patientResult)
|
|
2
|
179 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))
|
|
0
|
180 write.table(patientResult, file=paste(patient, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
181 colnames(patientResult) = colnamesBak
|
|
|
182
|
|
|
183 patientResult$Locus = factor(patientResult$Locus, Titles)
|
|
|
184 patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))
|
|
|
185
|
|
|
186 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "Both")])
|
|
|
187 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=Both), stat='identity', position="dodge", fill="#79c36a")
|
|
|
188 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
189 plt = plt + geom_text(aes(ymax=max(Both), x=cut_off_value,y=Both,label=Both), angle=90, hjust=0)
|
|
|
190 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in both")
|
|
|
191 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
|
192 png(paste(patient, "_", onShort, ".png", sep=""), width=1920, height=1080)
|
|
|
193 print(plt)
|
|
|
194 dev.off()
|
|
|
195 #(t,r,b,l)
|
|
|
196 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "percentage")])
|
|
|
197 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=percentage), stat='identity', position="dodge", fill="#79c36a")
|
|
|
198 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
199 plt = plt + geom_text(aes(ymax=max(percentage), x=cut_off_value,y=percentage,label=percentage), angle=90, hjust=0)
|
|
|
200 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("% clones in both left and right")
|
|
|
201 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
|
202 png(paste(patient, "_percent_", onShort, ".png", sep=""), width=1920, height=1080)
|
|
|
203 print(plt)
|
|
|
204 dev.off()
|
|
|
205
|
|
|
206 patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample)] ,id.vars=1:2)
|
|
|
207 patientResult$relativeValue = patientResult$value * 10
|
|
|
208 patientResult[patientResult$relativeValue == 0,]$relativeValue = 1
|
|
|
209 plt = ggplot(patientResult)
|
|
|
210 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge")
|
|
|
211 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
212 plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9)))
|
|
|
213 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)
|
|
|
214 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)
|
|
|
215 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle(paste("Number of clones in only ", oneSample, " and only ", twoSample, sep=""))
|
|
|
216 png(paste(patient, "_", onShort, "_both.png", sep=""), width=1920, height=1080)
|
|
|
217 print(plt)
|
|
|
218 dev.off()
|
|
|
219 }
|
|
|
220
|
|
3
|
221 cat("<tr><td>Starting Frequency analysis</td></tr>", file=logfile, append=T)
|
|
|
222
|
|
0
|
223 interval = intervalFreq
|
|
|
224 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
4
|
225 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)))
|
|
|
226 mclapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
|
|
0
|
227
|
|
3
|
228 cat("<tr><td>Starting Cell Count analysis</td></tr>", file=logfile, append=T)
|
|
|
229
|
|
0
|
230 interval = intervalReads
|
|
|
231 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
4
|
232 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)))
|
|
9
|
233 mclapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="normalized_read_count")
|
|
0
|
234
|
|
3
|
235 cat("</table></html>", file=logfile, append=T)
|
|
|
236
|
|
7
|
237
|
|
|
238 tripletAnalysis <- function(patient1, label1, patient2, label2, patient3, label3, product, interval, on, appendTriplets= FALSE){
|
|
|
239 onShort = "reads"
|
|
|
240 if(on == "Frequency"){
|
|
|
241 onShort = "freq"
|
|
|
242 }
|
|
|
243 type="triplet"
|
|
|
244
|
|
|
245 threshholdIndex = which(colnames(product) == "interval")
|
|
|
246 V_SegmentIndex = which(colnames(product) == "V_Segments")
|
|
|
247 J_SegmentIndex = which(colnames(product) == "J_Segments")
|
|
|
248 titleIndex = which(colnames(product) == "Titles")
|
|
|
249 sampleIndex = which(colnames(patient1) == "Sample")
|
|
|
250 patientIndex = which(colnames(patient1) == "Patient")
|
|
|
251 oneSample = paste(patient1[1,sampleIndex], sep="")
|
|
|
252 twoSample = paste(patient2[1,sampleIndex], sep="")
|
|
|
253 threeSample = paste(patient3[1,sampleIndex], sep="")
|
|
|
254
|
|
12
|
255 #patient1$merge = paste(patient1$V_Segment_Major_Gene, patient1$J_Segment_Major_Gene, patient1$CDR3_Sense_Sequence)
|
|
|
256 #patient2$merge = paste(patient2$V_Segment_Major_Gene, patient2$J_Segment_Major_Gene, patient2$CDR3_Sense_Sequence)
|
|
|
257 #patient3$merge = paste(patient3$V_Segment_Major_Gene, patient3$J_Segment_Major_Gene, patient3$CDR3_Sense_Sequence)
|
|
|
258
|
|
|
259 patient1$merge = paste(patient1$CDR3_Sense_Sequence)
|
|
|
260 patient2$merge = paste(patient2$CDR3_Sense_Sequence)
|
|
|
261 patient3$merge = paste(patient3$CDR3_Sense_Sequence)
|
|
9
|
262
|
|
|
263 patientMerge = merge(patient1, patient2, by="merge")
|
|
|
264 patientMerge = merge(patientMerge, patient3, by="merge")
|
|
|
265 colnames(patientMerge)[28:length(colnames(patientMerge))] = paste(colnames(patientMerge)[28:length(colnames(patientMerge))], ".z", sep="")
|
|
7
|
266 res1 = vector()
|
|
|
267 res2 = vector()
|
|
|
268 res3 = vector()
|
|
|
269 resAll = vector()
|
|
|
270 read1Count = vector()
|
|
|
271 read2Count = vector()
|
|
|
272 read3Count = vector()
|
|
|
273
|
|
|
274 if(appendTriplets){
|
|
9
|
275 cat(paste(label1, label2, label3, sep="\t"), file="triplets.txt", append=T, sep="", fill=3)
|
|
7
|
276 }
|
|
|
277 for(iter in 1:length(product[,1])){
|
|
|
278 threshhold = product[iter,threshholdIndex]
|
|
|
279 V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="")
|
|
|
280 J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="")
|
|
|
281 all = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge[,paste(on, ".x", sep="")] > threshhold & patientMerge[,paste(on, ".y", sep="")] > threshhold & patientMerge[,paste(on, ".z", sep="")] > threshhold)
|
|
10
|
282 one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$CDR3_Sense_Sequence %in% patientMerge[all,]$CDR3_Sense_Sequence.x))
|
|
|
283 two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$CDR3_Sense_Sequence %in% patientMerge[all,]$CDR3_Sense_Sequence.x))
|
|
|
284 three = (grepl(V_Segment, patient3$V_Segment_Major_Gene) & grepl(J_Segment, patient3$J_Segment_Major_Gene) & patient3[,on] > threshhold & !(patient3$CDR3_Sense_Sequence %in% patientMerge[all,]$CDR3_Sense_Sequence.x))
|
|
7
|
285
|
|
|
286 read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.x))
|
|
|
287 read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.y))
|
|
|
288 read3Count = append(read3Count, sum(patient3[three,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.z))
|
|
|
289 res1 = append(res1, sum(one))
|
|
|
290 res2 = append(res2, sum(two))
|
|
|
291 res3 = append(res3, sum(three))
|
|
|
292 resAll = append(resAll, sum(all))
|
|
|
293 #threshhold = 0
|
|
|
294 if(threshhold != 0){
|
|
|
295 if(sum(one) > 0){
|
|
9
|
296 dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "CDR3_Sense_Sequence", "Related_to_leukemia_clone")]
|
|
12
|
297 colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "CDR3 Sequence", "Related_to_leukemia_clone")
|
|
7
|
298 filenameOne = paste(label1, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
299 write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
300 }
|
|
|
301 if(sum(two) > 0){
|
|
9
|
302 dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "CDR3_Sense_Sequence", "Related_to_leukemia_clone")]
|
|
12
|
303 colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "CDR3 Sequence", "Related_to_leukemia_clone")
|
|
7
|
304 filenameTwo = paste(label2, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
305 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
306 }
|
|
|
307 if(sum(three) > 0){
|
|
9
|
308 dfThree = patient3[three,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "CDR3_Sense_Sequence", "Related_to_leukemia_clone")]
|
|
12
|
309 colnames(dfThree) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "CDR3 Sequence", "Related_to_leukemia_clone")
|
|
7
|
310 filenameThree = paste(label3, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
311 write.table(dfThree, file=paste(filenameThree, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
312 }
|
|
|
313 }
|
|
|
314 if(sum(all) > 0){
|
|
9
|
315 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", "CDR3_Sense_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")]
|
|
|
316 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),"CDR3_Sense_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))
|
|
7
|
317 filenameAll = paste(label1, "_", label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
318 write.table(dfAll, file=paste(filenameAll, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
319 }
|
|
|
320 }
|
|
|
321 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))
|
|
|
322 colnames(patientResult)[6] = oneSample
|
|
|
323 colnames(patientResult)[8] = twoSample
|
|
|
324 colnames(patientResult)[10] = threeSample
|
|
|
325
|
|
|
326 colnamesBak = colnames(patientResult)
|
|
|
327 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("Normalized Read Count", oneSample), paste("Number of sequences", twoSample), paste("Normalized Read Count", twoSample), paste("Number of sequences", threeSample), paste("Normalized Read Count", threeSample))
|
|
|
328 write.table(patientResult, file=paste(label1, "_", label2, "_", label3, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
329 colnames(patientResult) = colnamesBak
|
|
|
330
|
|
|
331 patientResult$Locus = factor(patientResult$Locus, Titles)
|
|
|
332 patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))
|
|
|
333
|
|
|
334 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "All")])
|
|
|
335 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=All), stat='identity', position="dodge", fill="#79c36a")
|
|
|
336 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
337 plt = plt + geom_text(aes(ymax=max(All), x=cut_off_value,y=All,label=All), angle=90, hjust=0)
|
|
|
338 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in All")
|
|
|
339 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
|
340 png(paste(label1, "_", label2, "_", label3, "_", onShort, "_total_all.png", sep=""), width=1920, height=1080)
|
|
|
341 print(plt)
|
|
|
342 dev.off()
|
|
|
343
|
|
|
344 fontSize = 4
|
|
|
345
|
|
|
346 bak = patientResult
|
|
|
347 patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample, threeSample)] ,id.vars=1:2)
|
|
|
348 patientResult$relativeValue = patientResult$value * 10
|
|
|
349 patientResult[patientResult$relativeValue == 0,]$relativeValue = 1
|
|
|
350 plt = ggplot(patientResult)
|
|
|
351 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge")
|
|
|
352 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
353 plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9)))
|
|
|
354 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)
|
|
|
355 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)
|
|
|
356 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)
|
|
|
357 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in only one sample")
|
|
|
358 png(paste(label1, "_", label2, "_", label3, "_", onShort, "_indiv_all.png", sep=""), width=1920, height=1080)
|
|
|
359 print(plt)
|
|
|
360 dev.off()
|
|
|
361 }
|
|
|
362
|
|
9
|
363
|
|
|
364 triplets$uniqueID = "ID"
|
|
|
365
|
|
|
366 triplets[grepl("16278_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
|
|
|
367 triplets[grepl("26402_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
|
|
|
368 triplets[grepl("26759_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
|
|
|
369
|
|
|
370 triplets[grepl("16278_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
|
|
|
371 triplets[grepl("26402_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
|
|
|
372 triplets[grepl("26759_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
|
|
|
373
|
|
|
374 triplets[grepl("14696", triplets$Patient),]$uniqueID = "14696"
|
|
|
375
|
|
|
376 triplets = data.frame(data.table(triplets)[, list(Patient=unique(.SD$uniqueID), Clone_Molecule_Count_From_Spikes=sum(.SD$Clone_Molecule_Count_From_Spikes), Log10_Frequency=sum(.SD$Log10_Frequency), Total_Read_Count=sum(.SD$Total_Read_Count), Related_to_leukemia_clone=any(.SD$Related_to_leukemia_clone)), by=c("Sample", "Cell_Count", "J_Segment_Major_Gene", "V_Segment_Major_Gene", "CDR3_Sense_Sequence")])
|
|
|
377
|
|
|
378 triplets$Frequency = (10^as.numeric(triplets$Log10_Frequency))*100
|
|
|
379 triplets$normalized_read_count = round(triplets$Clone_Molecule_Count_From_Spikes / triplets$Cell_Count * 1000000 / 2)
|
|
|
380
|
|
7
|
381 interval = intervalReads
|
|
|
382 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
|
383 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)))
|
|
|
384
|
|
9
|
385 one = triplets[triplets$Sample == "14696_reg_BM",]
|
|
|
386 two = triplets[triplets$Sample == "24536_reg_BM",]
|
|
|
387 three = triplets[triplets$Sample == "24062_reg_BM",]
|
|
8
|
388 tripletAnalysis(one, "14696_1", two, "14696_2", three, "14696_3", product=product, interval=interval, on="normalized_read_count", T)
|
|
7
|
389
|
|
9
|
390 one = triplets[triplets$Sample == "16278_Left",]
|
|
|
391 two = triplets[triplets$Sample == "26402_Left",]
|
|
|
392 three = triplets[triplets$Sample == "26759_Left",]
|
|
8
|
393 tripletAnalysis(one, "16278_Left", two, "26402_Left", three, "26759_Left", product=product, interval=interval, on="normalized_read_count", T)
|
|
7
|
394
|
|
9
|
395 one = triplets[triplets$Sample == "16278_Right",]
|
|
|
396 two = triplets[triplets$Sample == "26402_Right",]
|
|
|
397 three = triplets[triplets$Sample == "26759_Right",]
|
|
8
|
398 tripletAnalysis(one, "16278_Right", two, "26402_Right", three, "26759_Right", product=product, interval=interval, on="normalized_read_count", T)
|
|
7
|
399
|
|
|
400
|
|
|
401 interval = intervalFreq
|
|
|
402 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
|
403 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)))
|
|
|
404
|
|
9
|
405 one = triplets[triplets$Sample == "14696_reg_BM",]
|
|
|
406 two = triplets[triplets$Sample == "24536_reg_BM",]
|
|
|
407 three = triplets[triplets$Sample == "24062_reg_BM",]
|
|
8
|
408 tripletAnalysis(one, "14696_1", two, "14696_2", three, "14696_3", product=product, interval=interval, on="Frequency", F)
|
|
7
|
409
|
|
9
|
410 one = triplets[triplets$Sample == "16278_Left",]
|
|
|
411 two = triplets[triplets$Sample == "26402_Left",]
|
|
|
412 three = triplets[triplets$Sample == "26759_Left",]
|
|
8
|
413 tripletAnalysis(one, "16278_Left", two, "26402_Left", three, "26759_Left", product=product, interval=interval, on="Frequency", F)
|
|
7
|
414
|
|
9
|
415 one = triplets[triplets$Sample == "16278_Right",]
|
|
|
416 two = triplets[triplets$Sample == "26402_Right",]
|
|
|
417 three = triplets[triplets$Sample == "26759_Right",]
|
|
8
|
418 tripletAnalysis(one, "16278_Right", two, "26402_Right", three, "26759_Right", product=product, interval=interval, on="Frequency", F)
|
