Mercurial > repos > davidvanzessen > plotting_merged

comparison RScript.r @ 47:33ebd6f0d256 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Mon, 02 Dec 2013 10:56:40 -0500
parents a2db27f58627
children 3da85f1a7f42
comparison
equal deleted inserted replaced
46:a2db27f58627 47:33ebd6f0d256
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]
27 27
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
33 if("Replicate" %in% colnames(test))
34 {
35 test$SRID = do.call(paste, c(test[c("Sample", "Replicate")], sep = "-"))
36 }
37 32
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
72 Total = 0 67 Total = 0
73 Total = ddply(PRODFJ, .(Sample), function(x) data.frame(Total = sum(x$Length))) 68 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) 69 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)) 70 PRODFJ = ddply(PRODFJ, c("Sample", "Top.J.Gene"), summarise, relFreq= (Length*100 / Total))
76 71
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") 72 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) 73 tcV = textConnection(V)
79 Vchain = read.table(tcV, sep="\t", header=TRUE) 74 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) 75 PRODFV = merge(PRODFV, Vchain, by.x='Top.V.Gene', by.y='v.name', all.x=TRUE)
81 close(tcV) 76 close(tcV)
82 77
129 return() 124 return()
130 } 125 }
131 img = ggplot() + 126 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)) + 127 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)) + 128 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)) + 129 scale_fill_gradient(low="gold", high="blue", na.value="white") +
135 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 130 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) +
136 xlab("D genes") + 131 xlab("D genes") +
137 ylab("V Genes") 132 ylab("V Genes")
138 133
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))) 134 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() 159 return()
165 } 160 }
166 img = ggplot() + 161 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)) + 162 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)) + 163 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)) + 164 scale_fill_gradient(low="gold", high="blue", na.value="white") +
170 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 165 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) +
171 xlab("J genes") + 166 xlab("J genes") +
172 ylab("V Genes") 167 ylab("V Genes")
173 168
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))) 169 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() 191 return()
197 } 192 }
198 img = ggplot() + 193 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)) + 194 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)) + 195 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)) + 196 scale_fill_gradient(low="gold", high="blue", na.value="white") +
202 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 197 ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) +
203 xlab("J genes") + 198 xlab("J genes") +
204 ylab("D Genes") 199 ylab("D Genes")
205 200
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))) 201 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") 222 sampleFile <- file("samples.txt")
228 un = unique(test$Sample) 223 un = unique(test$Sample)
229 un = paste(un, sep="\n") 224 un = paste(un, sep="\n")
230 writeLines(un, sampleFile) 225 writeLines(un, sampleFile)
231 close(sampleFile) 226 close(sampleFile)
227
228
229 if("Replicate" %in% colnames(test))
230 {
231 clonalityFrame = PROD
232 clonalityFrame$ReplicateConcat = do.call(paste, c(clonalityFrame[c("VDJCDR3", "Sample", "Replicate")], sep = ":"))
233 clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ]
234 clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "VDJCDR3")])
235 clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")])
236 clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count
237 clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")])
238 clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample")
239
240 ct = c('Type\tWeight\n2\t1\n3\t3\n4\t6\n5\t10\n6\t15')
241 tcct = textConnection(ct)
242 CT = read.table(tcct, sep="\t", header=TRUE)
243 close(tcct)
244 clonalFreqCount = merge(clonalFreqCount, CT, by.x="Type", by.y="Type", all.x=T)
245 clonalFreqCount$WeightedCount = clonalFreqCount$Count * clonalFreqCount$Weight
246
247 ReplicateReads = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "Replicate", "VDJCDR3")])
248 ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")])
249 ReplicateReads$squared = ReplicateReads$Reads * ReplicateReads$Reads
250
251 ReplicatePrint <- function(dat){
252 write.table(dat[-1], paste("ReplicateReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
253 }
254
255 ReplicateSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"])
256 lapply(ReplicateSplit, FUN=ReplicatePrint)
257
258 ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(Reads), ReadsSquaredSum=sum(squared)), by=c("Sample")])
259 clonalFreqCount = merge(clonalFreqCount, ReplicateReads, by.x="Sample", by.y="Sample", all.x=T)
260
261
262 ReplicateSumPrint <- function(dat){
263 write.table(dat[-1], paste("ReplicateSumReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
264 }
265
266 ReplicateSumSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"])
267 lapply(ReplicateSumSplit, FUN=ReplicateSumPrint)
268
269 clonalFreqCountSum = data.frame(data.table(clonalFreqCount)[, list(Numerator=sum(WeightedCount, na.rm=T)), by=c("Sample")])
270 clonalFreqCount = merge(clonalFreqCount, clonalFreqCountSum, by.x="Sample", by.y="Sample", all.x=T)
271
272 clonalFreqCount$Denominator = (((clonalFreqCount$ReadsSum * clonalFreqCount$ReadsSum) - clonalFreqCount$ReadsSquaredSum) / 2)
273 clonalFreqCount$Result = (clonalFreqCount$Numerator + 1) / (clonalFreqCount$Denominator + 1)
274
275 ClonalityScorePrint <- function(dat){
276 write.table(dat$Result, paste("ClonalityScore_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
277 }
278
279 clonalityScore = clonalFreqCount[c("Sample", "Result")]
280 clonalityScore = unique(clonalityScore)
281
282 clonalityScoreSplit = split(clonalityScore, f=clonalityScore[,"Sample"])
283 lapply(clonalityScoreSplit, FUN=ClonalityScorePrint)
284
285 clonalityOverview = clonalFreqCount[c("Sample", "Type", "Count", "Weight", "WeightedCount")]
286
287
288
289 ClonalityOverviewPrint <- function(dat){
290 write.table(dat[-1], paste("ClonalityOverView_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
291 }
292
293 clonalityOverviewSplit = split(clonalityOverview, f=clonalityOverview$Sample)
294 lapply(clonalityOverviewSplit, FUN=ClonalityOverviewPrint)
295 }