Mercurial > repos > davidvanzessen > report_clonality_igg

comparison RScript.r @ 6:18ac92a69ef1 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Tue, 18 Mar 2014 09:03:54 -0400
parents 99834201251f
children e2972f0935e9
comparison
equal deleted inserted replaced
5:99834201251f 6:18ac92a69ef1
22 22
23 if (!("data.table" %in% rownames(installed.packages()))) { 23 if (!("data.table" %in% rownames(installed.packages()))) {
24 install.packages("data.table", repos="http://cran.xl-mirror.nl/") 24 install.packages("data.table", repos="http://cran.xl-mirror.nl/")
25 } 25 }
26 library(data.table) 26 library(data.table)
27
28 if (!("reshape2" %in% rownames(installed.packages()))) {
29 install.packages("reshape2", repos="http://cran.xl-mirror.nl/")
30 }
31 library(reshape2)
27 32
28 33
29 test = read.table(inFile, sep="\t", header=TRUE, fill=T, comment.char="") 34 test = read.table(inFile, sep="\t", header=TRUE, fill=T, comment.char="")
30 35
31 test = test[test$Sample != "",] 36 test = test[test$Sample != "",]
91 PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE) 96 PRODFJ = merge(PRODFJ, Jchain, by.x='Top.J.Gene', by.y='v.name', all.x=TRUE)
92 close(tcJ) 97 close(tcJ)
93 98
94 setwd(outDir) 99 setwd(outDir)
95 100
96 write.table(PRODF, "allUnique.tsv", sep="\t",quote=F,row.names=F,col.names=T) 101 write.table(PRODF, "allUnique.csv", sep=",",quote=F,row.names=F,col.names=T)
97 102
98 pV = ggplot(PRODFV) 103 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)) 104 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") 105 pV = pV + xlab("Summary of V gene") + ylab("Frequency") + ggtitle("Relative frequency of V gene usage")
106 write.table(x=PRODFV, file="VFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
101 107
102 png("VPlot.png",width = 1280, height = 720) 108 png("VPlot.png",width = 1280, height = 720)
103 pV 109 pV
104 dev.off(); 110 dev.off();
105 111
106 pD = ggplot(PRODFD) 112 pD = ggplot(PRODFD)
107 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)) 113 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))
108 pD = pD + xlab("Summary of D gene") + ylab("Frequency") + ggtitle("Relative frequency of D gene usage") 114 pD = pD + xlab("Summary of D gene") + ylab("Frequency") + ggtitle("Relative frequency of D gene usage")
115 write.table(x=PRODFD, file="DFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
109 116
110 png("DPlot.png",width = 800, height = 600) 117 png("DPlot.png",width = 800, height = 600)
111 pD 118 pD
112 dev.off(); 119 dev.off();
113 120
114 pJ = ggplot(PRODFJ) 121 pJ = ggplot(PRODFJ)
115 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)) 122 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))
116 pJ = pJ + xlab("Summary of J gene") + ylab("Frequency") + ggtitle("Relative frequency of J gene usage") 123 pJ = pJ + xlab("Summary of J gene") + ylab("Frequency") + ggtitle("Relative frequency of J gene usage")
124 write.table(x=PRODFJ, file="JFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
117 125
118 png("JPlot.png",width = 800, height = 600) 126 png("JPlot.png",width = 800, height = 600)
119 pJ 127 pJ
120 dev.off(); 128 dev.off();
129
130 VGenes = PRODF[,c("Sample", "Top.V.Gene")]
131 VGenes$Top.V.Gene = gsub("-.*", "", VGenes$Top.V.Gene)
132 VGenes = data.frame(data.table(VGenes)[, list(Count=.N), by=c("Sample", "Top.V.Gene")])
133 TotalPerSample = data.frame(data.table(VGenes)[, list(total=sum(.SD$Count)), by=Sample])
134 VGenes = merge(VGenes, TotalPerSample, by="Sample")
135 VGenes$Frequency = VGenes$Count * 100 / VGenes$total
136 VPlot = ggplot(VGenes)
137 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))
138 png("VFPlot.png")
139 VPlot
140 dev.off();
141 write.table(x=VGenes, file="VFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
142
143 DGenes = PRODF[,c("Sample", "Top.D.Gene")]
144 DGenes$Top.D.Gene = gsub("-.*", "", DGenes$Top.D.Gene)
145 DGenes = data.frame(data.table(DGenes)[, list(Count=.N), by=c("Sample", "Top.D.Gene")])
146 TotalPerSample = data.frame(data.table(DGenes)[, list(total=sum(.SD$Count)), by=Sample])
147 DGenes = merge(DGenes, TotalPerSample, by="Sample")
148 DGenes$Frequency = DGenes$Count * 100 / DGenes$total
149 DPlot = ggplot(DGenes)
150 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))
151 png("DFPlot.png")
152 DPlot
153 dev.off();
154 write.table(x=DGenes, file="DFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
155
156 JGenes = PRODF[,c("Sample", "Top.J.Gene")]
157 JGenes$Top.J.Gene = gsub("-.*", "", JGenes$Top.J.Gene)
158 JGenes = data.frame(data.table(JGenes)[, list(Count=.N), by=c("Sample", "Top.J.Gene")])
159 TotalPerSample = data.frame(data.table(JGenes)[, list(total=sum(.SD$Count)), by=Sample])
160 JGenes = merge(JGenes, TotalPerSample, by="Sample")
161 JGenes$Frequency = JGenes$Count * 100 / JGenes$total
162 JPlot = ggplot(JGenes)
163 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))
164 png("JFPlot.png")
165 JPlot
166 dev.off();
167 write.table(x=JGenes, file="JFFrequency.csv", sep=",",quote=F,row.names=F,col.names=T)
121 168
122 revVchain = Vchain 169 revVchain = Vchain
123 revDchain = Dchain 170 revDchain = Dchain
124 revVchain$chr.orderV = rev(revVchain$chr.orderV) 171 revVchain$chr.orderV = rev(revVchain$chr.orderV)
125 revDchain$chr.orderD = rev(revDchain$chr.orderD) 172 revDchain$chr.orderD = rev(revDchain$chr.orderD)
136 xlab("D genes") + 183 xlab("D genes") +
137 ylab("V Genes") 184 ylab("V Genes")
138 185
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))) 186 png(paste("HeatmapVD_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Dchain$v.name)), height=100+(15*length(Vchain$v.name)))
140 print(img) 187 print(img)
188
141 dev.off() 189 dev.off()
190 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)
142 } 191 }
143 192
144 VandDCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.V.Gene", "Top.D.Gene", "Sample")]) 193 VandDCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.V.Gene", "Top.D.Gene", "Sample")])
145 194
146 VandDCount$l = log(VandDCount$Length) 195 VandDCount$l = log(VandDCount$Length)
172 ylab("V Genes") 221 ylab("V Genes")
173 222
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))) 223 png(paste("HeatmapVJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Vchain$v.name)))
175 print(img) 224 print(img)
176 dev.off() 225 dev.off()
226 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)
177 } 227 }
178 228
179 VandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.V.Gene", "Top.J.Gene", "Sample")]) 229 VandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.V.Gene", "Top.J.Gene", "Sample")])
180 230
181 VandJCount$l = log(VandJCount$Length) 231 VandJCount$l = log(VandJCount$Length)
204 ylab("D Genes") 254 ylab("D Genes")
205 255
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))) 256 png(paste("HeatmapDJ_", unique(dat[3])[1,1] , ".png", sep=""), width=150+(15*length(Jchain$v.name)), height=100+(15*length(Dchain$v.name)))
207 print(img) 257 print(img)
208 dev.off() 258 dev.off()
259 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)
209 } 260 }
210 261
211 DandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.D.Gene", "Top.J.Gene", "Sample")]) 262 DandJCount = data.frame(data.table(PRODF)[, list(Length=.N), by=c("Top.D.Gene", "Top.J.Gene", "Sample")])
212 263
213 DandJCount$l = log(DandJCount$Length) 264 DandJCount$l = log(DandJCount$Length)
234 if("Replicate" %in% colnames(test)) 285 if("Replicate" %in% colnames(test))
235 { 286 {
236 clonalityFrame = PROD 287 clonalityFrame = PROD
237 clonalityFrame$ReplicateConcat = do.call(paste, c(clonalityFrame[c("VDJCDR3", "Sample", "Replicate")], sep = ":")) 288 clonalityFrame$ReplicateConcat = do.call(paste, c(clonalityFrame[c("VDJCDR3", "Sample", "Replicate")], sep = ":"))
238 clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ] 289 clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ]
239 write.table(clonalityFrame, "clonalityComplete.tsv", sep="\t",quote=F,row.names=F,col.names=T) 290 write.table(clonalityFrame, "clonalityComplete.csv", sep=",",quote=F,row.names=F,col.names=T)
240 291
241 ClonalitySampleReplicatePrint <- function(dat){ 292 ClonalitySampleReplicatePrint <- function(dat){
242 write.table(dat, paste("clonality_", unique(dat$Sample) , "_", unique(dat$Replicate), ".tsv", sep=""), sep="\t",quote=F,row.names=F,col.names=T) 293 write.table(dat, paste("clonality_", unique(dat$Sample) , "_", unique(dat$Replicate), ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=T)
243 } 294 }
244 295
245 clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,c("Sample", "Replicate")]) 296 clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,c("Sample", "Replicate")])
246 #lapply(clonalityFrameSplit, FUN=ClonalitySampleReplicatePrint) 297 #lapply(clonalityFrameSplit, FUN=ClonalitySampleReplicatePrint)
247 298
248 ClonalitySamplePrint <- function(dat){ 299 ClonalitySamplePrint <- function(dat){
249 write.table(dat, paste("clonality_", unique(dat$Sample) , ".tsv", sep=""), sep="\t",quote=F,row.names=F,col.names=T) 300 write.table(dat, paste("clonality_", unique(dat$Sample) , ".csv", sep=""), sep=",",quote=F,row.names=F,col.names=T)
250 } 301 }
251 302
252 clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,"Sample"]) 303 clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,"Sample"])
253 #lapply(clonalityFrameSplit, FUN=ClonalitySamplePrint) 304 #lapply(clonalityFrameSplit, FUN=ClonalitySamplePrint)
254 305