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)
|
13
|
18 dat = read.table(inFile, header=T, sep="\t", dec=".", fill=T, stringsAsFactors=F)
|
9
|
19 dat = dat[!is.na(dat$Patient),]
|
13
|
20 dat$Related_to_leukemia_clone = F
|
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
|
3
|
27 cat("<tr><td>Calculating Frequency</td></tr>", file=logfile, append=T)
|
12
|
28
|
13
|
29 dat$Frequency = ((10^dat$Log10_Frequency)*100)
|
2
|
30
|
3
|
31 dat = dat[dat$Frequency >= min_freq,]
|
|
32
|
13
|
33 triplets = dat[grepl("VanDongen_cALL_14696", dat$Patient) | grepl("(16278)|(26402)|(26759)", dat$Sample),]
|
|
34
|
|
35 cat("<tr><td>Normalizing to lowest cell count within locus</td></tr>", file=logfile, append=T)
|
|
36
|
|
37 dat$locus_V = substring(dat$V_Segment_Major_Gene, 0, 4)
|
|
38 dat$locus_J = substring(dat$J_Segment_Major_Gene, 0, 4)
|
|
39 min_cell_count = data.frame(data.table(dat)[, list(min_cell_count=min(.SD$Cell_Count)), by=c("Patient", "locus_V", "locus_J")])
|
|
40
|
|
41 dat$min_cell_paste = paste(dat$Patient, dat$locus_V, dat$locus_J)
|
|
42 min_cell_count$min_cell_paste = paste(min_cell_count$Patient, min_cell_count$locus_V, min_cell_count$locus_J)
|
|
43
|
|
44 min_cell_count = min_cell_count[,c("min_cell_paste", "min_cell_count")]
|
|
45
|
|
46 dat = merge(dat, min_cell_count, by="min_cell_paste")
|
|
47
|
|
48 dat$normalized_read_count = round(dat$Clone_Molecule_Count_From_Spikes / dat$Cell_Count * dat$min_cell_count / 2, digits=2) #??????????????????????????????????? wel of geen / 2
|
|
49
|
3
|
50 dat = dat[dat$normalized_read_count >= min_cells,]
|
13
|
51
|
|
52 dat$paste = paste(dat$Sample, dat$Clone_Sequence)
|
9
|
53
|
22
|
54 patients = split(dat, dat$Patient, drop=T)
|
9
|
55 intervalReads = rev(c(0,10,25,50,100,250,500,750,1000,10000))
|
6
|
56 intervalFreq = rev(c(0,0.01,0.05,0.1,0.5,1,5))
|
0
|
57 V_Segments = c(".*", "IGHV", "IGHD", "IGKV", "IGKV", "IgKINTR", "TRGV", "TRDV", "TRDD" , "TRBV")
|
|
58 J_Segments = c(".*", ".*", ".*", "IGKJ", "KDE", ".*", ".*", ".*", ".*", ".*")
|
|
59 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")
|
|
60 Titles = factor(Titles, levels=Titles)
|
|
61 TitlesOrder = data.frame("Title"=Titles, "TitlesOrder"=1:length(Titles))
|
|
62
|
|
63 patientCountOnColumn <- function(x, product, interval, on, appendtxt=F){
|
|
64 x$Sample = factor(x$Sample, levels=unique(x$Sample))
|
|
65 onShort = "reads"
|
|
66 if(on == "Frequency"){
|
|
67 onShort = "freq"
|
|
68 }
|
18
|
69 onx = paste(on, ".x", sep="")
|
|
70 ony = paste(on, ".y", sep="")
|
0
|
71 splt = split(x, x$Sample, drop=T)
|
4
|
72 type="pair"
|
2
|
73 if(length(splt) == 1){
|
3
|
74 print(paste(paste(x[1,which(colnames(x) == "Patient")]), "has one sample"))
|
4
|
75 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))
|
|
76 type="single"
|
2
|
77 }
|
0
|
78 patient1 = splt[[1]]
|
|
79 patient2 = splt[[2]]
|
|
80
|
|
81 threshholdIndex = which(colnames(product) == "interval")
|
|
82 V_SegmentIndex = which(colnames(product) == "V_Segments")
|
|
83 J_SegmentIndex = which(colnames(product) == "J_Segments")
|
|
84 titleIndex = which(colnames(product) == "Titles")
|
|
85 sampleIndex = which(colnames(x) == "Sample")
|
|
86 patientIndex = which(colnames(x) == "Patient")
|
|
87 oneSample = paste(patient1[1,sampleIndex], sep="")
|
|
88 twoSample = paste(patient2[1,sampleIndex], sep="")
|
|
89 patient = paste(x[1,patientIndex])
|
3
|
90
|
0
|
91 switched = F
|
|
92 if(length(grep(".*_Right$", twoSample)) == 1 || length(grep(".*_Dx_BM$", twoSample)) == 1 || length(grep(".*_Dx$", twoSample)) == 1 ){
|
|
93 tmp = twoSample
|
|
94 twoSample = oneSample
|
|
95 oneSample = tmp
|
|
96 tmp = patient1
|
|
97 patient1 = patient2
|
|
98 patient2 = tmp
|
|
99 switched = T
|
|
100 }
|
|
101 if(appendtxt){
|
4
|
102 cat(paste(patient, oneSample, twoSample, type, sep="\t"), file="patients.txt", append=T, sep="", fill=3)
|
0
|
103 }
|
3
|
104 cat(paste("<tr><td>", patient, "</td></tr>", sep=""), file=logfile, append=T)
|
9
|
105
|
12
|
106 #patient1$merge = paste(patient1$V_Segment_Major_Gene, patient1$J_Segment_Major_Gene, patient1$CDR3_Sense_Sequence)
|
|
107 #patient2$merge = paste(patient2$V_Segment_Major_Gene, patient2$J_Segment_Major_Gene, patient2$CDR3_Sense_Sequence)
|
13
|
108 patient1$merge = paste(patient1$Clone_Sequence)
|
|
109 patient2$merge = paste(patient2$Clone_Sequence)
|
9
|
110
|
12
|
111 #patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge")
|
9
|
112 patientMerge = merge(patient1, patient2, by.x="merge", by.y="merge")
|
18
|
113 patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony])
|
0
|
114 res1 = vector()
|
|
115 res2 = vector()
|
|
116 resBoth = vector()
|
|
117 read1Count = vector()
|
|
118 read2Count = vector()
|
2
|
119 locussum1 = vector()
|
|
120 locussum2 = vector()
|
9
|
121
|
|
122 print(patient)
|
0
|
123 #for(iter in 1){
|
|
124 for(iter in 1:length(product[,1])){
|
|
125 threshhold = product[iter,threshholdIndex]
|
|
126 V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="")
|
|
127 J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="")
|
18
|
128 #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
|
|
129 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 higher than threshold
|
19
|
130 one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$Clone_Sequence %in% patientMerge[both,]$merge))
|
|
131 two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$Clone_Sequence %in% patientMerge[both,]$merge))
|
14
|
132 read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count))
|
|
133 read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count))
|
0
|
134 res1 = append(res1, sum(one))
|
2
|
135 res2 = append(res2, sum(two))
|
0
|
136 resBoth = append(resBoth, sum(both))
|
2
|
137 locussum1 = append(locussum1, sum(patient1[(grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene)),]$normalized_read_count))
|
|
138 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
|
139 #threshhold = 0
|
|
140 if(threshhold != 0){
|
|
141 if(sum(one) > 0){
|
15
|
142 dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
|
|
143 colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone Sequence", "Related_to_leukemia_clone")
|
0
|
144 filenameOne = paste(oneSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
145 write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
146 }
|
|
147 if(sum(two) > 0){
|
15
|
148 dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
|
|
149 colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone Sequence", "Related_to_leukemia_clone")
|
0
|
150 filenameTwo = paste(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
151 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
152 }
|
|
153 }
|
|
154 if(sum(both) > 0){
|
15
|
155 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")]
|
|
156 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))
|
0
|
157 filenameBoth = paste(oneSample, "_", twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
158 write.table(dfBoth, file=paste(filenameBoth, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
159 }
|
|
160 }
|
2
|
161 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
|
162 if(sum(is.na(patientResult$percentage)) > 0){
|
|
163 patientResult[is.na(patientResult$percentage),]$percentage = 0
|
|
164 }
|
|
165 colnames(patientResult)[6] = oneSample
|
|
166 colnames(patientResult)[8] = twoSample
|
|
167 colnamesBak = colnames(patientResult)
|
2
|
168 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
|
169 write.table(patientResult, file=paste(patient, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
170 colnames(patientResult) = colnamesBak
|
|
171
|
|
172 patientResult$Locus = factor(patientResult$Locus, Titles)
|
|
173 patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))
|
|
174
|
|
175 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "Both")])
|
|
176 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=Both), stat='identity', position="dodge", fill="#79c36a")
|
|
177 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
178 plt = plt + geom_text(aes(ymax=max(Both), x=cut_off_value,y=Both,label=Both), angle=90, hjust=0)
|
|
179 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in both")
|
|
180 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
181 png(paste(patient, "_", onShort, ".png", sep=""), width=1920, height=1080)
|
|
182 print(plt)
|
|
183 dev.off()
|
|
184 #(t,r,b,l)
|
|
185 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "percentage")])
|
|
186 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=percentage), stat='identity', position="dodge", fill="#79c36a")
|
|
187 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
188 plt = plt + geom_text(aes(ymax=max(percentage), x=cut_off_value,y=percentage,label=percentage), angle=90, hjust=0)
|
|
189 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("% clones in both left and right")
|
|
190 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
191 png(paste(patient, "_percent_", onShort, ".png", sep=""), width=1920, height=1080)
|
|
192 print(plt)
|
|
193 dev.off()
|
|
194
|
|
195 patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample)] ,id.vars=1:2)
|
|
196 patientResult$relativeValue = patientResult$value * 10
|
|
197 patientResult[patientResult$relativeValue == 0,]$relativeValue = 1
|
|
198 plt = ggplot(patientResult)
|
|
199 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge")
|
|
200 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
201 plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9)))
|
|
202 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)
|
|
203 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)
|
|
204 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle(paste("Number of clones in only ", oneSample, " and only ", twoSample, sep=""))
|
|
205 png(paste(patient, "_", onShort, "_both.png", sep=""), width=1920, height=1080)
|
|
206 print(plt)
|
|
207 dev.off()
|
|
208 }
|
|
209
|
3
|
210 cat("<tr><td>Starting Frequency analysis</td></tr>", file=logfile, append=T)
|
|
211
|
0
|
212 interval = intervalFreq
|
|
213 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
4
|
214 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)))
|
|
215 mclapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
|
0
|
216
|
3
|
217 cat("<tr><td>Starting Cell Count analysis</td></tr>", file=logfile, append=T)
|
|
218
|
0
|
219 interval = intervalReads
|
|
220 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
4
|
221 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
|
222 mclapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="normalized_read_count")
|
0
|
223
|
3
|
224 cat("</table></html>", file=logfile, append=T)
|
|
225
|
7
|
226
|
13
|
227
|
7
|
228 tripletAnalysis <- function(patient1, label1, patient2, label2, patient3, label3, product, interval, on, appendTriplets= FALSE){
|
|
229 onShort = "reads"
|
|
230 if(on == "Frequency"){
|
|
231 onShort = "freq"
|
|
232 }
|
18
|
233 onx = paste(on, ".x", sep="")
|
|
234 ony = paste(on, ".y", sep="")
|
|
235 onz = paste(on, ".z", sep="")
|
7
|
236 type="triplet"
|
|
237
|
|
238 threshholdIndex = which(colnames(product) == "interval")
|
|
239 V_SegmentIndex = which(colnames(product) == "V_Segments")
|
|
240 J_SegmentIndex = which(colnames(product) == "J_Segments")
|
|
241 titleIndex = which(colnames(product) == "Titles")
|
|
242 sampleIndex = which(colnames(patient1) == "Sample")
|
|
243 patientIndex = which(colnames(patient1) == "Patient")
|
|
244 oneSample = paste(patient1[1,sampleIndex], sep="")
|
|
245 twoSample = paste(patient2[1,sampleIndex], sep="")
|
|
246 threeSample = paste(patient3[1,sampleIndex], sep="")
|
|
247
|
12
|
248 #patient1$merge = paste(patient1$V_Segment_Major_Gene, patient1$J_Segment_Major_Gene, patient1$CDR3_Sense_Sequence)
|
|
249 #patient2$merge = paste(patient2$V_Segment_Major_Gene, patient2$J_Segment_Major_Gene, patient2$CDR3_Sense_Sequence)
|
|
250 #patient3$merge = paste(patient3$V_Segment_Major_Gene, patient3$J_Segment_Major_Gene, patient3$CDR3_Sense_Sequence)
|
|
251
|
15
|
252 patient1$merge = paste(patient1$Clone_Sequence)
|
|
253 patient2$merge = paste(patient2$Clone_Sequence)
|
|
254 patient3$merge = paste(patient3$Clone_Sequence)
|
9
|
255
|
|
256 patientMerge = merge(patient1, patient2, by="merge")
|
|
257 patientMerge = merge(patientMerge, patient3, by="merge")
|
13
|
258 colnames(patientMerge)[30:length(colnames(patientMerge))] = paste(colnames(patientMerge)[30:length(colnames(patientMerge))], ".z", sep="")
|
18
|
259 patientMerge$thresholdValue = pmax(patientMerge[,onx], patientMerge[,ony], patientMerge[,onz])
|
20
|
260
|
|
261 patientMerge12 = merge(patient1, patient2, by="merge")
|
|
262 patientMerge12$thresholdValue = pmax(patientMerge12[,onx], patientMerge12[,ony])
|
|
263 patientMerge13 = merge(patient1, patient3, by="merge")
|
|
264 patientMerge13$thresholdValue = pmax(patientMerge13[,onx], patientMerge13[,ony])
|
|
265 patientMerge23 = merge(patient2, patient3, by="merge")
|
|
266 patientMerge23$thresholdValue = pmax(patientMerge23[,onx], patientMerge23[,ony])
|
|
267
|
24
|
268 scatterplot_data_columns = c("Clone_Sequence", "Frequency", "normalized_read_count", "V_Segment_Major_Gene", "J_Segment_Major_Gene")
|
20
|
269 scatterplot_data = rbind(patient1[,scatterplot_data_columns], patient2[,scatterplot_data_columns], patient3[,scatterplot_data_columns])
|
|
270 scatterplot_data = scatterplot_data[!duplicated(scatterplot_data$Clone_Sequence),]
|
|
271 scatterplot_data$type = factor(x="single", levels=c("In one", "In two", "In three"))
|
|
272
|
7
|
273 res1 = vector()
|
|
274 res2 = vector()
|
|
275 res3 = vector()
|
20
|
276 res12 = vector()
|
|
277 res13 = vector()
|
|
278 res23 = vector()
|
7
|
279 resAll = vector()
|
|
280 read1Count = vector()
|
|
281 read2Count = vector()
|
|
282 read3Count = vector()
|
|
283
|
|
284 if(appendTriplets){
|
9
|
285 cat(paste(label1, label2, label3, sep="\t"), file="triplets.txt", append=T, sep="", fill=3)
|
7
|
286 }
|
|
287 for(iter in 1:length(product[,1])){
|
|
288 threshhold = product[iter,threshholdIndex]
|
|
289 V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="")
|
|
290 J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="")
|
18
|
291 #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)
|
|
292 all = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge$thresholdValue > threshhold)
|
7
|
293
|
25
|
294 one_two = (grepl(V_Segment, patientMerge12$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge12$J_Segment_Major_Gene.x) & patientMerge12$thresholdValue > threshhold & !(patientMerge12$Clone_Sequence.x %in% patientMerge[all,]$merge))
|
|
295 one_three = (grepl(V_Segment, patientMerge13$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge13$J_Segment_Major_Gene.x) & patientMerge13$thresholdValue > threshhold & !(patientMerge13$Clone_Sequence.x %in% patientMerge[all,]$merge))
|
|
296 two_three = (grepl(V_Segment, patientMerge23$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge23$J_Segment_Major_Gene.x) & patientMerge23$thresholdValue > threshhold & !(patientMerge23$Clone_Sequence.x %in% patientMerge[all,]$merge))
|
24
|
297
|
25
|
298 one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$Clone_Sequence %in% patientMerge[all,]$merge) & !(patient1$Clone_Sequence %in% patientMerge12[one_two,]$merge) & !(patient1$Clone_Sequence %in% patientMerge13[one_three,]$merge))
|
|
299 two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$Clone_Sequence %in% patientMerge[all,]$merge) & !(patient2$Clone_Sequence %in% patientMerge12[one_two,]$merge) & !(patient2$Clone_Sequence %in% patientMerge23[two_three,]$merge))
|
|
300 three = (grepl(V_Segment, patient3$V_Segment_Major_Gene) & grepl(J_Segment, patient3$J_Segment_Major_Gene) & patient3[,on] > threshhold & !(patient3$Clone_Sequence %in% patientMerge[all,]$merge) & !(patient3$Clone_Sequence %in% patientMerge13[one_three,]$merge) & !(patient3$Clone_Sequence %in% patientMerge23[two_three,]$merge))
|
20
|
301
|
18
|
302 read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.x))
|
|
303 read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.y))
|
|
304 read3Count = append(read3Count, sum(patient3[three,]$normalized_read_count) + sum(patientMerge[all,]$normalized_read_count.z))
|
7
|
305 res1 = append(res1, sum(one))
|
|
306 res2 = append(res2, sum(two))
|
|
307 res3 = append(res3, sum(three))
|
|
308 resAll = append(resAll, sum(all))
|
20
|
309 res12 = append(res12, sum(one_two))
|
|
310 res13 = append(res13, sum(one_three))
|
|
311 res23 = append(res23, sum(two_three))
|
7
|
312 #threshhold = 0
|
|
313 if(threshhold != 0){
|
|
314 if(sum(one) > 0){
|
20
|
315 dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
|
|
316 colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")
|
7
|
317 filenameOne = paste(label1, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
318 write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
319 }
|
|
320 if(sum(two) > 0){
|
20
|
321 dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
|
|
322 colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")
|
7
|
323 filenameTwo = paste(label2, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
324 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
325 }
|
|
326 if(sum(three) > 0){
|
20
|
327 dfThree = patient3[three,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")]
|
|
328 colnames(dfThree) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Clone_Sequence", "Related_to_leukemia_clone")
|
7
|
329 filenameThree = paste(label3, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
330 write.table(dfThree, file=paste(filenameThree, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
331 }
|
20
|
332 if(sum(one_two) > 0){
|
|
333 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")]
|
|
334 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))
|
|
335 filenameOne_two = paste(label1, "_", label2, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="")
|
|
336 write.table(dfOne_two, file=paste(filenameOne_two, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
337 }
|
|
338 if(sum(one_three) > 0){
|
|
339 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")]
|
|
340 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))
|
|
341 filenameOne_three = paste(label1, "_", label3, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="")
|
|
342 write.table(dfOne_three, file=paste(filenameOne_three, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
343 }
|
|
344 if(sum(two_three) > 0){
|
|
345 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")]
|
|
346 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))
|
|
347 filenameTwo_three = paste(label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, onShort, sep="")
|
|
348 write.table(dfTwo_three, file=paste(filenameTwo_three, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
349 }
|
|
350 } else { #scatterplot data
|
24
|
351 scatterplot_locus_data = scatterplot_data[grepl(V_Segment, scatterplot_data$V_Segment_Major_Gene) & grepl(J_Segment, scatterplot_data$J_Segment_Major_Gene),]
|
20
|
352 scatterplot_locus_data$type = ifelse(scatterplot_locus_data$Clone_Sequence %in% patientMerge12[one_two,]$Clone_Sequence.x, "In two", "In one")
|
|
353 scatterplot_locus_data$type = ifelse(scatterplot_locus_data$Clone_Sequence %in% patientMerge13[one_three,]$Clone_Sequence.x, "In two", "In one")
|
|
354 scatterplot_locus_data$type = ifelse(scatterplot_locus_data$Clone_Sequence %in% patientMerge23[two_three,]$Clone_Sequence.x, "In two", "In one")
|
|
355 scatterplot_locus_data$type = ifelse(scatterplot_locus_data$type == "In two", ifelse(scatterplot_locus_data$Clone_Sequence %in% patientMerge[all,]$Clone_Sequence.x, "In three", "In two"), "In one")
|
|
356 scatterplot_locus_data$type = ifelse(scatterplot_locus_data$type == "In one", "In one", "In multiple")
|
|
357
|
|
358 p = NULL
|
|
359 if(nrow(scatterplot_locus_data) != 0){
|
|
360 if(on == "normalized_read_count"){
|
|
361 scales = 10^(0:ceiling(log10(max(scatterplot_locus_data$normalized_read_count))))
|
|
362 p = ggplot(scatterplot_locus_data, aes(type, normalized_read_count)) + scale_y_log10(breaks=scales,labels=scales)
|
|
363 } else {
|
|
364 p = ggplot(scatterplot_locus_data, aes(type, Frequency))
|
|
365 }
|
|
366 p = p + geom_point(aes(colour=type), position="jitter")
|
25
|
367 p = p + xlab("In one or in multiple samples") + ylab(onShort) + ggtitle(paste(label1, label2, label3, onShort, product[iter, titleIndex]))
|
20
|
368 } else {
|
|
369 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]))
|
|
370 }
|
|
371 png(paste(label1, "_", label2, "_", label3, "_", onShort, "_", product[iter, titleIndex],"_scatter.png", sep=""))
|
|
372 print(p)
|
|
373 dev.off()
|
|
374 }
|
7
|
375 if(sum(all) > 0){
|
20
|
376 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")]
|
|
377 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))
|
7
|
378 filenameAll = paste(label1, "_", label2, "_", label3, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
379 write.table(dfAll, file=paste(filenameAll, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
380 }
|
|
381 }
|
20
|
382 #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))
|
|
383 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)
|
7
|
384 colnames(patientResult)[6] = oneSample
|
20
|
385 colnames(patientResult)[7] = twoSample
|
|
386 colnames(patientResult)[8] = threeSample
|
|
387 colnames(patientResult)[9] = paste(oneSample, twoSample, sep="_")
|
|
388 colnames(patientResult)[10] = paste(oneSample, twoSample, sep="_")
|
|
389 colnames(patientResult)[11] = paste(oneSample, twoSample, sep="_")
|
7
|
390
|
|
391 colnamesBak = colnames(patientResult)
|
20
|
392 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))
|
7
|
393 write.table(patientResult, file=paste(label1, "_", label2, "_", label3, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
394 colnames(patientResult) = colnamesBak
|
|
395
|
|
396 patientResult$Locus = factor(patientResult$Locus, Titles)
|
|
397 patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))
|
|
398
|
|
399 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "All")])
|
|
400 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=All), stat='identity', position="dodge", fill="#79c36a")
|
|
401 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
402 plt = plt + geom_text(aes(ymax=max(All), x=cut_off_value,y=All,label=All), angle=90, hjust=0)
|
|
403 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in All")
|
|
404 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
405 png(paste(label1, "_", label2, "_", label3, "_", onShort, "_total_all.png", sep=""), width=1920, height=1080)
|
|
406 print(plt)
|
|
407 dev.off()
|
|
408
|
|
409 fontSize = 4
|
|
410
|
|
411 bak = patientResult
|
|
412 patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample, threeSample)] ,id.vars=1:2)
|
|
413 patientResult$relativeValue = patientResult$value * 10
|
|
414 patientResult[patientResult$relativeValue == 0,]$relativeValue = 1
|
|
415 plt = ggplot(patientResult)
|
|
416 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge")
|
|
417 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
418 plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9)))
|
|
419 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)
|
|
420 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)
|
|
421 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)
|
|
422 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in only one sample")
|
|
423 png(paste(label1, "_", label2, "_", label3, "_", onShort, "_indiv_all.png", sep=""), width=1920, height=1080)
|
|
424 print(plt)
|
|
425 dev.off()
|
|
426 }
|
|
427
|
9
|
428 triplets$uniqueID = "ID"
|
|
429
|
|
430 triplets[grepl("16278_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
|
|
431 triplets[grepl("26402_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
|
|
432 triplets[grepl("26759_Left", triplets$Sample),]$uniqueID = "16278_26402_26759_Left"
|
|
433
|
|
434 triplets[grepl("16278_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
|
|
435 triplets[grepl("26402_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
|
|
436 triplets[grepl("26759_Right", triplets$Sample),]$uniqueID = "16278_26402_26759_Right"
|
|
437
|
|
438 triplets[grepl("14696", triplets$Patient),]$uniqueID = "14696"
|
|
439
|
13
|
440 triplets$locus_V = substring(triplets$V_Segment_Major_Gene, 0, 4)
|
|
441 triplets$locus_J = substring(triplets$J_Segment_Major_Gene, 0, 4)
|
|
442 min_cell_count = data.frame(data.table(triplets)[, list(min_cell_count=min(.SD$Cell_Count)), by=c("uniqueID", "locus_V", "locus_J")])
|
|
443
|
|
444 triplets$min_cell_paste = paste(triplets$uniqueID, triplets$locus_V, triplets$locus_J)
|
|
445 min_cell_count$min_cell_paste = paste(min_cell_count$uniqueID, min_cell_count$locus_V, min_cell_count$locus_J)
|
|
446
|
|
447 min_cell_count = min_cell_count[,c("min_cell_paste", "min_cell_count")]
|
9
|
448
|
13
|
449 triplets = merge(triplets, min_cell_count, by="min_cell_paste")
|
|
450
|
|
451 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
|
|
452
|
|
453 triplets = triplets[triplets$normalized_read_count >= min_cells,]
|
|
454
|
|
455 column_drops = c("locus_V", "locus_J", "min_cell_count", "min_cell_paste")
|
|
456
|
|
457 triplets = triplets[,!(colnames(triplets) %in% column_drops)]
|
9
|
458
|
7
|
459 interval = intervalReads
|
|
460 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
461 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)))
|
|
462
|
9
|
463 one = triplets[triplets$Sample == "14696_reg_BM",]
|
|
464 two = triplets[triplets$Sample == "24536_reg_BM",]
|
|
465 three = triplets[triplets$Sample == "24062_reg_BM",]
|
8
|
466 tripletAnalysis(one, "14696_1", two, "14696_2", three, "14696_3", product=product, interval=interval, on="normalized_read_count", T)
|
7
|
467
|
9
|
468 one = triplets[triplets$Sample == "16278_Left",]
|
|
469 two = triplets[triplets$Sample == "26402_Left",]
|
|
470 three = triplets[triplets$Sample == "26759_Left",]
|
8
|
471 tripletAnalysis(one, "16278_Left", two, "26402_Left", three, "26759_Left", product=product, interval=interval, on="normalized_read_count", T)
|
7
|
472
|
9
|
473 one = triplets[triplets$Sample == "16278_Right",]
|
|
474 two = triplets[triplets$Sample == "26402_Right",]
|
|
475 three = triplets[triplets$Sample == "26759_Right",]
|
8
|
476 tripletAnalysis(one, "16278_Right", two, "26402_Right", three, "26759_Right", product=product, interval=interval, on="normalized_read_count", T)
|
7
|
477
|
|
478
|
|
479 interval = intervalFreq
|
|
480 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
481 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)))
|
|
482
|
9
|
483 one = triplets[triplets$Sample == "14696_reg_BM",]
|
|
484 two = triplets[triplets$Sample == "24536_reg_BM",]
|
|
485 three = triplets[triplets$Sample == "24062_reg_BM",]
|
8
|
486 tripletAnalysis(one, "14696_1", two, "14696_2", three, "14696_3", product=product, interval=interval, on="Frequency", F)
|
7
|
487
|
9
|
488 one = triplets[triplets$Sample == "16278_Left",]
|
|
489 two = triplets[triplets$Sample == "26402_Left",]
|
|
490 three = triplets[triplets$Sample == "26759_Left",]
|
8
|
491 tripletAnalysis(one, "16278_Left", two, "26402_Left", three, "26759_Left", product=product, interval=interval, on="Frequency", F)
|
7
|
492
|
9
|
493 one = triplets[triplets$Sample == "16278_Right",]
|
|
494 two = triplets[triplets$Sample == "26402_Right",]
|
|
495 three = triplets[triplets$Sample == "26759_Right",]
|
8
|
496 tripletAnalysis(one, "16278_Right", two, "26402_Right", three, "26759_Right", product=product, interval=interval, on="Frequency", F)
|