Mercurial > repos > davidvanzessen > argalaxy_tools

comparison report_clonality/RScript.r @ 58:a073fa12ef98 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Fri, 18 Mar 2016 08:02:22 -0400
parents
children 11ec9edfefee
comparison
equal deleted inserted replaced
57:16c7fc1c4bf8 58:a073fa12ef98
1 # ---------------------- load/install packages ----------------------
2
3 if (!("gridExtra" %in% rownames(installed.packages()))) {
4 install.packages("gridExtra", repos="http://cran.xl-mirror.nl/")
5 }
6 library(gridExtra)
7 if (!("ggplot2" %in% rownames(installed.packages()))) {
8 install.packages("ggplot2", repos="http://cran.xl-mirror.nl/")
9 }
10 library(ggplot2)
11 if (!("plyr" %in% rownames(installed.packages()))) {
12 install.packages("plyr", repos="http://cran.xl-mirror.nl/")
13 }
14 library(plyr)
15
16 if (!("data.table" %in% rownames(installed.packages()))) {
17 install.packages("data.table", repos="http://cran.xl-mirror.nl/")
18 }
19 library(data.table)
20
21 if (!("reshape2" %in% rownames(installed.packages()))) {
22 install.packages("reshape2", repos="http://cran.xl-mirror.nl/")
23 }
24 library(reshape2)
25
26 if (!("lymphclon" %in% rownames(installed.packages()))) {
27 install.packages("lymphclon", repos="http://cran.xl-mirror.nl/")
28 }
29 library(lymphclon)
30
31 # ---------------------- parameters ----------------------
32
33 args <- commandArgs(trailingOnly = TRUE)
34
35 infile = args[1] #path to input file
36 outfile = args[2] #path to output file
37 outdir = args[3] #path to output folder (html/images/data)
38 clonaltype = args[4] #clonaltype definition, or 'none' for no unique filtering
39 ct = unlist(strsplit(clonaltype, ","))
40 species = args[5] #human or mouse
41 locus = args[6] # IGH, IGK, IGL, TRB, TRA, TRG or TRD
42 filterproductive = ifelse(args[7] == "yes", T, F) #should unproductive sequences be filtered out? (yes/no)
43 clonality_method = args[8]
44
45
46 # ---------------------- Data preperation ----------------------
47
48 print("Report Clonality - Data preperation")
49
50 inputdata = read.table(infile, sep="\t", header=TRUE, fill=T, comment.char="")
51
52 setwd(outdir)
53
54 # remove weird rows
55 inputdata = inputdata[inputdata$Sample != "",]
56
57 #remove the allele from the V,D and J genes
58 inputdata$Top.V.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.V.Gene)
59 inputdata$Top.D.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.D.Gene)
60 inputdata$Top.J.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.J.Gene)
61
62 #filter uniques
63 inputdata.removed = inputdata[NULL,]
64
65 inputdata$clonaltype = 1:nrow(inputdata)
66
67 PRODF = inputdata
68 UNPROD = inputdata
69 if(filterproductive){
70 if("Functionality" %in% colnames(inputdata)) { # "Functionality" is an IMGT column
71 PRODF = inputdata[inputdata$Functionality == "productive" | inputdata$Functionality == "productive (see comment)", ]
72 UNPROD = inputdata[!(inputdata$Functionality == "productive" | inputdata$Functionality == "productive (see comment)"), ]
73 } else {
74 PRODF = inputdata[inputdata$VDJ.Frame != "In-frame with stop codon" & inputdata$VDJ.Frame != "Out-of-frame" & inputdata$CDR3.Found.How != "NOT_FOUND" , ]
75 UNPROD = inputdata[!(inputdata$VDJ.Frame != "In-frame with stop codon" & inputdata$VDJ.Frame != "Out-of-frame" & inputdata$CDR3.Found.How != "NOT_FOUND" ), ]
76 }
77 }
78
79 clonalityFrame = PRODF
80
81 #remove duplicates based on the clonaltype
82 if(clonaltype != "none"){
83 clonaltype = paste(clonaltype, ",Sample", sep="") #add sample column to clonaltype, unique within samples
84 PRODF$clonaltype = do.call(paste, c(PRODF[unlist(strsplit(clonaltype, ","))], sep = ":"))
85 PRODF = PRODF[!duplicated(PRODF$clonaltype), ]
86
87 UNPROD$clonaltype = do.call(paste, c(UNPROD[unlist(strsplit(clonaltype, ","))], sep = ":"))
88 UNPROD = UNPROD[!duplicated(UNPROD$clonaltype), ]
89
90 #again for clonalityFrame but with sample+replicate
91 clonalityFrame$clonaltype = do.call(paste, c(clonalityFrame[unlist(strsplit(clonaltype, ","))], sep = ":"))
92 clonalityFrame$clonality_clonaltype = do.call(paste, c(clonalityFrame[unlist(strsplit(paste(clonaltype, ",Replicate", sep=""), ","))], sep = ":"))
93 clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$clonality_clonaltype), ]
94 }
95
96 PRODF$freq = 1
97
98 if(any(grepl(pattern="_", x=PRODF$ID))){ #the frequency can be stored in the ID with the pattern ".*_freq_.*"
99 PRODF$freq = gsub("^[0-9]+_", "", PRODF$ID)
100 PRODF$freq = gsub("_.*", "", PRODF$freq)
101 PRODF$freq = as.numeric(PRODF$freq)
102 if(any(is.na(PRODF$freq))){ #if there was an "_" in the ID, but not the frequency, go back to frequency of 1 for every sequence
103 PRODF$freq = 1
104 }
105 }
106
107
108
109 #write the complete dataset that is left over, will be the input if 'none' for clonaltype and 'no' for filterproductive
110 write.table(PRODF, "allUnique.txt", sep=",",quote=F,row.names=F,col.names=T)
111 write.table(PRODF, "allUnique.csv", sep="\t",quote=F,row.names=F,col.names=T)
112 write.table(UNPROD, "allUnproductive.csv", sep=",",quote=F,row.names=F,col.names=T)
113
114 #write the samples to a file
115 sampleFile <- file("samples.txt")
116 un = unique(inputdata$Sample)
117 un = paste(un, sep="\n")
118 writeLines(un, sampleFile)
119 close(sampleFile)
120
121 # ---------------------- Counting the productive/unproductive and unique sequences ----------------------
122
123 print("Report Clonality - counting productive/unproductive/unique")
124
125 if(!("Functionality" %in% inputdata)){ #add a functionality column to the igblast data
126 inputdata$Functionality = "unproductive"
127 search = (inputdata$VDJ.Frame != "In-frame with stop codon" & inputdata$VDJ.Frame != "Out-of-frame" & inputdata$CDR3.Found.How != "NOT_FOUND")
128 if(sum(search) > 0){
129 inputdata[search,]$Functionality = "productive"
130 }
131 }
132
133 inputdata.dt = data.table(inputdata) #for speed
134
135 if(clonaltype == "none"){
136 ct = c("clonaltype")
137 }
138
139 inputdata.dt$samples_replicates = paste(inputdata.dt$Sample, inputdata.dt$Replicate, sep="_")
140 samples_replicates = c(unique(inputdata.dt$samples_replicates), unique(as.character(inputdata.dt$Sample)))
141 frequency_table = data.frame(ID = samples_replicates[order(samples_replicates)])
142
143
144 sample_productive_count = inputdata.dt[, list(All=.N,
145 Productive = nrow(.SD[.SD$Functionality == "productive" | .SD$Functionality == "productive (see comment)",]),
146 perc_prod = 1,
147 Productive_unique = nrow(.SD[.SD$Functionality == "productive" | .SD$Functionality == "productive (see comment)",list(count=.N),by=ct]),
148 perc_prod_un = 1,
149 Unproductive= nrow(.SD[.SD$Functionality != "productive" & .SD$Functionality != "productive (see comment)",]),
150 perc_unprod = 1,
151 Unproductive_unique =nrow(.SD[.SD$Functionality != "productive" & .SD$Functionality != "productive (see comment)",list(count=.N),by=ct]),
152 perc_unprod_un = 1),
153 by=c("Sample")]
154
155 sample_productive_count$perc_prod = round(sample_productive_count$Productive / sample_productive_count$All * 100)
156 sample_productive_count$perc_prod_un = round(sample_productive_count$Productive_unique / sample_productive_count$All * 100)
157
158 sample_productive_count$perc_unprod = round(sample_productive_count$Unproductive / sample_productive_count$All * 100)
159 sample_productive_count$perc_unprod_un = round(sample_productive_count$Unproductive_unique / sample_productive_count$All * 100)
160
161 sample_replicate_productive_count = inputdata.dt[, list(All=.N,
162 Productive = nrow(.SD[.SD$Functionality == "productive" | .SD$Functionality == "productive (see comment)",]),
163 perc_prod = 1,
164 Productive_unique = nrow(.SD[.SD$Functionality == "productive" | .SD$Functionality == "productive (see comment)",list(count=.N),by=ct]),
165 perc_prod_un = 1,
166 Unproductive= nrow(.SD[.SD$Functionality != "productive" & .SD$Functionality != "productive (see comment)",]),
167 perc_unprod = 1,
168 Unproductive_unique =nrow(.SD[.SD$Functionality != "productive" & .SD$Functionality != "productive (see comment)",list(count=.N),by=ct]),
169 perc_unprod_un = 1),
170 by=c("samples_replicates")]
171
172 sample_replicate_productive_count$perc_prod = round(sample_replicate_productive_count$Productive / sample_replicate_productive_count$All * 100)
173 sample_replicate_productive_count$perc_prod_un = round(sample_replicate_productive_count$Productive_unique / sample_replicate_productive_count$All * 100)
174
175 sample_replicate_productive_count$perc_unprod = round(sample_replicate_productive_count$Unproductive / sample_replicate_productive_count$All * 100)
176 sample_replicate_productive_count$perc_unprod_un = round(sample_replicate_productive_count$Unproductive_unique / sample_replicate_productive_count$All * 100)
177
178 setnames(sample_replicate_productive_count, colnames(sample_productive_count))
179
180 counts = rbind(sample_replicate_productive_count, sample_productive_count)
181 counts = counts[order(counts$Sample),]
182
183 write.table(x=counts, file="productive_counting.txt", sep=",",quote=F,row.names=F,col.names=F)
184
185 # ---------------------- Frequency calculation for V, D and J ----------------------
186
187 print("Report Clonality - frequency calculation V, D and J")
188
189 PRODFV = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Sample", "Top.V.Gene")])
190 Total = ddply(PRODFV, .(Sample), function(x) data.frame(Total = sum(x$Length)))
191 PRODFV = merge(PRODFV, Total, by.x='Sample', by.y='Sample', all.x=TRUE)
192 PRODFV = ddply(PRODFV, c("Sample", "Top.V.Gene"), summarise, relFreq= (Length*100 / Total))
193
194 PRODFD = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Sample", "Top.D.Gene")])
195 Total = ddply(PRODFD, .(Sample), function(x) data.frame(Total = sum(x$Length)))
196 PRODFD = merge(PRODFD, Total, by.x='Sample', by.y='Sample', all.x=TRUE)
197 PRODFD = ddply(PRODFD, c("Sample", "Top.D.Gene"), summarise, relFreq= (Length*100 / Total))
198
199 PRODFJ = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Sample", "Top.J.Gene")])
200 Total = ddply(PRODFJ, .(Sample), function(x) data.frame(Total = sum(x$Length)))
201 PRODFJ = merge(PRODFJ, Total, by.x='Sample', by.y='Sample', all.x=TRUE)
202 PRODFJ = ddply(PRODFJ, c("Sample", "Top.J.Gene"), summarise, relFreq= (Length*100 / Total))
203
204 # ---------------------- Setting up the gene names for the different species/loci ----------------------
205
206 print("Report Clonality - getting genes for species/loci")
207
208 Vchain = ""
209 Dchain = ""
210 Jchain = ""
211
212 if(species == "custom"){
213 print("Custom genes: ")
214 splt = unlist(strsplit(locus, ";"))
215 print(paste("V:", splt[1]))
216 print(paste("D:", splt[2]))
217 print(paste("J:", splt[3]))
218
219 Vchain = unlist(strsplit(splt[1], ","))
220 Vchain = data.frame(v.name = Vchain, chr.orderV = 1:length(Vchain))
221
222 Dchain = unlist(strsplit(splt[2], ","))
223 if(length(Dchain) > 0){
224 Dchain = data.frame(v.name = Dchain, chr.orderD = 1:length(Dchain))
225 } else {
226 Dchain = data.frame(v.name = character(0), chr.orderD = numeric(0))
227 }
228
229 Jchain = unlist(strsplit(splt[3], ","))
230 Jchain = data.frame(v.name = Jchain, chr.orderJ = 1:length(Jchain))
231
232 } else {
233 genes = read.table("genes.txt", sep="\t", header=TRUE, fill=T, comment.char="")
234
235 Vchain = genes[grepl(species, genes$Species) & genes$locus == locus & genes$region == "V",c("IMGT.GENE.DB", "chr.order")]
236 colnames(Vchain) = c("v.name", "chr.orderV")
237 Dchain = genes[grepl(species, genes$Species) & genes$locus == locus & genes$region == "D",c("IMGT.GENE.DB", "chr.order")]
238 colnames(Dchain) = c("v.name", "chr.orderD")
239 Jchain = genes[grepl(species, genes$Species) & genes$locus == locus & genes$region == "J",c("IMGT.GENE.DB", "chr.order")]
240 colnames(Jchain) = c("v.name", "chr.orderJ")
241 }
242 useD = TRUE
243 if(nrow(Dchain) == 0){
244 useD = FALSE
245 cat("No D Genes in this species/locus")
246 }
247 print(paste(nrow(Vchain), "genes in V"))
248 print(paste(nrow(Dchain), "genes in D"))
249 print(paste(nrow(Jchain), "genes in J"))
250
251 # ---------------------- merge with the frequency count ----------------------
252
253 PRODFV = merge(PRODFV, Vchain, by.x='Top.V.Gene', by.y='v.name', all.x=TRUE)
254
255 PRODFD = merge(PRODFD, Dchain, by.x='Top.D.Gene', by.y='v.name', all.x=TRUE)
256
257 PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE)
258
259 # ---------------------- Create the V, D and J frequency plots and write the data.frame for every plot to a file ----------------------
260
261 print("Report Clonality - V, D and J frequency plots")
262
263 pV = ggplot(PRODFV)
264 pV = pV + geom_bar( aes( x=factor(reorder(Top.V.Gene, chr.orderV)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1))
265 pV = pV + xlab("Summary of V gene") + ylab("Frequency") + ggtitle("Relative frequency of V gene usage")
266 write.table(x=PRODFV, file="VFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
267
268 png("VPlot.png",width = 1280, height = 720)
269 pV
270 dev.off();
271
272 if(useD){
273 pD = ggplot(PRODFD)
274 pD = pD + geom_bar( aes( x=factor(reorder(Top.D.Gene, chr.orderD)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1))
275 pD = pD + xlab("Summary of D gene") + ylab("Frequency") + ggtitle("Relative frequency of D gene usage")
276 write.table(x=PRODFD, file="DFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
277
278 png("DPlot.png",width = 800, height = 600)
279 print(pD)
280 dev.off();
281 }
282
283 pJ = ggplot(PRODFJ)
284 pJ = pJ + geom_bar( aes( x=factor(reorder(Top.J.Gene, chr.orderJ)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1))
285 pJ = pJ + xlab("Summary of J gene") + ylab("Frequency") + ggtitle("Relative frequency of J gene usage")
286 write.table(x=PRODFJ, file="JFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
287
288 png("JPlot.png",width = 800, height = 600)
289 pJ
290 dev.off();
291
292 pJ = ggplot(PRODFJ)
293 pJ = pJ + geom_bar( aes( x=factor(reorder(Top.J.Gene, chr.orderJ)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1))
294 pJ = pJ + xlab("Summary of J gene") + ylab("Frequency") + ggtitle("Relative frequency of J gene usage")
295 write.table(x=PRODFJ, file="JFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
296
297 png("JPlot.png",width = 800, height = 600)
298 pJ
299 dev.off();
300
301 # ---------------------- Now the frequency plots of the V, D and J families ----------------------
302
303 print("Report Clonality - V, D and J family plots")
304
305 VGenes = PRODF[,c("Sample", "Top.V.Gene")]
306 VGenes$Top.V.Gene = gsub("-.*", "", VGenes$Top.V.Gene)
307 VGenes = data.frame(data.table(VGenes)[, list(Count=.N), by=c("Sample", "Top.V.Gene")])
308 TotalPerSample = data.frame(data.table(VGenes)[, list(total=sum(.SD$Count)), by=Sample])
309 VGenes = merge(VGenes, TotalPerSample, by="Sample")
310 VGenes$Frequency = VGenes$Count * 100 / VGenes$total
311 VPlot = ggplot(VGenes)
312 VPlot = VPlot + geom_bar(aes( x = Top.V.Gene, y = Frequency, fill = Sample), stat='identity', position='dodge' ) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
313 ggtitle("Distribution of V gene families") +
314 ylab("Percentage of sequences")
315 png("VFPlot.png")
316 VPlot
317 dev.off();
318 write.table(x=VGenes, file="VFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
319
320 if(useD){
321 DGenes = PRODF[,c("Sample", "Top.D.Gene")]
322 DGenes$Top.D.Gene = gsub("-.*", "", DGenes$Top.D.Gene)
323 DGenes = data.frame(data.table(DGenes)[, list(Count=.N), by=c("Sample", "Top.D.Gene")])
324 TotalPerSample = data.frame(data.table(DGenes)[, list(total=sum(.SD$Count)), by=Sample])
325 DGenes = merge(DGenes, TotalPerSample, by="Sample")
326 DGenes$Frequency = DGenes$Count * 100 / DGenes$total
327 DPlot = ggplot(DGenes)
328 DPlot = DPlot + geom_bar(aes( x = Top.D.Gene, y = Frequency, fill = Sample), stat='identity', position='dodge' ) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
329 ggtitle("Distribution of D gene families") +
330 ylab("Percentage of sequences")
331 png("DFPlot.png")
332 print(DPlot)
333 dev.off();
334 write.table(x=DGenes, file="DFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
335 }
336
337 JGenes = PRODF[,c("Sample", "Top.J.Gene")]
338 JGenes$Top.J.Gene = gsub("-.*", "", JGenes$Top.J.Gene)
339 JGenes = data.frame(data.table(JGenes)[, list(Count=.N), by=c("Sample", "Top.J.Gene")])
340 TotalPerSample = data.frame(data.table(JGenes)[, list(total=sum(.SD$Count)), by=Sample])
341 JGenes = merge(JGenes, TotalPerSample, by="Sample")
342 JGenes$Frequency = JGenes$Count * 100 / JGenes$total
343 JPlot = ggplot(JGenes)
344 JPlot = JPlot + geom_bar(aes( x = Top.J.Gene, y = Frequency, fill = Sample), stat='identity', position='dodge' ) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
345 ggtitle("Distribution of J gene families") +
346 ylab("Percentage of sequences")
347 png("JFPlot.png")
348 JPlot
349 dev.off();
350 write.table(x=JGenes, file="JFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
351
352 # ---------------------- Plotting the cdr3 length ----------------------
353
354 print("Report Clonality - CDR3 length plot")
355
356 CDR3Length = data.frame(data.table(PRODF)[, list(Count=.N), by=c("Sample", "CDR3.Length.DNA")])
357 TotalPerSample = data.frame(data.table(CDR3Length)[, list(total=sum(.SD$Count)), by=Sample])
358 CDR3Length = merge(CDR3Length, TotalPerSample, by="Sample")
359 CDR3Length$Frequency = CDR3Length$Count * 100 / CDR3Length$total
360 CDR3LengthPlot = ggplot(CDR3Length)
361 CDR3LengthPlot = CDR3LengthPlot + geom_bar(aes( x = CDR3.Length.DNA, y = Frequency, fill = Sample), stat='identity', position='dodge' ) + theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
362 ggtitle("Length distribution of CDR3") +
363 xlab("CDR3 Length") +
364 ylab("Percentage of sequences")
365 png("CDR3LengthPlot.png",width = 1280, height = 720)
366 CDR3LengthPlot
367 dev.off()
368 write.table(x=CDR3Length, file="CDR3LengthPlot.csv", sep=",",quote=F,row.names=F,col.names=T)
369
370 # ---------------------- Plot the heatmaps ----------------------
371
372 #get the reverse order for the V and D genes
373 revVchain = Vchain
374 revDchain = Dchain
375 revVchain$chr.orderV = rev(revVchain$chr.orderV)
376 revDchain$chr.orderD = rev(revDchain$chr.orderD)
377
378 if(useD){
379 print("Report Clonality - Heatmaps VD")
380 plotVD <- function(dat){
381 if(length(dat[,1]) == 0){
382 return()
383 }
384 img = ggplot() +
385 geom_tile(data=dat, aes(x=factor(reorder(Top.D.Gene, chr.orderD)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) +
386 theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
387 scale_fill_gradient(low="gold", high="blue", na.value="white") +
388 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) +
389 xlab("D genes") +
390 ylab("V Genes")
391
392 png(paste("HeatmapVD_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Dchain$v.name)), height=100+(15*length(Vchain$v.name)))
393 print(img)
394 dev.off()
395 write.table(x=acast(dat, Top.V.Gene~Top.D.Gene, value.var="Length"), file=paste("HeatmapVD_", unique(dat[3])[1,1], ".csv", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
396 }
397
398 VandDCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.V.Gene", "Top.D.Gene", "Sample")])
399
400 VandDCount$l = log(VandDCount$Length)
401 maxVD = data.frame(data.table(VandDCount)[, list(max=max(l)), by=c("Sample")])
402 VandDCount = merge(VandDCount, maxVD, by.x="Sample", by.y="Sample", all.x=T)
403 VandDCount$relLength = VandDCount$l / VandDCount$max
404
405 cartegianProductVD = expand.grid(Top.V.Gene = Vchain$v.name, Top.D.Gene = Dchain$v.name, Sample = unique(inputdata$Sample))
406
407 completeVD = merge(VandDCount, cartegianProductVD, all.y=TRUE)
408 completeVD = merge(completeVD, revVchain, by.x="Top.V.Gene", by.y="v.name", all.x=TRUE)
409 completeVD = merge(completeVD, Dchain, by.x="Top.D.Gene", by.y="v.name", all.x=TRUE)
410
411 fltr = is.nan(completeVD$relLength)
412 if(any(fltr)){
413 completeVD[fltr,"relLength"] = 1
414 }
415
416 VDList = split(completeVD, f=completeVD[,"Sample"])
417 lapply(VDList, FUN=plotVD)
418 }
419
420 print("Report Clonality - Heatmaps VJ")
421
422 plotVJ <- function(dat){
423 if(length(dat[,1]) == 0){
424 return()
425 }
426 cat(paste(unique(dat[3])[1,1]))
427 img = ggplot() +
428 geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) +
429 theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
430 scale_fill_gradient(low="gold", high="blue", na.value="white") +
431 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) +
432 xlab("J genes") +
433 ylab("V Genes")
434
435 png(paste("HeatmapVJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Vchain$v.name)))
436 print(img)
437 dev.off()
438 write.table(x=acast(dat, Top.V.Gene~Top.J.Gene, value.var="Length"), file=paste("HeatmapVJ_", unique(dat[3])[1,1], ".csv", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
439 }
440
441 VandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.V.Gene", "Top.J.Gene", "Sample")])
442
443 VandJCount$l = log(VandJCount$Length)
444 maxVJ = data.frame(data.table(VandJCount)[, list(max=max(l)), by=c("Sample")])
445 VandJCount = merge(VandJCount, maxVJ, by.x="Sample", by.y="Sample", all.x=T)
446 VandJCount$relLength = VandJCount$l / VandJCount$max
447
448 cartegianProductVJ = expand.grid(Top.V.Gene = Vchain$v.name, Top.J.Gene = Jchain$v.name, Sample = unique(inputdata$Sample))
449
450 completeVJ = merge(VandJCount, cartegianProductVJ, all.y=TRUE)
451 completeVJ = merge(completeVJ, revVchain, by.x="Top.V.Gene", by.y="v.name", all.x=TRUE)
452 completeVJ = merge(completeVJ, Jchain, by.x="Top.J.Gene", by.y="v.name", all.x=TRUE)
453
454 fltr = is.nan(completeVJ$relLength)
455 if(any(fltr)){
456 completeVJ[fltr,"relLength"] = 1
457 }
458
459 VJList = split(completeVJ, f=completeVJ[,"Sample"])
460 lapply(VJList, FUN=plotVJ)
461
462
463
464 if(useD){
465 print("Report Clonality - Heatmaps DJ")
466 plotDJ <- function(dat){
467 if(length(dat[,1]) == 0){
468 return()
469 }
470 img = ggplot() +
471 geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.D.Gene, chr.orderD)), fill=relLength)) +
472 theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
473 scale_fill_gradient(low="gold", high="blue", na.value="white") +
474 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) +
475 xlab("J genes") +
476 ylab("D Genes")
477
478 png(paste("HeatmapDJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Dchain$v.name)))
479 print(img)
480 dev.off()
481 write.table(x=acast(dat, Top.D.Gene~Top.J.Gene, value.var="Length"), file=paste("HeatmapDJ_", unique(dat[3])[1,1], ".csv", sep=""), sep=",",quote=F,row.names=T,col.names=NA)
482 }
483
484
485 DandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.D.Gene", "Top.J.Gene", "Sample")])
486
487 DandJCount$l = log(DandJCount$Length)
488 maxDJ = data.frame(data.table(DandJCount)[, list(max=max(l)), by=c("Sample")])
489 DandJCount = merge(DandJCount, maxDJ, by.x="Sample", by.y="Sample", all.x=T)
490 DandJCount$relLength = DandJCount$l / DandJCount$max
491
492 cartegianProductDJ = expand.grid(Top.D.Gene = Dchain$v.name, Top.J.Gene = Jchain$v.name, Sample = unique(inputdata$Sample))
493
494 completeDJ = merge(DandJCount, cartegianProductDJ, all.y=TRUE)
495 completeDJ = merge(completeDJ, revDchain, by.x="Top.D.Gene", by.y="v.name", all.x=TRUE)
496 completeDJ = merge(completeDJ, Jchain, by.x="Top.J.Gene", by.y="v.name", all.x=TRUE)
497
498 fltr = is.nan(completeDJ$relLength)
499 if(any(fltr)){
500 completeDJ[fltr, "relLength"] = 1
501 }
502
503 DJList = split(completeDJ, f=completeDJ[,"Sample"])
504 lapply(DJList, FUN=plotDJ)
505 }
506
507
508 # ---------------------- output tables for the circos plots ----------------------
509
510 print("Report Clonality - Circos data")
511
512 for(smpl in unique(PRODF$Sample)){
513 PRODF.sample = PRODF[PRODF$Sample == smpl,]
514
515 fltr = PRODF.sample$Top.V.Gene == ""
516 if(any(fltr, na.rm=T)){
517 PRODF.sample[fltr, "Top.V.Gene"] = "NA"
518 }
519
520 fltr = PRODF.sample$Top.D.Gene == ""
521 if(any(fltr, na.rm=T)){
522 PRODF.sample[fltr, "Top.D.Gene"] = "NA"
523 }
524
525 fltr = PRODF.sample$Top.J.Gene == ""
526 if(any(fltr, na.rm=T)){
527 PRODF.sample[fltr, "Top.J.Gene"] = "NA"
528 }
529
530 v.d = table(PRODF.sample$Top.V.Gene, PRODF.sample$Top.D.Gene)
531 v.j = table(PRODF.sample$Top.V.Gene, PRODF.sample$Top.J.Gene)
532 d.j = table(PRODF.sample$Top.D.Gene, PRODF.sample$Top.J.Gene)
533
534 write.table(v.d, file=paste(smpl, "_VD_circos.txt", sep=""), sep="\t", quote=F, row.names=T, col.names=NA)
535 write.table(v.j, file=paste(smpl, "_VJ_circos.txt", sep=""), sep="\t", quote=F, row.names=T, col.names=NA)
536 write.table(d.j, file=paste(smpl, "_DJ_circos.txt", sep=""), sep="\t", quote=F, row.names=T, col.names=NA)
537 }
538
539 # ---------------------- calculating the clonality score ----------------------
540
541 if("Replicate" %in% colnames(inputdata)) #can only calculate clonality score when replicate information is available
542 {
543 print("Report Clonality - Clonality")
544 if(clonality_method == "boyd"){
545 samples = split(clonalityFrame, clonalityFrame$Sample, drop=T)
546
547 for (sample in samples){
548 res = data.frame(paste=character(0))
549 sample_id = unique(sample$Sample)[[1]]
550 for(replicate in unique(sample$Replicate)){
551 tmp = sample[sample$Replicate == replicate,]
552 clone_table = data.frame(table(tmp$clonaltype))
553 clone_col_name = paste("V", replicate, sep="")
554 colnames(clone_table) = c("paste", clone_col_name)
555 res = merge(res, clone_table, by="paste", all=T)
556 }
557
558 res[is.na(res)] = 0
559 infer.result = infer.clonality(as.matrix(res[,2:ncol(res)]))
560
561 write.table(data.table(infer.result[[12]]), file=paste("lymphclon_clonality_", sample_id, ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=F)
562
563 res$type = rowSums(res[,2:ncol(res)])
564
565 coincidence.table = data.frame(table(res$type))
566 colnames(coincidence.table) = c("Coincidence Type", "Raw Coincidence Freq")
567 write.table(coincidence.table, file=paste("lymphclon_coincidences_", sample_id, ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=T)
568 }
569 } else {
570 write.table(clonalityFrame, "clonalityComplete.csv", sep=",",quote=F,row.names=F,col.names=T)
571
572 clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "clonaltype")])
573 clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")])
574 clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count
575 clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")])
576 clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample")
577
578 ct = c('Type\tWeight\n2\t1\n3\t3\n4\t6\n5\t10\n6\t15')
579 tcct = textConnection(ct)
580 CT = read.table(tcct, sep="\t", header=TRUE)
581 close(tcct)
582 clonalFreqCount = merge(clonalFreqCount, CT, by.x="Type", by.y="Type", all.x=T)
583 clonalFreqCount$WeightedCount = clonalFreqCount$Count * clonalFreqCount$Weight
584
585 ReplicateReads = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "Replicate", "clonaltype")])
586 ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")])
587 clonalFreqCount$Reads = as.numeric(clonalFreqCount$Reads)
588 ReplicateReads$Reads = as.numeric(ReplicateReads$Reads)
589 ReplicateReads$squared = as.numeric(ReplicateReads$Reads * ReplicateReads$Reads)
590
591 ReplicatePrint <- function(dat){
592 write.table(dat[-1], paste("ReplicateReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
593 }
594
595 ReplicateSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"])
596 lapply(ReplicateSplit, FUN=ReplicatePrint)
597
598 ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(as.numeric(Reads)), ReadsSquaredSum=sum(as.numeric(squared))), by=c("Sample")])
599 clonalFreqCount = merge(clonalFreqCount, ReplicateReads, by.x="Sample", by.y="Sample", all.x=T)
600
601 ReplicateSumPrint <- function(dat){
602 write.table(dat[-1], paste("ReplicateSumReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
603 }
604
605 ReplicateSumSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"])
606 lapply(ReplicateSumSplit, FUN=ReplicateSumPrint)
607
608 clonalFreqCountSum = data.frame(data.table(clonalFreqCount)[, list(Numerator=sum(WeightedCount, na.rm=T)), by=c("Sample")])
609 clonalFreqCount = merge(clonalFreqCount, clonalFreqCountSum, by.x="Sample", by.y="Sample", all.x=T)
610 clonalFreqCount$ReadsSum = as.numeric(clonalFreqCount$ReadsSum) #prevent integer overflow
611 clonalFreqCount$Denominator = (((clonalFreqCount$ReadsSum * clonalFreqCount$ReadsSum) - clonalFreqCount$ReadsSquaredSum) / 2)
612 clonalFreqCount$Result = (clonalFreqCount$Numerator + 1) / (clonalFreqCount$Denominator + 1)
613
614 ClonalityScorePrint <- function(dat){
615 write.table(dat$Result, paste("ClonalityScore_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
616 }
617
618 clonalityScore = clonalFreqCount[c("Sample", "Result")]
619 clonalityScore = unique(clonalityScore)
620
621 clonalityScoreSplit = split(clonalityScore, f=clonalityScore[,"Sample"])
622 lapply(clonalityScoreSplit, FUN=ClonalityScorePrint)
623
624 clonalityOverview = clonalFreqCount[c("Sample", "Type", "Count", "Weight", "WeightedCount")]
625
626
627
628 ClonalityOverviewPrint <- function(dat){
629 write.table(dat[-1], paste("ClonalityOverView_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
630 }
631
632 clonalityOverviewSplit = split(clonalityOverview, f=clonalityOverview$Sample)
633 lapply(clonalityOverviewSplit, FUN=ClonalityOverviewPrint)
634 }
635 }
636
637 imgtcolumns = c("X3V.REGION.trimmed.nt.nb","P3V.nt.nb", "N1.REGION.nt.nb", "P5D.nt.nb", "X5D.REGION.trimmed.nt.nb", "X3D.REGION.trimmed.nt.nb", "P3D.nt.nb", "N2.REGION.nt.nb", "P5J.nt.nb", "X5J.REGION.trimmed.nt.nb", "X3V.REGION.trimmed.nt.nb", "X5D.REGION.trimmed.nt.nb", "X3D.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb", "N1.REGION.nt.nb", "N2.REGION.nt.nb", "P3V.nt.nb", "P5D.nt.nb", "P3D.nt.nb", "P5J.nt.nb")
638 if(all(imgtcolumns %in% colnames(inputdata)))
639 {
640 print("found IMGT columns, running junction analysis")
641
642 if(locus %in% c("IGK","IGL", "TRA", "TRG")){
643 print("VJ recombination, using N column for junction analysis")
644 print(names(PRODF))
645 print(head(PRODF$N.REGION.nt.nb, 30))
646 PRODF$N1.REGION.nt.nb = PRODF$N.REGION.nt.nb
647 }
648
649 num_median = function(x, na.rm) { as.numeric(median(x, na.rm=na.rm)) }
650
651 newData = data.frame(data.table(PRODF)[,list(unique=.N,
652 VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),
653 P1=mean(.SD$P3V.nt.nb, na.rm=T),
654 N1=mean(.SD$N1.REGION.nt.nb, na.rm=T),
655 P2=mean(.SD$P5D.nt.nb, na.rm=T),
656 DEL.DH=mean(.SD$X5D.REGION.trimmed.nt.nb, na.rm=T),
657 DH.DEL=mean(.SD$X3D.REGION.trimmed.nt.nb, na.rm=T),
658 P3=mean(.SD$P3D.nt.nb, na.rm=T),
659 N2=mean(.SD$N2.REGION.nt.nb, na.rm=T),
660 P4=mean(.SD$P5J.nt.nb, na.rm=T),
661 DEL.JH=mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T),
662 Total.Del=( mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T) +
663 mean(.SD$X5D.REGION.trimmed.nt.nb, na.rm=T) +
664 mean(.SD$X3D.REGION.trimmed.nt.nb, na.rm=T) +
665 mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T)),
666
667 Total.N=( mean(.SD$N1.REGION.nt.nb, na.rm=T) +
668 mean(.SD$N2.REGION.nt.nb, na.rm=T)),
669
670 Total.P=( mean(.SD$P3V.nt.nb, na.rm=T) +
671 mean(.SD$P5D.nt.nb, na.rm=T) +
672 mean(.SD$P3D.nt.nb, na.rm=T) +
673 mean(.SD$P5J.nt.nb, na.rm=T))),
674 by=c("Sample")])
675 newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1)
676 write.table(newData, "junctionAnalysisProd_mean.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F)
677
678 newData = data.frame(data.table(PRODF)[,list(unique=.N,
679 VH.DEL=num_median(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),
680 P1=num_median(.SD$P3V.nt.nb, na.rm=T),
681 N1=num_median(.SD$N1.REGION.nt.nb, na.rm=T),
682 P2=num_median(.SD$P5D.nt.nb, na.rm=T),
683 DEL.DH=num_median(.SD$X5D.REGION.trimmed.nt.nb, na.rm=T),
684 DH.DEL=num_median(.SD$X3D.REGION.trimmed.nt.nb, na.rm=T),
685 P3=num_median(.SD$P3D.nt.nb, na.rm=T),
686 N2=num_median(.SD$N2.REGION.nt.nb, na.rm=T),
687 P4=num_median(.SD$P5J.nt.nb, na.rm=T),
688 DEL.JH=num_median(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T),
689 Total.Del=num_median(c(.SD$X3V.REGION.trimmed.nt.nb,
690 .SD$X5D.REGION.trimmed.nt.nb,
691 .SD$X3D.REGION.trimmed.nt.nb,
692 .SD$X5J.REGION.trimmed.nt.nb), na.rm=T),
693 Total.N=num_median( c(.SD$N1.REGION.nt.nb,
694 .SD$N2.REGION.nt.nb), na.rm=T),
695 Total.P=num_median(c(.SD$P3V.nt.nb,
696 .SD$P5D.nt.nb,
697 .SD$P3D.nt.nb,
698 .SD$P5J.nt.nb), na.rm=T)),
699 by=c("Sample")])
700 newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1)
701 write.table(newData, "junctionAnalysisProd_median.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F)
702
703 newData = data.frame(data.table(UNPROD)[,list(unique=.N,
704 VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),
705 P1=mean(.SD$P3V.nt.nb, na.rm=T),
706 N1=mean(.SD$N1.REGION.nt.nb, na.rm=T),
707 P2=mean(.SD$P5D.nt.nb, na.rm=T),
708 DEL.DH=mean(.SD$X5D.REGION.trimmed.nt.nb, na.rm=T),
709 DH.DEL=mean(.SD$X3D.REGION.trimmed.nt.nb, na.rm=T),
710 P3=mean(.SD$P3D.nt.nb, na.rm=T),
711 N2=mean(.SD$N2.REGION.nt.nb, na.rm=T),
712 P4=mean(.SD$P5J.nt.nb, na.rm=T),
713 DEL.JH=mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T),
714 Total.Del=(mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T) +
715 mean(.SD$X5D.REGION.trimmed.nt.nb, na.rm=T) +
716 mean(.SD$X3D.REGION.trimmed.nt.nb, na.rm=T) +
717 mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T)),
718 Total.N=( mean(.SD$N1.REGION.nt.nb, na.rm=T) +
719 mean(.SD$N2.REGION.nt.nb, na.rm=T)),
720 Total.P=( mean(.SD$P3V.nt.nb, na.rm=T) +
721 mean(.SD$P5D.nt.nb, na.rm=T) +
722 mean(.SD$P3D.nt.nb, na.rm=T) +
723 mean(.SD$P5J.nt.nb, na.rm=T))),
724 by=c("Sample")])
725 newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1)
726 write.table(newData, "junctionAnalysisUnProd_mean.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F)
727
728 newData = data.frame(data.table(UNPROD)[,list(unique=.N,
729 VH.DEL=num_median(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),
730 P1=num_median(.SD$P3V.nt.nb, na.rm=T),
731 N1=num_median(.SD$N1.REGION.nt.nb, na.rm=T),
732 P2=num_median(.SD$P5D.nt.nb, na.rm=T),
733 DEL.DH=num_median(.SD$X5D.REGION.trimmed.nt.nb, na.rm=T),
734 DH.DEL=num_median(.SD$X3D.REGION.trimmed.nt.nb, na.rm=T),
735 P3=num_median(.SD$P3D.nt.nb, na.rm=T),
736 N2=num_median(.SD$N2.REGION.nt.nb, na.rm=T),
737 P4=num_median(.SD$P5J.nt.nb, na.rm=T),
738 DEL.JH=num_median(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T),
739 Total.Del=num_median(c(.SD$X3V.REGION.trimmed.nt.nb,
740 .SD$X5D.REGION.trimmed.nt.nb,
741 .SD$X3D.REGION.trimmed.nt.nb,
742 .SD$X5J.REGION.trimmed.nt.nb), na.rm=T),
743 Total.N=num_median( c(.SD$N1.REGION.nt.nb,
744 .SD$N2.REGION.nt.nb), na.rm=T),
745 Total.P=num_median(c(.SD$P3V.nt.nb,
746 .SD$P5D.nt.nb,
747 .SD$P3D.nt.nb,
748 .SD$P5J.nt.nb), na.rm=T)),
749 by=c("Sample")])
750
751 newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1)
752 write.table(newData, "junctionAnalysisUnProd_median.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F)
753 }
754
755 # ---------------------- AA composition in CDR3 ----------------------
756
757 AACDR3 = PRODF[,c("Sample", "CDR3.Seq")]
758
759 TotalPerSample = data.frame(data.table(AACDR3)[, list(total=sum(nchar(as.character(.SD$CDR3.Seq)))), by=Sample])
760
761 AAfreq = list()
762
763 for(i in 1:nrow(TotalPerSample)){
764 sample = TotalPerSample$Sample[i]
765 AAfreq[[i]] = data.frame(table(unlist(strsplit(as.character(AACDR3[AACDR3$Sample == sample,c("CDR3.Seq")]), ""))))
766 AAfreq[[i]]$Sample = sample
767 }
768
769 AAfreq = ldply(AAfreq, data.frame)
770 AAfreq = merge(AAfreq, TotalPerSample, by="Sample", all.x = T)
771 AAfreq$freq_perc = as.numeric(AAfreq$Freq / AAfreq$total * 100)
772
773
774 AAorder = read.table(sep="\t", header=TRUE, text="order.aa\tAA\n1\tR\n2\tK\n3\tN\n4\tD\n5\tQ\n6\tE\n7\tH\n8\tP\n9\tY\n10\tW\n11\tS\n12\tT\n13\tG\n14\tA\n15\tM\n16\tC\n17\tF\n18\tL\n19\tV\n20\tI")
775 AAfreq = merge(AAfreq, AAorder, by.x='Var1', by.y='AA', all.x=TRUE)
776
777 AAfreq = AAfreq[!is.na(AAfreq$order.aa),]
778
779 AAfreqplot = ggplot(AAfreq)
780 AAfreqplot = AAfreqplot + geom_bar(aes( x=factor(reorder(Var1, order.aa)), y = freq_perc, fill = Sample), stat='identity', position='dodge' )
781 AAfreqplot = AAfreqplot + annotate("rect", xmin = 0.5, xmax = 2.5, ymin = 0, ymax = Inf, fill = "red", alpha = 0.2)
782 AAfreqplot = AAfreqplot + annotate("rect", xmin = 3.5, xmax = 4.5, ymin = 0, ymax = Inf, fill = "blue", alpha = 0.2)
783 AAfreqplot = AAfreqplot + annotate("rect", xmin = 5.5, xmax = 6.5, ymin = 0, ymax = Inf, fill = "blue", alpha = 0.2)
784 AAfreqplot = AAfreqplot + annotate("rect", xmin = 6.5, xmax = 7.5, ymin = 0, ymax = Inf, fill = "red", alpha = 0.2)
785 AAfreqplot = AAfreqplot + ggtitle("Amino Acid Composition in the CDR3") + xlab("Amino Acid, from Hydrophilic (left) to Hydrophobic (right)") + ylab("Percentage")
786
787 png("AAComposition.png",width = 1280, height = 720)
788 AAfreqplot
789 dev.off()
790 write.table(AAfreq, "AAComposition.csv" , sep=",",quote=F,na="-",row.names=F,col.names=T)
791
792