Mercurial > repos > davidvanzessen > clonal_sequences_in_paired_samples

comparison RScript.r @ 0:c5ac9a871b26 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Thu, 19 Jun 2014 10:39:41 -0400
parents
children 8d562506f4f9
comparison
equal deleted inserted replaced
-1:000000000000 0:c5ac9a871b26
1 args <- commandArgs(trailingOnly = TRUE)
2
3 inFile = args[1]
4 outDir = args[2]
5
6 require(ggplot2)
7 require(reshape2)
8 require(data.table)
9 require(grid)
10 #require(xtable)
11 #dat = read.csv(inFile)
12 dat = data.frame(fread(inFile)) #faster but with a dep
13 setwd(outDir)
14 dat$Frequency = ((10^dat$Log10_Frequency)*100)
15 patients = split(dat, dat$Patient, drop=T)
16 intervalReads = rev(c(0,2,10,100,1000,10000))
17 intervalFreq = rev(c(0,0.01,0.1,0.5,1,5))
18 V_Segments = c(".*", "IGHV", "IGHD", "IGKV", "IGKV", "IgKINTR", "TRGV", "TRDV", "TRDD" , "TRBV")
19 J_Segments = c(".*", ".*", ".*", "IGKJ", "KDE", ".*", ".*", ".*", ".*", ".*")
20 Titles = c("Total", "IGH-Vh-Jh", "IGH-Dh-Jh", "Vk-Jk", "Vk-Kde" , "Intron-Kde", "TCRG", "TCRD-Vd-Dd", "TCRD-Dd-Dd", "TCRB-Vb-Jb")
21 Titles = factor(Titles, levels=Titles)
22 TitlesOrder = data.frame("Title"=Titles, "TitlesOrder"=1:length(Titles))
23
24 patientCountOnColumn <- function(x, product, interval, on, appendtxt=F){
25 x$Sample = factor(x$Sample, levels=unique(x$Sample))
26 onShort = "reads"
27 if(on == "Frequency"){
28 onShort = "freq"
29 }
30 splt = split(x, x$Sample, drop=T)
31 patient1 = splt[[1]]
32 patient2 = splt[[2]]
33
34 threshholdIndex = which(colnames(product) == "interval")
35 V_SegmentIndex = which(colnames(product) == "V_Segments")
36 J_SegmentIndex = which(colnames(product) == "J_Segments")
37 titleIndex = which(colnames(product) == "Titles")
38 sampleIndex = which(colnames(x) == "Sample")
39 patientIndex = which(colnames(x) == "Patient")
40 oneSample = paste(patient1[1,sampleIndex], sep="")
41 twoSample = paste(patient2[1,sampleIndex], sep="")
42 patient = paste(x[1,patientIndex])
43
44 #print(c(length(grep(".*_Right$", twoSample)) == 1, "Right"))
45 #print(c(length(grep(".*_Dx_BM$", twoSample)) == 1, "Dx_BM"))
46 #print(c(length(grep(".*_Dx$", twoSample)) == 1, "Dx"))
47 switched = F
48 if(length(grep(".*_Right$", twoSample)) == 1 || length(grep(".*_Dx_BM$", twoSample)) == 1 || length(grep(".*_Dx$", twoSample)) == 1 ){
49 tmp = twoSample
50 twoSample = oneSample
51 oneSample = tmp
52 tmp = patient1
53 patient1 = patient2
54 patient2 = tmp
55 switched = T
56 }
57 if(appendtxt){
58 cat(paste(patient, oneSample, twoSample, sep="\t"), file="patients.txt", append=T, sep="", fill=3)
59 }
60 patientMerge = merge(patient1, patient2, by="Clone_Sequence")
61 res1 = vector()
62 res2 = vector()
63 resBoth = vector()
64 read1Count = vector()
65 read2Count = vector()
66
67 print(c(patient, " ", on))
68 #for(iter in 1){
69 for(iter in 1:length(product[,1])){
70 threshhold = product[iter,threshholdIndex]
71 V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="")
72 J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="")
73 both = (grepl(V_Segment, patientMerge$V_Segment_Major_Gene.x) & grepl(J_Segment, patientMerge$J_Segment_Major_Gene.x) & patientMerge[,paste(on, ".x", sep="")] > threshhold & patientMerge[,paste(on, ".y", sep="")] > threshhold)
74 one = (grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene) & patient1[,on] > threshhold & !(patient1$Clone_Sequence %in% patientMerge[both,]$Clone_Sequence))
75 two = (grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene) & patient2[,on] > threshhold & !(patient2$Clone_Sequence %in% patientMerge[both,]$Clone_Sequence))
76 read1Count = append(read1Count, sum(patient1[one,]$Clone_Molecule_Count_From_Spikes) + sum(patientMerge[both,]$Clone_Molecule_Count_From_Spikes.x))
77 read2Count = append(read2Count, sum(patient2[two,]$Clone_Molecule_Count_From_Spikes) + sum(patientMerge[both,]$Clone_Molecule_Count_From_Spikes.y))
78 res1 = append(res1, sum(one))
79 res2= append(res2, sum(two))
80 resBoth = append(resBoth, sum(both))
81 #threshhold = 0
82 if(threshhold != 0){
83 if(sum(one) > 0){
84 dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "Clone_Molecule_Count_From_Spikes", "Frequency", "Clone_Sequence")]
85 colnames(dfOne) = c("Proximal segment", "Distal segment", "Clone_Molecule_Count_From_Spikes", "Frequency", "Sequence")
86 filenameOne = paste(oneSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
87 #print(xtable(dfOne),type="html", file=paste(filenameOne, ".html", sep=""), include.rownames=FALSE)
88 write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
89 }
90 if(sum(two) > 0){
91 dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "Clone_Molecule_Count_From_Spikes", "Frequency", "Clone_Sequence")]
92 colnames(dfTwo) = c("Proximal segment", "Distal segment", "Clone_Molecule_Count_From_Spikes", "Frequency", "Sequence")
93 filenameTwo = paste(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
94 #print(xtable(dfTwo),type="html", file=paste(filenameTwo, ".html", sep=""), include.rownames=FALSE)
95 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
96 }
97 }
98 if(sum(both) > 0){
99 dfBoth = patientMerge[both,c("V_Segment_Major_Gene.x", "J_Segment_Major_Gene.x", "Clone_Molecule_Count_From_Spikes.x", "Frequency.x", "Clone_Sequence", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "Clone_Molecule_Count_From_Spikes.y", "Frequency.y")]
100 colnames(dfBoth) = c(paste("Proximal segment", oneSample), paste("Distal segment", oneSample), paste("Read Count", oneSample), paste("Frequency", oneSample), "Sequence", paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Read_Count", twoSample), paste("Frequency", twoSample))
101 filenameBoth = paste(oneSample, "_", twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
102 #print(xtable(dfBoth),type="html", file=paste(filenameBoth, ".html", sep=""), include.rownames=FALSE)
103 write.table(dfBoth, file=paste(filenameBoth, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
104 }
105 }
106 patientResult = data.frame("Locus"=product$Titles, "J_Segment"=product$J_Segments, "V_Segment"=product$V_Segments, "cut_off_value"=paste(">", product$interval, sep=""), "Both"=resBoth, "tmp1"=res1, "read_count1" = round(read1Count), "tmp2"=res2, "read_count2"= round(read2Count), "Sum"=res1 + res2 + resBoth, "percentage" = round((resBoth/(res1 + res2 + resBoth)) * 100, digits=2))
107 if(sum(is.na(patientResult$percentage)) > 0){
108 patientResult[is.na(patientResult$percentage),]$percentage = 0
109 }
110 colnames(patientResult)[6] = oneSample
111 colnames(patientResult)[8] = twoSample
112 colnamesBak = colnames(patientResult)
113 colnames(patientResult) = c("Ig/TCR gene rearrangement type", "Distal Gene segment", "Proximal gene segment", "cut_off_value", paste("Number of sequences ", patient, "_Both", sep=""), paste("Number of sequences", oneSample, sep=""), paste("Read Count", oneSample), paste("Number of sequences", twoSample, sep=""), paste("Read Count", twoSample), paste("Sum number of sequences", patient), paste("Percentage of sequences ", patient, "_Both", sep=""))
114 write.table(patientResult, file=paste(patient, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
115 colnames(patientResult) = colnamesBak
116
117 patientResult$Locus = factor(patientResult$Locus, Titles)
118 patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))
119
120 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "Both")])
121 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=Both), stat='identity', position="dodge", fill="#79c36a")
122 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
123 plt = plt + geom_text(aes(ymax=max(Both), x=cut_off_value,y=Both,label=Both), angle=90, hjust=0)
124 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("Number of clones in both")
125 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
126 png(paste(patient, "_", onShort, ".png", sep=""), width=1920, height=1080)
127 print(plt)
128 dev.off()
129 #(t,r,b,l)
130 plt = ggplot(patientResult[,c("Locus", "cut_off_value", "percentage")])
131 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=percentage), stat='identity', position="dodge", fill="#79c36a")
132 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
133 plt = plt + geom_text(aes(ymax=max(percentage), x=cut_off_value,y=percentage,label=percentage), angle=90, hjust=0)
134 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle("% clones in both left and right")
135 plt = plt + theme(plot.margin = unit(c(1,8.8,0.5,1.5), "lines"))
136 png(paste(patient, "_percent_", onShort, ".png", sep=""), width=1920, height=1080)
137 print(plt)
138 dev.off()
139
140 patientResult = melt(patientResult[,c('Locus','cut_off_value', oneSample, twoSample)] ,id.vars=1:2)
141 patientResult$relativeValue = patientResult$value * 10
142 patientResult[patientResult$relativeValue == 0,]$relativeValue = 1
143 plt = ggplot(patientResult)
144 plt = plt + geom_bar( aes( x=factor(cut_off_value), y=relativeValue, fill=variable), stat='identity', position="dodge")
145 plt = plt + facet_grid(.~Locus) + theme(axis.text.x = element_text(angle = 45, hjust = 1))
146 plt = plt + scale_y_continuous(trans="log", breaks=10^c(0:10), labels=c(0, 10^c(0:9)))
147 plt = plt + geom_text(data=patientResult[patientResult$variable == oneSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=-0.2)
148 plt = plt + geom_text(data=patientResult[patientResult$variable == twoSample,], aes(ymax=max(value), x=cut_off_value,y=relativeValue,label=value), angle=90, position=position_dodge(width=0.9), hjust=0, vjust=0.8)
149 plt = plt + xlab("Reads per locus") + ylab("Count") + ggtitle(paste("Number of clones in only ", oneSample, " and only ", twoSample, sep=""))
150 png(paste(patient, "_", onShort, "_both.png", sep=""), width=1920, height=1080)
151 print(plt)
152 dev.off()
153 }
154
155 interval = intervalFreq
156 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
157 product = data.frame("Titles"=rep(Titles, each=6), "interval"=rep(interval, times=10), "V_Segments"=rep(V_Segments, each=6), "J_Segments"=rep(J_Segments, each=6))
158 #patientFrequencyCount(patient1)
159 #lapply(patients[c(5,6,10)], FUN=patientFrequencyCount)
160 #lapply(patients[c(5,6,7,8,13)], FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
161 #lapply(patients[c(6,7,8)], FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
162 #lapply(patients[c(6)], FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
163 lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Frequency", appendtxt=T)
164
165 interval = intervalReads
166 intervalOrder = data.frame("interval"=paste(">", interval, sep=""), "intervalOrder"=1:length(interval))
167 product = data.frame("Titles"=rep(Titles, each=6), "interval"=rep(interval, times=10), "V_Segments"=rep(V_Segments, each=6), "J_Segments"=rep(J_Segments, each=6))
168 #patientResult = patientReadCount(patient1)
169 #lapply(patients[c(5,6,10)], FUN=patientReadCount)
170 #lapply(patients[c(5,6,7,8,13)], FUN=patientCountOnColumn, product = product, interval=interval, on="Clone_Molecule_Count_From_Spikes")
171 #lapply(patients[c(6)], FUN=patientCountOnColumn, product = product, interval=interval, on="Clone_Molecule_Count_From_Spikes")
172 lapply(patients, FUN=patientCountOnColumn, product = product, interval=interval, on="Clone_Molecule_Count_From_Spikes")
173