Mercurial > repos > davidvanzessen > combined_immune_repertoire_pipeline

comparison RScript.r @ 9:8d83319a0f3d draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Tue, 10 Dec 2013 05:53:08 -0500
parents 00d432c66fb8
children a5c224bb0be5
comparison
equal deleted inserted replaced
8:00d432c66fb8 9:8d83319a0f3d
1 options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } ) 1 #options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
2 2
3 args <- commandArgs(trailingOnly = TRUE) 3 args <- commandArgs(trailingOnly = TRUE)
4 4
5 inFile = args[1] 5 inFile = args[1]
6 outFile = args[2] 6 outFile = args[2]
28 28
29 test = read.table(inFile, sep="\t", header=TRUE, fill=T) 29 test = read.table(inFile, sep="\t", header=TRUE, fill=T)
30 30
31 test = test[test$Sample != "",] 31 test = test[test$Sample != "",]
32 32
33 if("Replicate" %in% colnames(test))
34 {
35 test$SRID = do.call(paste, c(test[c("Sample", "Replicate")], sep = "-"))
36 }
37
38 test$Top.V.Gene = gsub("[*]([0-9]+)", "", test$Top.V.Gene) 33 test$Top.V.Gene = gsub("[*]([0-9]+)", "", test$Top.V.Gene)
39 test$Top.D.Gene = gsub("[*]([0-9]+)", "", test$Top.D.Gene) 34 test$Top.D.Gene = gsub("[*]([0-9]+)", "", test$Top.D.Gene)
40 test$Top.J.Gene = gsub("[*]([0-9]+)", "", test$Top.J.Gene) 35 test$Top.J.Gene = gsub("[*]([0-9]+)", "", test$Top.J.Gene)
41 36
42 #test$VDJCDR3 = do.call(paste, c(test[c("Top.V.Gene", "Top.D.Gene", "Top.J.Gene","CDR3.Seq.DNA")], sep = ":")) 37 #test$VDJCDR3 = do.call(paste, c(test[c("Top.V.Gene", "Top.D.Gene", "Top.J.Gene","CDR3.Seq.DNA")], sep = ":"))
43 test$VDJCDR3 = do.call(paste, c(test[unlist(strsplit(clonalType, ","))], sep = ":")) 38 test$VDJCDR3 = do.call(paste, c(test[unlist(strsplit(clonalType, ","))], sep = ":"))
44 39
45 PROD = test[test$VDJ.Frame != "In-frame with stop codon" & test$VDJ.Frame != "Out-of-frame" & test$CDR3.Found.How != "NOT_FOUND" , ] 40 PROD = test[test$VDJ.Frame != "In-frame with stop codon" & test$VDJ.Frame != "Out-of-frame" & test$CDR3.Found.How != "NOT_FOUND" , ]
41 if("Functionality" %in% colnames(test)) {
42 PROD = test[test$Functionality == "productive" | test$Functionality == "productive (see comment)", ]
43 }
46 44
47 NONPROD = test[test$VDJ.Frame == "In-frame with stop codon" | test$VDJ.Frame == "Out-of-frame" | test$CDR3.Found.How == "NOT_FOUND" , ] 45 NONPROD = test[test$VDJ.Frame == "In-frame with stop codon" | test$VDJ.Frame == "Out-of-frame" | test$CDR3.Found.How == "NOT_FOUND" , ]
48 46
49 #PRODF = PROD[ -1] 47 #PRODF = PROD[ -1]
50 48
72 Total = 0 70 Total = 0
73 Total = ddply(PRODFJ, .(Sample), function(x) data.frame(Total = sum(x$Length))) 71 Total = ddply(PRODFJ, .(Sample), function(x) data.frame(Total = sum(x$Length)))
74 PRODFJ = merge(PRODFJ, Total, by.x='Sample', by.y='Sample', all.x=TRUE) 72 PRODFJ = merge(PRODFJ, Total, by.x='Sample', by.y='Sample', all.x=TRUE)
75 PRODFJ = ddply(PRODFJ, c("Sample", "Top.J.Gene"), summarise, relFreq= (Length*100 / Total)) 73 PRODFJ = ddply(PRODFJ, c("Sample", "Top.J.Gene"), summarise, relFreq= (Length*100 / Total))
76 74
77 V = c("v.name\tchr.orderV\nIGHV7-81\t1\nIGHV3-74\t2\nIGHV3-73\t3\nIGHV3-72\t4\nIGHV3-71\t5\nIGHV2-70\t6\nIGHV1-69\t7\nIGHV3-66\t8\nIGHV3-64\t9\nIGHV4-61\t10\nIGHV4-59\t11\nIGHV1-58\t12\nIGHV3-53\t13\nIGHV3-52\t14\nIGHV5-a\t15\nIGHV5-51\t16\nIGHV3-49\t17\nIGHV3-48\t18\nIGHV3-47\t19\nIGHV1-46\t20\nIGHV1-45\t21\nIGHV3-43\t22\nIGHV4-39\t23\nIGHV3-35\t24\nIGHV4-34\t25\nIGHV3-33\t26\nIGHV4-31\t27\nIGHV4-30-4\t28\nIGHV4-30-2\t29\nIGHV3-30-3\t30\nIGHV3-30\t31\nIGHV4-28\t32\nIGHV2-26\t33\nIGHV1-24\t34\nIGHV3-23\t35\nIGHV3-22\t36\nIGHV3-21\t37\nIGHV3-20\t38\nIGHV3-19\t39\nIGHV1-18\t40\nIGHV3-15\t41\nIGHV3-13\t42\nIGHV3-11\t43\nIGHV3-9\t44\nIGHV1-8\t45\nIGHV3-7\t46\nIGHV2-5\t47\nIGHV7-4-1\t48\nIGHV4-4\t49\nIGHV4-b\t50\nIGHV1-3\t51\nIGHV1-2\t52\nIGHV6-1\t53") 75 V = c("v.name\tchr.orderV\nIGHV7-81\t1\nIGHV3-74\t2\nIGHV3-73\t3\nIGHV3-72\t4\nIGHV2-70\t6\nIGHV1-69\t7\nIGHV3-66\t8\nIGHV3-64\t9\nIGHV4-61\t10\nIGHV4-59\t11\nIGHV1-58\t12\nIGHV3-53\t13\nIGHV5-a\t15\nIGHV5-51\t16\nIGHV3-49\t17\nIGHV3-48\t18\nIGHV1-46\t20\nIGHV1-45\t21\nIGHV3-43\t22\nIGHV4-39\t23\nIGHV3-35\t24\nIGHV4-34\t25\nIGHV3-33\t26\nIGHV4-31\t27\nIGHV4-30-4\t28\nIGHV4-30-2\t29\nIGHV3-30-3\t30\nIGHV3-30\t31\nIGHV4-28\t32\nIGHV2-26\t33\nIGHV1-24\t34\nIGHV3-23\t35\nIGHV3-21\t37\nIGHV3-20\t38\nIGHV1-18\t40\nIGHV3-15\t41\nIGHV3-13\t42\nIGHV3-11\t43\nIGHV3-9\t44\nIGHV1-8\t45\nIGHV3-7\t46\nIGHV2-5\t47\nIGHV7-4-1\t48\nIGHV4-4\t49\nIGHV4-b\t50\nIGHV1-3\t51\nIGHV1-2\t52\nIGHV6-1\t53")
78 tcV = textConnection(V) 76 tcV = textConnection(V)
79 Vchain = read.table(tcV, sep="\t", header=TRUE) 77 Vchain = read.table(tcV, sep="\t", header=TRUE)
80 PRODFV = merge(PRODFV, Vchain, by.x='Top.V.Gene', by.y='v.name', all.x=TRUE) 78 PRODFV = merge(PRODFV, Vchain, by.x='Top.V.Gene', by.y='v.name', all.x=TRUE)
81 close(tcV) 79 close(tcV)
82 80
93 PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE) 91 PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE)
94 close(tcJ) 92 close(tcJ)
95 93
96 setwd(outDir) 94 setwd(outDir)
97 95
96 write.table(PRODF, "allUnique.tsv", sep="\t",quote=F,row.names=T,col.names=T)
97
98 pV = ggplot(PRODFV) 98 pV = ggplot(PRODFV)
99 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)) 99 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))
100 pV = pV + xlab("Summary of V gene") + ylab("Frequency") + ggtitle("Relative frequency of V gene usage") 100 pV = pV + xlab("Summary of V gene") + ylab("Frequency") + ggtitle("Relative frequency of V gene usage")
101 101
102 png("VPlot.png",width = 1280, height = 720) 102 png("VPlot.png",width = 1280, height = 720)
129 return() 129 return()
130 } 130 }
131 img = ggplot() + 131 img = ggplot() +
132 geom_tile(data=dat, aes(x=factor(reorder(Top.D.Gene, chr.orderD)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) + 132 geom_tile(data=dat, aes(x=factor(reorder(Top.D.Gene, chr.orderD)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) +
133 theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 133 theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
134 scale_fill_gradient(low="gold", high="blue", na.value="white", limits=c(0,1)) + 134 scale_fill_gradient(low="gold", high="blue", na.value="white") +
135 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 135 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) +
136 xlab("D genes") + 136 xlab("D genes") +
137 ylab("V Genes") 137 ylab("V Genes")
138 138
139 png(paste("HeatmapVD_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Dchain$v.name)), height=100+(15*length(Vchain$v.name))) 139 png(paste("HeatmapVD_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Dchain$v.name)), height=100+(15*length(Vchain$v.name)))
164 return() 164 return()
165 } 165 }
166 img = ggplot() + 166 img = ggplot() +
167 geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) + 167 geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) +
168 theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 168 theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
169 scale_fill_gradient(low="gold", high="blue", na.value="white", limits=c(0,1)) + 169 scale_fill_gradient(low="gold", high="blue", na.value="white") +
170 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 170 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) +
171 xlab("J genes") + 171 xlab("J genes") +
172 ylab("V Genes") 172 ylab("V Genes")
173 173
174 png(paste("HeatmapVJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Vchain$v.name))) 174 png(paste("HeatmapVJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Vchain$v.name)))
196 return() 196 return()
197 } 197 }
198 img = ggplot() + 198 img = ggplot() +
199 geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.D.Gene, chr.orderD)), fill=relLength)) + 199 geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.D.Gene, chr.orderD)), fill=relLength)) +
200 theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 200 theme(axis.text.x = element_text(angle = 90, hjust = 1)) +
201 scale_fill_gradient(low="gold", high="blue", na.value="white", limits=c(0,1)) + 201 scale_fill_gradient(low="gold", high="blue", na.value="white") +
202 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 202 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) +
203 xlab("J genes") + 203 xlab("J genes") +
204 ylab("D Genes") 204 ylab("D Genes")
205 205
206 png(paste("HeatmapDJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Dchain$v.name))) 206 png(paste("HeatmapDJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Dchain$v.name)))
227 sampleFile <- file("samples.txt") 227 sampleFile <- file("samples.txt")
228 un = unique(test$Sample) 228 un = unique(test$Sample)
229 un = paste(un, sep="\n") 229 un = paste(un, sep="\n")
230 writeLines(un, sampleFile) 230 writeLines(un, sampleFile)
231 close(sampleFile) 231 close(sampleFile)
232
233
234 if("Replicate" %in% colnames(test))
235 {
236 clonalityFrame = PROD
237 clonalityFrame$ReplicateConcat = do.call(paste, c(clonalityFrame[c("VDJCDR3", "Sample", "Replicate")], sep = ":"))
238 clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ]
239 write.table(clonalityFrame, "clonalityComplete.tsv", sep="\t",quote=F,row.names=T,col.names=T)
240
241 ClonalitySampleReplicatePrint <- function(dat){
242 write.table(dat, paste("clonality_", unique(dat$Sample) , "_", unique(dat$Replicate), ".tsv", sep=""), sep="\t",quote=F,row.names=T,col.names=T)
243 }
244
245 clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,c("Sample", "Replicate")])
246 lapply(clonalityFrameSplit, FUN=ClonalitySampleReplicatePrint)
247
248 ClonalitySamplePrint <- function(dat){
249 write.table(dat, paste("clonality_", unique(dat$Sample) , ".tsv", sep=""), sep="\t",quote=F,row.names=T,col.names=T)
250 }
251
252 clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,"Sample"])
253 lapply(clonalityFrameSplit, FUN=ClonalitySamplePrint)
254
255 clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "VDJCDR3")])
256 clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")])
257 clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count
258 clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")])
259 clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample")
260
261 ct = c('Type\tWeight\n2\t1\n3\t3\n4\t6\n5\t10\n6\t15')
262 tcct = textConnection(ct)
263 CT = read.table(tcct, sep="\t", header=TRUE)
264 close(tcct)
265 clonalFreqCount = merge(clonalFreqCount, CT, by.x="Type", by.y="Type", all.x=T)
266 clonalFreqCount$WeightedCount = clonalFreqCount$Count * clonalFreqCount$Weight
267
268 ReplicateReads = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "Replicate", "VDJCDR3")])
269 ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")])
270 ReplicateReads$squared = ReplicateReads$Reads * ReplicateReads$Reads
271
272 ReplicatePrint <- function(dat){
273 write.table(dat[-1], paste("ReplicateReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
274 }
275
276 ReplicateSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"])
277 lapply(ReplicateSplit, FUN=ReplicatePrint)
278
279 ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(Reads), ReadsSquaredSum=sum(squared)), by=c("Sample")])
280 clonalFreqCount = merge(clonalFreqCount, ReplicateReads, by.x="Sample", by.y="Sample", all.x=T)
281
282
283 ReplicateSumPrint <- function(dat){
284 write.table(dat[-1], paste("ReplicateSumReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
285 }
286
287 ReplicateSumSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"])
288 lapply(ReplicateSumSplit, FUN=ReplicateSumPrint)
289
290 clonalFreqCountSum = data.frame(data.table(clonalFreqCount)[, list(Numerator=sum(WeightedCount, na.rm=T)), by=c("Sample")])
291 clonalFreqCount = merge(clonalFreqCount, clonalFreqCountSum, by.x="Sample", by.y="Sample", all.x=T)
292
293 clonalFreqCount$Denominator = (((clonalFreqCount$ReadsSum * clonalFreqCount$ReadsSum) - clonalFreqCount$ReadsSquaredSum) / 2)
294 clonalFreqCount$Result = (clonalFreqCount$Numerator + 1) / (clonalFreqCount$Denominator + 1)
295
296 ClonalityScorePrint <- function(dat){
297 write.table(dat$Result, paste("ClonalityScore_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
298 }
299
300 clonalityScore = clonalFreqCount[c("Sample", "Result")]
301 clonalityScore = unique(clonalityScore)
302
303 clonalityScoreSplit = split(clonalityScore, f=clonalityScore[,"Sample"])
304 lapply(clonalityScoreSplit, FUN=ClonalityScorePrint)
305
306 clonalityOverview = clonalFreqCount[c("Sample", "Type", "Count", "Weight", "WeightedCount")]
307
308
309
310 ClonalityOverviewPrint <- function(dat){
311 write.table(dat[-1], paste("ClonalityOverView_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
312 }
313
314 clonalityOverviewSplit = split(clonalityOverview, f=clonalityOverview$Sample)
315 lapply(clonalityOverviewSplit, FUN=ClonalityOverviewPrint)
316 }