Mercurial > repos > davidvanzessen > plotting_merged

comparison RScript.r @ 49:3da85f1a7f42 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Mon, 09 Dec 2013 05:42:56 -0500
parents 33ebd6f0d256
children 8f5ab5290c88
comparison
equal deleted inserted replaced
48:dea2d3353a42 49:3da85f1a7f42
87 Jchain = read.table(tcJ, sep="\t", header=TRUE) 87 Jchain = read.table(tcJ, sep="\t", header=TRUE)
88 PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE) 88 PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE)
89 close(tcJ) 89 close(tcJ)
90 90
91 setwd(outDir) 91 setwd(outDir)
92
93 write.table(PRODF, "allUnique.tsv", sep="\t",quote=F,row.names=T,col.names=T)
92 94
93 pV = ggplot(PRODFV) 95 pV = ggplot(PRODFV)
94 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)) 96 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))
95 pV = pV + xlab("Summary of V gene") + ylab("Frequency") + ggtitle("Relative frequency of V gene usage") 97 pV = pV + xlab("Summary of V gene") + ylab("Frequency") + ggtitle("Relative frequency of V gene usage")
96 98
229 if("Replicate" %in% colnames(test)) 231 if("Replicate" %in% colnames(test))
230 { 232 {
231 clonalityFrame = PROD 233 clonalityFrame = PROD
232 clonalityFrame$ReplicateConcat = do.call(paste, c(clonalityFrame[c("VDJCDR3", "Sample", "Replicate")], sep = ":")) 234 clonalityFrame$ReplicateConcat = do.call(paste, c(clonalityFrame[c("VDJCDR3", "Sample", "Replicate")], sep = ":"))
233 clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ] 235 clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ]
236 write.table(clonalityFrame, "clonalityComplete.tsv", sep="\t",quote=F,row.names=T,col.names=T)
237
238 ClonalitySampleReplicatePrint <- function(dat){
239 write.table(dat, paste("clonality_", unique(dat$Sample) , "_", unique(dat$Replicate), ".tsv", sep=""), sep="\t",quote=F,row.names=T,col.names=T)
240 }
241
242 clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,c("Sample", "Replicate")])
243 lapply(clonalityFrameSplit, FUN=ClonalitySampleReplicatePrint)
244
245 ClonalitySamplePrint <- function(dat){
246 write.table(dat, paste("clonality_", unique(dat$Sample) , ".tsv", sep=""), sep="\t",quote=F,row.names=T,col.names=T)
247 }
248
249 clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,"Sample"])
250 lapply(clonalityFrameSplit, FUN=ClonalitySamplePrint)
251
234 clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "VDJCDR3")]) 252 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")]) 253 clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")])
236 clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count 254 clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count
237 clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")]) 255 clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")])
238 clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample") 256 clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample")