Mercurial > repos > davidvanzessen > complete_immunerepertoire_igg
comparison RScript_t.r @ 1:778a9d130904 draft
Uploaded
author | davidvanzessen |
---|---|
date | Thu, 04 Sep 2014 07:46:23 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
0:7d97fa9a0423 | 1:778a9d130904 |
---|---|
1 #options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) | |
2 | |
3 args <- commandArgs(trailingOnly = TRUE) | |
4 | |
5 inFile = args[1] | |
6 outFile = args[2] | |
7 outDir = args[3] | |
8 clonalType = args[4] | |
9 species = args[5] | |
10 locus = args[6] | |
11 selection = args[7] | |
12 | |
13 if (!("gridExtra" %in% rownames(installed.packages()))) { | |
14 install.packages("gridExtra", repos="http://cran.xl-mirror.nl/") | |
15 } | |
16 library(gridExtra) | |
17 if (!("ggplot2" %in% rownames(installed.packages()))) { | |
18 install.packages("ggplot2", repos="http://cran.xl-mirror.nl/") | |
19 } | |
20 library(ggplot2) | |
21 if (!("plyr" %in% rownames(installed.packages()))) { | |
22 install.packages("plyr", repos="http://cran.xl-mirror.nl/") | |
23 } | |
24 library(plyr) | |
25 | |
26 if (!("data.table" %in% rownames(installed.packages()))) { | |
27 install.packages("data.table", repos="http://cran.xl-mirror.nl/") | |
28 } | |
29 library(data.table) | |
30 | |
31 if (!("reshape2" %in% rownames(installed.packages()))) { | |
32 install.packages("reshape2", repos="http://cran.xl-mirror.nl/") | |
33 } | |
34 library(reshape2) | |
35 | |
36 | |
37 test = read.table(inFile, sep="\t", header=TRUE, fill=T, comment.char="") | |
38 | |
39 test = test[test$Sample != "",] | |
40 | |
41 print("test1\n") | |
42 | |
43 test$Top.V.Gene = gsub("[*]([0-9]+)", "", test$Top.V.Gene) | |
44 test$Top.D.Gene = gsub("[*]([0-9]+)", "", test$Top.D.Gene) | |
45 test$Top.J.Gene = gsub("[*]([0-9]+)", "", test$Top.J.Gene) | |
46 | |
47 #test$VDJCDR3 = do.call(paste, c(test[c("Top.V.Gene", "Top.D.Gene", "Top.J.Gene","CDR3.Seq.DNA")], sep = ":")) | |
48 test$VDJCDR3 = do.call(paste, c(test[unlist(strsplit(clonalType, ","))], sep = ":")) | |
49 | |
50 PROD = test[test$VDJ.Frame != "In-frame with stop codon" & test$VDJ.Frame != "Out-of-frame" & test$CDR3.Found.How != "NOT_FOUND" , ] | |
51 if("Functionality" %in% colnames(test)) { | |
52 PROD = test[test$Functionality == "productive" | test$Functionality == "productive (see comment)", ] | |
53 } | |
54 | |
55 NONPROD = test[test$VDJ.Frame == "In-frame with stop codon" | test$VDJ.Frame == "Out-of-frame" | test$CDR3.Found.How == "NOT_FOUND" , ] | |
56 | |
57 #PRODF = PROD[ -1] | |
58 | |
59 PRODF = PROD | |
60 print("test2\n") | |
61 #PRODF = unique(PRODF) | |
62 if(any(grepl(pattern="_", x=PRODF$ID))){ #dumb and way to simple | |
63 PRODF$freq = gsub("^[0-9]+_", "", PRODF$ID) | |
64 PRODF$freq = gsub("_.*", "", PRODF$freq) | |
65 PRODF$freq = as.numeric(PRODF$freq) | |
66 if(any(is.na(PRODF$freq))){ #fix the dumbness if it fails | |
67 PRODF$freq = 1 | |
68 if(selection == "unique"){ | |
69 PRODF = PRODF[!duplicated(PRODF$VDJCDR3), ] | |
70 } | |
71 } | |
72 } else { | |
73 PRODF$freq = 1 | |
74 if(selection == "unique"){ | |
75 PRODF = PRODF[!duplicated(PRODF$VDJCDR3), ] | |
76 } | |
77 } | |
78 | |
79 PRODFV = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Sample", "Top.V.Gene")]) | |
80 PRODFV$Length = as.numeric(PRODFV$Length) | |
81 Total = 0 | |
82 Total = ddply(PRODFV, .(Sample), function(x) data.frame(Total = sum(x$Length))) | |
83 PRODFV = merge(PRODFV, Total, by.x='Sample', by.y='Sample', all.x=TRUE) | |
84 PRODFV = ddply(PRODFV, c("Sample", "Top.V.Gene"), summarise, relFreq= (Length*100 / Total)) | |
85 | |
86 PRODFD = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Sample", "Top.D.Gene")]) | |
87 PRODFD$Length = as.numeric(PRODFD$Length) | |
88 Total = 0 | |
89 Total = ddply(PRODFD, .(Sample), function(x) data.frame(Total = sum(x$Length))) | |
90 PRODFD = merge(PRODFD, Total, by.x='Sample', by.y='Sample', all.x=TRUE) | |
91 PRODFD = ddply(PRODFD, c("Sample", "Top.D.Gene"), summarise, relFreq= (Length*100 / Total)) | |
92 | |
93 PRODFJ = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Sample", "Top.J.Gene")]) | |
94 PRODFJ$Length = as.numeric(PRODFJ$Length) | |
95 Total = 0 | |
96 Total = ddply(PRODFJ, .(Sample), function(x) data.frame(Total = sum(x$Length))) | |
97 PRODFJ = merge(PRODFJ, Total, by.x='Sample', by.y='Sample', all.x=TRUE) | |
98 PRODFJ = ddply(PRODFJ, c("Sample", "Top.J.Gene"), summarise, relFreq= (Length*100 / Total)) | |
99 | |
100 V = c("v.name\tchr.orderV\n") | |
101 D = c("v.name\tchr.orderD\n") | |
102 J = c("v.name\tchr.orderJ\n") | |
103 | |
104 print("test3\n") | |
105 | |
106 if(species == "human"){ | |
107 if(locus == "trb"){ | |
108 V = c("v.name\tchr.orderV\nTRBV2\t1\nTRBV3-1\t2\nTRBV4-1\t3\nTRBV5-1\t4\nTRBV6-1\t5\nTRBV4-2\t6\nTRBV6-2\t7\nTRBV4-3\t8\nTRBV6-3\t9\nTRBV7-2\t10\nTRBV6-4\t11\nTRBV7-3\t12\nTRBV9\t13\nTRBV10-1\t14\nTRBV11-1\t15\nTRBV10-2\t16\nTRBV11-2\t17\nTRBV6-5\t18\nTRBV7-4\t19\nTRBV5-4\t20\nTRBV6-6\t21\nTRBV5-5\t22\nTRBV7-6\t23\nTRBV5-6\t24\nTRBV6-8\t25\nTRBV7-7\t26\nTRBV6-9\t27\nTRBV7-8\t28\nTRBV5-8\t29\nTRBV7-9\t30\nTRBV13\t31\nTRBV10-3\t32\nTRBV11-3\t33\nTRBV12-3\t34\nTRBV12-4\t35\nTRBV12-5\t36\nTRBV14\t37\nTRBV15\t38\nTRBV16\t39\nTRBV18\t40\nTRBV19\t41\nTRBV20-1\t42\nTRBV24-1\t43\nTRBV25-1\t44\nTRBV27\t45\nTRBV28\t46\nTRBV29-1\t47\nTRBV30\t48") | |
109 D = c("v.name\tchr.orderD\nTRBD1\t1\nTRBD2\t2\n") | |
110 J = c("v.name\tchr.orderJ\nTRBJ1-1\t1\nTRBJ1-2\t2\nTRBJ1-3\t3\nTRBJ1-4\t4\nTRBJ1-5\t5\nTRBJ1-6\t6\nTRBJ2-1\t7\nTRBJ2-2\t8\nTRBJ2-3\t9\nTRBJ2-4\t10\nTRBJ2-5\t11\nTRBJ2-6\t12\nTRBJ2-7\t13") | |
111 } else if (locus == "tra"){ | |
112 V = c("v.name\tchr.orderVTRAV1-1\t1\nTRAV1-2\t2\nTRAV2\t3\nTRAV3\t4\nTRAV4\t5\nTRAV5\t6\nTRAV6\t7\nTRAV7\t8\nTRAV8-1\t9\nTRAV9-1\t10\nTRAV10\t11\nTRAV12-1\t12\nTRAV8-2\t13\nTRAV8-3\t14\nTRAV13-1\t15\nTRAV12-2\t16\nTRAV8-4\t17\nTRAV13-2\t18\nTRAV14/DV4\t19\nTRAV9-2\t20\nTRAV12-3\t21\nTRAV8-6\t22\nTRAV16\t23\nTRAV17\t24\nTRAV18\t25\nTRAV19\t26\nTRAV20\t27\nTRAV21\t28\nTRAV22\t29\nTRAV23/DV6\t30\nTRAV24\t31\nTRAV25\t32\nTRAV26-1\t33\nTRAV27\t34\nTRAV29/DV5\t35\nTRAV30\t36\nTRAV26-2\t37\nTRAV34\t38\nTRAV35\t39\nTRAV36/DV7\t40\nTRAV38-1\t41\nTRAV38-2/DV8\t42\nTRAV39\t43\nTRAV40\t44\nTRAV41\t45\n") | |
113 D = c("v.name\tchr.orderD\n") | |
114 J = c("v.name\tchr.orderJ\nTRAJ57\t1\nTRAJ56\t2\nTRAJ54\t3\nTRAJ53\t4\nTRAJ52\t5\nTRAJ50\t6\nTRAJ49\t7\nTRAJ48\t8\nTRAJ47\t9\nTRAJ46\t10\nTRAJ45\t11\nTRAJ44\t12\nTRAJ43\t13\nTRAJ42\t14\nTRAJ41\t15\nTRAJ40\t16\nTRAJ39\t17\nTRAJ38\t18\nTRAJ37\t19\nTRAJ36\t20\nTRAJ34\t21\nTRAJ33\t22\nTRAJ32\t23\nTRAJ31\t24\nTRAJ30\t25\nTRAJ29\t26\nTRAJ28\t27\nTRAJ27\t28\nTRAJ26\t29\nTRAJ24\t30\nTRAJ23\t31\nTRAJ22\t32\nTRAJ21\t33\nTRAJ20\t34\nTRAJ18\t35\nTRAJ17\t36\nTRAJ16\t37\nTRAJ15\t38\nTRAJ14\t39\nTRAJ13\t40\nTRAJ12\t41\nTRAJ11\t42\nTRAJ10\t43\nTRAJ9\t44\nTRAJ8\t45\nTRAJ7\t46\nTRAJ6\t47\nTRAJ5\t48\nTRAJ4\t49\nTRAJ3\t50") | |
115 } else if (locus == "trg"){ | |
116 V = c("v.name\tchr.orderV\nTRGV9\t1\nTRGV8\t2\nTRGV5\t3\nTRGV4\t4\nTRGV3\t5\nTRGV2\t6") | |
117 D = c("v.name\tchr.orderD\n") | |
118 J = c("v.name\tchr.orderJ\nTRGJ2\t1\nTRGJP2\t2\nTRGJ1\t3\nTRGJP1\t4") | |
119 } else if (locus == "trd"){ | |
120 V = c("v.name\tchr.orderV\nTRDV1\t1\nTRDV2\t2\nTRDV3\t3") | |
121 D = c("v.name\tchr.orderD\nTRDD1\t1\nTRDD2\t2\nTRDD3\t3") | |
122 J = c("v.name\tchr.orderJ\nTRDJ1\t1\nTRDJ4\t2\nTRDJ2\t3\nTRDJ3\t4") | |
123 } | |
124 } else if (species == "mouse"){ | |
125 if(locus == "trb"){ | |
126 cat("mouse trb not yet implemented") | |
127 } else if (locus == "tra"){ | |
128 cat("mouse tra not yet implemented") | |
129 } else if (locus == "trg"){ | |
130 cat("mouse trg not yet implemented") | |
131 } else if (locus == "trd"){ | |
132 cat("mouse trd not yet implemented") | |
133 } | |
134 } | |
135 useD = TRUE | |
136 if(species == "human" && locus == "tra"){ | |
137 useD = FALSE | |
138 cat("No D Genes in this species/locus") | |
139 } | |
140 | |
141 print("test4\n") | |
142 | |
143 tcV = textConnection(V) | |
144 Vchain = read.table(tcV, sep="\t", header=TRUE) | |
145 PRODFV = merge(PRODFV, Vchain, by.x='Top.V.Gene', by.y='v.name', all.x=TRUE) | |
146 close(tcV) | |
147 | |
148 | |
149 tcD = textConnection(D) | |
150 Dchain = read.table(tcD, sep="\t", header=TRUE) | |
151 PRODFD = merge(PRODFD, Dchain, by.x='Top.D.Gene', by.y='v.name', all.x=TRUE) | |
152 close(tcD) | |
153 | |
154 | |
155 | |
156 tcJ = textConnection(J) | |
157 Jchain = read.table(tcJ, sep="\t", header=TRUE) | |
158 PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE) | |
159 close(tcJ) | |
160 | |
161 setwd(outDir) | |
162 | |
163 write.table(PRODF, "allUnique.csv", sep=",",quote=F,row.names=F,col.names=T) | |
164 | |
165 pV = ggplot(PRODFV) | |
166 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)) | |
167 pV = pV + xlab("Summary of V gene") + ylab("Frequency") + ggtitle("Relative frequency of V gene usage") | |
168 write.table(x=PRODFV, file="VFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) | |
169 | |
170 png("VPlot.png",width = 1280, height = 720) | |
171 pV | |
172 dev.off(); | |
173 | |
174 pD = ggplot(PRODFD) | |
175 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)) | |
176 pD = pD + xlab("Summary of D gene") + ylab("Frequency") + ggtitle("Relative frequency of D gene usage") | |
177 write.table(x=PRODFD, file="DFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) | |
178 | |
179 png("DPlot.png",width = 800, height = 600) | |
180 pD | |
181 dev.off(); | |
182 | |
183 pJ = ggplot(PRODFJ) | |
184 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)) | |
185 pJ = pJ + xlab("Summary of J gene") + ylab("Frequency") + ggtitle("Relative frequency of J gene usage") | |
186 write.table(x=PRODFJ, file="JFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) | |
187 | |
188 png("JPlot.png",width = 800, height = 600) | |
189 pJ | |
190 dev.off(); | |
191 | |
192 print("test5\n") | |
193 | |
194 VGenes = PRODF[,c("Sample", "Top.V.Gene", "freq")] | |
195 VGenes$Top.V.Gene = gsub("-.*", "", VGenes$Top.V.Gene) | |
196 VGenes = data.frame(data.table(VGenes)[, list(Count=sum(freq)), by=c("Sample", "Top.V.Gene")]) | |
197 TotalPerSample = data.frame(data.table(VGenes)[, list(total=sum(.SD$Count)), by=Sample]) | |
198 VGenes = merge(VGenes, TotalPerSample, by="Sample") | |
199 VGenes$Frequency = VGenes$Count * 100 / VGenes$total | |
200 VPlot = ggplot(VGenes) | |
201 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)) + | |
202 ggtitle("Distribution of V gene families") + | |
203 ylab("Percentage of sequences") | |
204 png("VFPlot.png") | |
205 VPlot | |
206 dev.off(); | |
207 write.table(x=VGenes, file="VFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) | |
208 | |
209 DGenes = PRODF[,c("Sample", "Top.D.Gene", "freq")] | |
210 DGenes$Top.D.Gene = gsub("-.*", "", DGenes$Top.D.Gene) | |
211 DGenes = data.frame(data.table(DGenes)[, list(Count=sum(freq)), by=c("Sample", "Top.D.Gene")]) | |
212 TotalPerSample = data.frame(data.table(DGenes)[, list(total=sum(.SD$Count)), by=Sample]) | |
213 DGenes = merge(DGenes, TotalPerSample, by="Sample") | |
214 DGenes$Frequency = DGenes$Count * 100 / DGenes$total | |
215 DPlot = ggplot(DGenes) | |
216 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)) + | |
217 ggtitle("Distribution of D gene families") + | |
218 ylab("Percentage of sequences") | |
219 png("DFPlot.png") | |
220 DPlot | |
221 dev.off(); | |
222 write.table(x=DGenes, file="DFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) | |
223 | |
224 JGenes = PRODF[,c("Sample", "Top.J.Gene", "freq")] | |
225 JGenes$Top.J.Gene = gsub("-.*", "", JGenes$Top.J.Gene) | |
226 JGenes = data.frame(data.table(JGenes)[, list(Count=sum(freq)), by=c("Sample", "Top.J.Gene")]) | |
227 TotalPerSample = data.frame(data.table(JGenes)[, list(total=sum(.SD$Count)), by=Sample]) | |
228 JGenes = merge(JGenes, TotalPerSample, by="Sample") | |
229 JGenes$Frequency = JGenes$Count * 100 / JGenes$total | |
230 JPlot = ggplot(JGenes) | |
231 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)) + | |
232 ggtitle("Distribution of J gene families") + | |
233 ylab("Percentage of sequences") | |
234 png("JFPlot.png") | |
235 JPlot | |
236 dev.off(); | |
237 write.table(x=JGenes, file="JFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T) | |
238 | |
239 CDR3Length = data.frame(data.table(PRODF)[, list(Count=sum(freq)), by=c("Sample", "CDR3.Length.DNA")]) | |
240 TotalPerSample = data.frame(data.table(CDR3Length)[, list(total=sum(.SD$Count)), by=Sample]) | |
241 CDR3Length = merge(CDR3Length, TotalPerSample, by="Sample") | |
242 CDR3Length$Frequency = CDR3Length$Count * 100 / CDR3Length$total | |
243 CDR3LengthPlot = ggplot(CDR3Length) | |
244 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)) + | |
245 ggtitle("Length distribution of CDR3") + | |
246 xlab("CDR3 Length") + | |
247 ylab("Percentage of sequences") | |
248 png("CDR3LengthPlot.png",width = 1280, height = 720) | |
249 CDR3LengthPlot | |
250 dev.off() | |
251 write.table(x=CDR3Length, file="CDR3LengthPlot.csv", sep=",",quote=F,row.names=F,col.names=T) | |
252 | |
253 revVchain = Vchain | |
254 revDchain = Dchain | |
255 revVchain$chr.orderV = rev(revVchain$chr.orderV) | |
256 revDchain$chr.orderD = rev(revDchain$chr.orderD) | |
257 | |
258 print("test6\n") | |
259 | |
260 plotVD <- function(dat){ | |
261 if(length(dat[,1]) == 0){ | |
262 return() | |
263 } | |
264 img = ggplot() + | |
265 geom_tile(data=dat, aes(x=factor(reorder(Top.D.Gene, chr.orderD)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) + | |
266 theme(axis.text.x = element_text(angle = 90, hjust = 1)) + | |
267 scale_fill_gradient(low="gold", high="blue", na.value="white") + | |
268 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + | |
269 xlab("D genes") + | |
270 ylab("V Genes") | |
271 | |
272 png(paste("HeatmapVD_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Dchain$v.name)), height=100+(15*length(Vchain$v.name))) | |
273 print(img) | |
274 | |
275 dev.off() | |
276 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) | |
277 } | |
278 | |
279 VandDCount = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Top.V.Gene", "Top.D.Gene", "Sample")]) | |
280 | |
281 VandDCount$l = log(VandDCount$Length) | |
282 maxVD = data.frame(data.table(VandDCount)[, list(max=max(l)), by=c("Sample")]) | |
283 VandDCount = merge(VandDCount, maxVD, by.x="Sample", by.y="Sample", all.x=T) | |
284 VandDCount$relLength = VandDCount$l / VandDCount$max | |
285 | |
286 cartegianProductVD = expand.grid(Top.V.Gene = Vchain$v.name, Top.D.Gene = Dchain$v.name, Sample = unique(test$Sample)) | |
287 | |
288 completeVD = merge(VandDCount, cartegianProductVD, all.y=TRUE) | |
289 completeVD = merge(completeVD, revVchain, by.x="Top.V.Gene", by.y="v.name", all.x=TRUE) | |
290 completeVD = merge(completeVD, Dchain, by.x="Top.D.Gene", by.y="v.name", all.x=TRUE) | |
291 VDList = split(completeVD, f=completeVD[,"Sample"]) | |
292 | |
293 lapply(VDList, FUN=plotVD) | |
294 | |
295 plotVJ <- function(dat){ | |
296 if(length(dat[,1]) == 0){ | |
297 return() | |
298 } | |
299 img = ggplot() + | |
300 geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) + | |
301 theme(axis.text.x = element_text(angle = 90, hjust = 1)) + | |
302 scale_fill_gradient(low="gold", high="blue", na.value="white") + | |
303 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + | |
304 xlab("J genes") + | |
305 ylab("V Genes") | |
306 | |
307 png(paste("HeatmapVJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Vchain$v.name))) | |
308 print(img) | |
309 dev.off() | |
310 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) | |
311 } | |
312 | |
313 VandJCount = data.frame(data.table(PRODF)[, list(Length=sum(freq)), by=c("Top.V.Gene", "Top.J.Gene", "Sample")]) | |
314 | |
315 VandJCount$l = log(VandJCount$Length) | |
316 maxVJ = data.frame(data.table(VandJCount)[, list(max=max(l)), by=c("Sample")]) | |
317 VandJCount = merge(VandJCount, maxVJ, by.x="Sample", by.y="Sample", all.x=T) | |
318 VandJCount$relLength = VandJCount$l / VandJCount$max | |
319 | |
320 cartegianProductVJ = expand.grid(Top.V.Gene = Vchain$v.name, Top.J.Gene = Jchain$v.name, Sample = unique(test$Sample)) | |
321 | |
322 completeVJ = merge(VandJCount, cartegianProductVJ, all.y=TRUE) | |
323 completeVJ = merge(completeVJ, revVchain, by.x="Top.V.Gene", by.y="v.name", all.x=TRUE) | |
324 completeVJ = merge(completeVJ, Jchain, by.x="Top.J.Gene", by.y="v.name", all.x=TRUE) | |
325 VJList = split(completeVJ, f=completeVJ[,"Sample"]) | |
326 lapply(VJList, FUN=plotVJ) | |
327 | |
328 plotDJ <- function(dat){ | |
329 if(length(dat[,1]) == 0){ | |
330 return() | |
331 } | |
332 img = ggplot() + | |
333 geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.D.Gene, chr.orderD)), fill=relLength)) + | |
334 theme(axis.text.x = element_text(angle = 90, hjust = 1)) + | |
335 scale_fill_gradient(low="gold", high="blue", na.value="white") + | |
336 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + | |
337 xlab("J genes") + | |
338 ylab("D Genes") | |
339 | |
340 png(paste("HeatmapDJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Dchain$v.name))) | |
341 print(img) | |
342 dev.off() | |
343 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) | |
344 } | |
345 | |
346 DandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.D.Gene", "Top.J.Gene", "Sample")]) | |
347 | |
348 DandJCount$l = log(DandJCount$Length) | |
349 maxDJ = data.frame(data.table(DandJCount)[, list(max=max(l)), by=c("Sample")]) | |
350 DandJCount = merge(DandJCount, maxDJ, by.x="Sample", by.y="Sample", all.x=T) | |
351 DandJCount$relLength = DandJCount$l / DandJCount$max | |
352 | |
353 cartegianProductDJ = expand.grid(Top.D.Gene = Dchain$v.name, Top.J.Gene = Jchain$v.name, Sample = unique(test$Sample)) | |
354 | |
355 completeDJ = merge(DandJCount, cartegianProductDJ, all.y=TRUE) | |
356 completeDJ = merge(completeDJ, revDchain, by.x="Top.D.Gene", by.y="v.name", all.x=TRUE) | |
357 completeDJ = merge(completeDJ, Jchain, by.x="Top.J.Gene", by.y="v.name", all.x=TRUE) | |
358 DJList = split(completeDJ, f=completeDJ[,"Sample"]) | |
359 lapply(DJList, FUN=plotDJ) | |
360 | |
361 sampleFile <- file("samples.txt") | |
362 un = unique(test$Sample) | |
363 un = paste(un, sep="\n") | |
364 writeLines(un, sampleFile) | |
365 close(sampleFile) | |
366 | |
367 print("test7\n") | |
368 | |
369 if("Replicate" %in% colnames(test)) | |
370 { | |
371 clonalityFrame = PROD | |
372 clonalityFrame$ReplicateConcat = do.call(paste, c(clonalityFrame[c("VDJCDR3", "Sample", "Replicate")], sep = ":")) | |
373 clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ] | |
374 write.table(clonalityFrame, "clonalityComplete.csv", sep=",",quote=F,row.names=F,col.names=T) | |
375 | |
376 ClonalitySampleReplicatePrint <- function(dat){ | |
377 write.table(dat, paste("clonality_", unique(dat$Sample) , "_", unique(dat$Replicate), ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=T) | |
378 } | |
379 | |
380 clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,c("Sample", "Replicate")]) | |
381 #lapply(clonalityFrameSplit, FUN=ClonalitySampleReplicatePrint) | |
382 | |
383 ClonalitySamplePrint <- function(dat){ | |
384 write.table(dat, paste("clonality_", unique(dat$Sample) , ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=T) | |
385 } | |
386 | |
387 clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,"Sample"]) | |
388 #lapply(clonalityFrameSplit, FUN=ClonalitySamplePrint) | |
389 | |
390 clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "VDJCDR3")]) | |
391 clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")]) | |
392 clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count | |
393 clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")]) | |
394 clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample") | |
395 | |
396 ct = c('Type\tWeight\n2\t1\n3\t3\n4\t6\n5\t10\n6\t15') | |
397 tcct = textConnection(ct) | |
398 CT = read.table(tcct, sep="\t", header=TRUE) | |
399 close(tcct) | |
400 clonalFreqCount = merge(clonalFreqCount, CT, by.x="Type", by.y="Type", all.x=T) | |
401 clonalFreqCount$WeightedCount = clonalFreqCount$Count * clonalFreqCount$Weight | |
402 | |
403 ReplicateReads = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "Replicate", "VDJCDR3")]) | |
404 ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")]) | |
405 clonalFreqCount$Reads = as.numeric(clonalFreqCount$Reads) | |
406 ReplicateReads$squared = ReplicateReads$Reads * ReplicateReads$Reads | |
407 | |
408 ReplicatePrint <- function(dat){ | |
409 write.table(dat[-1], paste("ReplicateReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) | |
410 } | |
411 | |
412 ReplicateSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"]) | |
413 lapply(ReplicateSplit, FUN=ReplicatePrint) | |
414 | |
415 ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(Reads), ReadsSquaredSum=sum(squared)), by=c("Sample")]) | |
416 clonalFreqCount = merge(clonalFreqCount, ReplicateReads, by.x="Sample", by.y="Sample", all.x=T) | |
417 | |
418 | |
419 ReplicateSumPrint <- function(dat){ | |
420 write.table(dat[-1], paste("ReplicateSumReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) | |
421 } | |
422 | |
423 ReplicateSumSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"]) | |
424 lapply(ReplicateSumSplit, FUN=ReplicateSumPrint) | |
425 | |
426 writeClonalitySequences <- function(dat){ | |
427 for(i in c(2,3,4,5,6)){ | |
428 fltr = dat[dat$Type == i,] | |
429 if(length(fltr[,1]) == 0){ | |
430 next | |
431 } | |
432 tmp = clonalityFrame[clonalityFrame$Sample == fltr$Sample[1] & clonalityFrame$VDJCDR3 %in% fltr$VDJCDR3,] | |
433 tmp = tmp[order(tmp$VDJCDR3),] | |
434 write.table(tmp, paste("ClonalitySequences_", unique(dat[1])[1,1] , "_", i, ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=T) | |
435 } | |
436 } | |
437 freqsplt = split(clonalFreq[clonalFreq$Type > 1,], clonalFreq[clonalFreq$Type > 1,]$Sample) | |
438 lapply(freqsplt, FUN=writeClonalitySequences) | |
439 | |
440 clonalFreqCountSum = data.frame(data.table(clonalFreqCount)[, list(Numerator=sum(WeightedCount, na.rm=T)), by=c("Sample")]) | |
441 clonalFreqCount = merge(clonalFreqCount, clonalFreqCountSum, by.x="Sample", by.y="Sample", all.x=T) | |
442 clonalFreqCount$ReadsSum = as.numeric(clonalFreqCount$ReadsSum) #prevent integer overflow | |
443 clonalFreqCount$Denominator = (((clonalFreqCount$ReadsSum * clonalFreqCount$ReadsSum) - clonalFreqCount$ReadsSquaredSum) / 2) | |
444 clonalFreqCount$Result = (clonalFreqCount$Numerator + 1) / (clonalFreqCount$Denominator + 1) | |
445 | |
446 ClonalityScorePrint <- function(dat){ | |
447 write.table(dat$Result, paste("ClonalityScore_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) | |
448 } | |
449 | |
450 clonalityScore = clonalFreqCount[c("Sample", "Result")] | |
451 clonalityScore = unique(clonalityScore) | |
452 | |
453 clonalityScoreSplit = split(clonalityScore, f=clonalityScore[,"Sample"]) | |
454 lapply(clonalityScoreSplit, FUN=ClonalityScorePrint) | |
455 | |
456 clonalityOverview = clonalFreqCount[c("Sample", "Type", "Count", "Weight", "WeightedCount")] | |
457 | |
458 | |
459 | |
460 ClonalityOverviewPrint <- function(dat){ | |
461 write.table(dat[-1], paste("ClonalityOverView_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F) | |
462 } | |
463 | |
464 clonalityOverviewSplit = split(clonalityOverview, f=clonalityOverview$Sample) | |
465 lapply(clonalityOverviewSplit, FUN=ClonalityOverviewPrint) | |
466 } | |
467 | |
468 print("test8\n") | |
469 | |
470 if("Functionality" %in% colnames(test)) | |
471 { | |
472 newData = data.frame(data.table(PROD)[,list(unique=.N, | |
473 VH.DEL=mean(X3V.REGION.trimmed.nt.nb), | |
474 P1=mean(P3V.nt.nb), | |
475 N1=mean(N1.REGION.nt.nb), | |
476 P2=mean(P5D.nt.nb), | |
477 DEL.DH=mean(X5D.REGION.trimmed.nt.nb), | |
478 DH.DEL=mean(X3D.REGION.trimmed.nt.nb), | |
479 P3=mean(P3D.nt.nb), | |
480 N2=mean(N2.REGION.nt.nb), | |
481 P4=mean(P5J.nt.nb), | |
482 DEL.JH=mean(X5J.REGION.trimmed.nt.nb), | |
483 Total.Del=( mean(X3V.REGION.trimmed.nt.nb) + | |
484 mean(X5D.REGION.trimmed.nt.nb) + | |
485 mean(X3D.REGION.trimmed.nt.nb) + | |
486 mean(X5J.REGION.trimmed.nt.nb)), | |
487 | |
488 Total.N=( mean(N1.REGION.nt.nb) + | |
489 mean(N2.REGION.nt.nb)), | |
490 | |
491 Total.P=( mean(P3V.nt.nb) + | |
492 mean(P5D.nt.nb) + | |
493 mean(P3D.nt.nb) + | |
494 mean(P5J.nt.nb))), | |
495 by=c("Sample")]) | |
496 write.table(newData, "junctionAnalysis.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F) | |
497 } | |
498 | |
499 print("test9\n") |