|
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
|
|
|
11 require(ggplot2)
|
|
|
12 require(reshape2)
|
|
|
13 require(data.table)
|
|
|
14 require(grid)
|
|
|
15 #require(xtable)
|
|
3
|
16 cat("<tr><td>Reading input</td></tr>", file=logfile, append=T)
|
|
2
|
17 dat = read.csv(inFile, sep="\t")
|
|
|
18 #dat = data.frame(fread(inFile)) #faster but with a dep
|
|
0
|
19 setwd(outDir)
|
|
3
|
20 cat("<tr><td>Selecting first V/J Genes</td></tr>", file=logfile, append=T)
|
|
2
|
21 dat$V_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$V_Segment_Major_Gene), "; "), "[[", 1)))
|
|
|
22 dat$J_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$J_Segment_Major_Gene), "; "), "[[", 1)))
|
|
|
23
|
|
3
|
24 cat("<tr><td>Calculating Frequency</td></tr>", file=logfile, append=T)
|
|
0
|
25 dat$Frequency = ((10^dat$Log10_Frequency)*100)
|
|
2
|
26
|
|
3
|
27 dat = dat[dat$Frequency >= min_freq,]
|
|
|
28
|
|
|
29 cat("<tr><td>Normalizing cell count to 1.000.000</td></tr>", file=logfile, append=T)
|
|
2
|
30 dat$normalized_read_count = round(dat$Clone_Molecule_Count_From_Spikes / dat$Cell_Count * 1000000 / 2)
|
|
3
|
31 dat = dat[dat$normalized_read_count >= min_cells,]
|
|
2
|
32 dat$paste = paste(dat$Sample, dat$V_Segment_Major_Gene, dat$J_Segment_Major_Gene, dat$CDR3_Sense_Sequence)
|
|
3
|
33 cat("<tr><td>Removing duplicates</td></tr>", file=logfile, append=T)
|
|
2
|
34 dat = dat[!duplicated(dat$paste),]
|
|
0
|
35 patients = split(dat, dat$Patient, drop=T)
|
|
|
36 intervalReads = rev(c(0,2,10,100,1000,10000))
|
|
|
37 intervalFreq = rev(c(0,0.01,0.1,0.5,1,5))
|
|
|
38 V_Segments = c(".*", "IGHV", "IGHD", "IGKV", "IGKV", "IgKINTR", "TRGV", "TRDV", "TRDD" , "TRBV")
|
|
|
39 J_Segments = c(".*", ".*", ".*", "IGKJ", "KDE", ".*", ".*", ".*", ".*", ".*")
|
|
|
40 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")
|
|
|
41 Titles = factor(Titles, levels=Titles)
|
|
|
42 TitlesOrder = data.frame("Title"=Titles, "TitlesOrder"=1:length(Titles))
|
|
|
43
|
|
|
44 patientCountOnColumn <- function(x, product, interval, on, appendtxt=F){
|
|
|
45 x$Sample = factor(x$Sample, levels=unique(x$Sample))
|
|
|
46 onShort = "reads"
|
|
|
47 if(on == "Frequency"){
|
|
|
48 onShort = "freq"
|
|
|
49 }
|
|
|
50 splt = split(x, x$Sample, drop=T)
|
|
2
|
51 if(length(splt) == 1){
|
|
3
|
52 print(paste(paste(x[1,which(colnames(x) == "Patient")]), "has one sample"))
|
|
|
53 splt[[2]] = data.frame("Patient" = 'NA', "Receptor" = 'NA', "Sample" = 'NA', "Cell_Count" = 100, "Clone_Molecule_Count_From_Spikes" = 10, "Log10_Frequency" = 1, "Total_Read_Count" = 100, "dsMol_per_1e6_cells" = 100, "J_Segment_Major_Gene" = 'NA', "V_Segment_Major_Gene" = 'NA', "Clone_Sequence" = 'NA', "CDR3_Sense_Sequence" = 'NA', "Related_to_leukemia_clone" = FALSE, "Frequency"= 0, "normalized_read_count" = 0, "paste" = 'a')
|
|
2
|
54 }
|
|
0
|
55 patient1 = splt[[1]]
|
|
|
56 patient2 = splt[[2]]
|
|
|
57
|
|
|
58 threshholdIndex = which(colnames(product) == "interval")
|
|
|
59 V_SegmentIndex = which(colnames(product) == "V_Segments")
|
|
|
60 J_SegmentIndex = which(colnames(product) == "J_Segments")
|
|
|
61 titleIndex = which(colnames(product) == "Titles")
|
|
|
62 sampleIndex = which(colnames(x) == "Sample")
|
|
|
63 patientIndex = which(colnames(x) == "Patient")
|
|
|
64 oneSample = paste(patient1[1,sampleIndex], sep="")
|
|
|
65 twoSample = paste(patient2[1,sampleIndex], sep="")
|
|
|
66 patient = paste(x[1,patientIndex])
|
|
3
|
67
|
|
0
|
68 switched = F
|
|
|
69 if(length(grep(".*_Right$", twoSample)) == 1 || length(grep(".*_Dx_BM$", twoSample)) == 1 || length(grep(".*_Dx$", twoSample)) == 1 ){
|
|
|
70 tmp = twoSample
|
|
|
71 twoSample = oneSample
|
|
|
72 oneSample = tmp
|
|
|
73 tmp = patient1
|
|
|
74 patient1 = patient2
|
|
|
75 patient2 = tmp
|
|
|
76 switched = T
|
|
|
77 }
|
|
|
78 if(appendtxt){
|
|
|
79 cat(paste(patient, oneSample, twoSample, sep="\t"), file="patients.txt", append=T, sep="", fill=3)
|
|
|
80 }
|
|
3
|
81 cat(paste("<tr><td>", patient, "</td></tr>", sep=""), file=logfile, append=T)
|
|
0
|
82 patientMerge = merge(patient1, patient2, by="Clone_Sequence")
|
|
|
83 res1 = vector()
|
|
|
84 res2 = vector()
|
|
|
85 resBoth = vector()
|
|
|
86 read1Count = vector()
|
|
|
87 read2Count = vector()
|
|
2
|
88 locussum1 = vector()
|
|
|
89 locussum2 = vector()
|
|
0
|
90
|
|
2
|
91 print(patient)
|
|
0
|
92 #for(iter in 1){
|
|
|
93 for(iter in 1:length(product[,1])){
|
|
|
94 threshhold = product[iter,threshholdIndex]
|
|
|
95 V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="")
|
|
|
96 J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="")
|
|
|
97 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)
|
|
|
98 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,]$Clone_Sequence))
|
|
|
99 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,]$Clone_Sequence))
|
|
2
|
100 read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count) + sum(patientMerge[both,]$normalized_read_count.x))
|
|
|
101 read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count) + sum(patientMerge[both,]$normalized_read_count.y))
|
|
0
|
102 res1 = append(res1, sum(one))
|
|
2
|
103 res2 = append(res2, sum(two))
|
|
0
|
104 resBoth = append(resBoth, sum(both))
|
|
2
|
105 locussum1 = append(locussum1, sum(patient1[(grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene)),]$normalized_read_count))
|
|
|
106 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
|
107 #threshhold = 0
|
|
|
108 if(threshhold != 0){
|
|
|
109 if(sum(one) > 0){
|
|
2
|
110 dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence")]
|
|
|
111 colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Sequence")
|
|
0
|
112 filenameOne = paste(oneSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
113 write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
114 }
|
|
|
115 if(sum(two) > 0){
|
|
2
|
116 dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence")]
|
|
|
117 colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Sequence")
|
|
0
|
118 filenameTwo = paste(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
119 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
120 }
|
|
|
121 }
|
|
|
122 if(sum(both) > 0){
|
|
2
|
123 dfBoth = patientMerge[both,c("V_Segment_Major_Gene.x", "J_Segment_Major_Gene.x", "normalized_read_count.x", "Frequency.x", "Clone_Sequence", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "normalized_read_count.y", "Frequency.y")]
|
|
|
124 colnames(dfBoth) = c(paste("Proximal segment", oneSample), paste("Distal segment", oneSample), paste("Normalized_Read_Count", oneSample), paste("Frequency", oneSample), "Sequence", paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Normalized_Read_Count", twoSample), paste("Frequency", twoSample))
|
|
0
|
125 filenameBoth = paste(oneSample, "_", twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
|
|
|
126 write.table(dfBoth, file=paste(filenameBoth, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
127 }
|
|
|
128 }
|
|
2
|
129 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
|
130 if(sum(is.na(patientResult$percentage)) > 0){
|
|
|
131 patientResult[is.na(patientResult$percentage),]$percentage = 0
|
|
|
132 }
|
|
|
133 colnames(patientResult)[6] = oneSample
|
|
|
134 colnames(patientResult)[8] = twoSample
|
|
|
135 colnamesBak = colnames(patientResult)
|
|
2
|
136 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
|
137 write.table(patientResult, file=paste(patient, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
|
|
|
138 colnames(patientResult) = colnamesBak
|
|
|
139
|
|
|
140 patientResult$Locus = factor(patientResult$Locus, Titles)
|
|
|
141 patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))
|
|
|
142
|
|
|
143 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "Both")])
|
|
|
144 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=Both), stat='identity', position="dodge", fill="#79c36a")
|
|
|
145 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
146 plt = plt + geom_text(aes(ymax=max(Both), x=cut_off_value,y=Both,label=Both), angle=90, hjust=0)
|
|
|
147 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in both")
|
|
|
148 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
|
149 png(paste(patient, "_", onShort, ".png", sep=""), width=1920, height=1080)
|
|
|
150 print(plt)
|
|
|
151 dev.off()
|
|
|
152 #(t,r,b,l)
|
|
|
153 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "percentage")])
|
|
|
154 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=percentage), stat='identity', position="dodge", fill="#79c36a")
|
|
|
155 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
156 plt = plt + geom_text(aes(ymax=max(percentage), x=cut_off_value,y=percentage,label=percentage), angle=90, hjust=0)
|
|
|
157 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("% clones in both left and right")
|
|
|
158 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
|
|
|
159 png(paste(patient, "_percent_", onShort, ".png", sep=""), width=1920, height=1080)
|
|
|
160 print(plt)
|
|
|
161 dev.off()
|
|
|
162
|
|
|
163 patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample)] ,id.vars=1:2)
|
|
|
164 patientResult$relativeValue = patientResult$value * 10
|
|
|
165 patientResult[patientResult$relativeValue == 0,]$relativeValue = 1
|
|
|
166 plt = ggplot(patientResult)
|
|
|
167 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge")
|
|
|
168 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
|
|
|
169 plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9)))
|
|
|
170 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)
|
|
|
171 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)
|
|
|
172 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle(paste("Number of clones in only ", oneSample, " and only ", twoSample, sep=""))
|
|
|
173 png(paste(patient, "_", onShort, "_both.png", sep=""), width=1920, height=1080)
|
|
|
174 print(plt)
|
|
|
175 dev.off()
|
|
|
176 }
|
|
|
177
|
|
3
|
178 cat("<tr><td>Starting Frequency analysis</td></tr>", file=logfile, append=T)
|
|
|
179
|
|
0
|
180 interval = intervalFreq
|
|
|
181 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
|
182 product = data.frame("Titles"=rep(Titles, each=6), "interval"=rep(interval, times=10), "V_Segments"=rep(V_Segments, each=6), "J_Segments"=rep(J_Segments, each=6))
|
|
|
183 #patientFrequencyCount(patient1)
|
|
|
184 #lapply(patients[c(5,6,10)], FUN=patientFrequencyCount)
|
|
|
185 #lapply(patients[c(5,6,7,8,13)], FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
|
|
|
186 #lapply(patients[c(6,7,8)], FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
|
|
|
187 #lapply(patients[c(6)], FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
|
|
|
188 lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
|
|
|
189
|
|
3
|
190 cat("<tr><td>Starting Cell Count analysis</td></tr>", file=logfile, append=T)
|
|
|
191
|
|
0
|
192 interval = intervalReads
|
|
|
193 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
|
|
|
194 product = data.frame("Titles"=rep(Titles, each=6), "interval"=rep(interval, times=10), "V_Segments"=rep(V_Segments, each=6), "J_Segments"=rep(J_Segments, each=6))
|
|
|
195 #patientResult = patientReadCount(patient1)
|
|
|
196 #lapply(patients[c(5,6,10)], FUN=patientReadCount)
|
|
|
197 #lapply(patients[c(5,6,7,8,13)], FUN=patientCountOnColumn, product = product, interval=interval, on="Clone_Molecule_Count_From_Spikes")
|
|
|
198 #lapply(patients[c(6)], FUN=patientCountOnColumn, product = product, interval=interval, on="Clone_Molecule_Count_From_Spikes")
|
|
|
199 lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Clone_Molecule_Count_From_Spikes")
|
|
|
200
|
|
3
|
201 cat("</table></html>", file=logfile, append=T)
|
|
|
202
|
