Mercurial > repos > davidvanzessen > clonal_sequences_in_paired_samples

comparison RScript.r @ 2:8d562506f4f9 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Mon, 25 Aug 2014 03:38:07 -0400
parents c5ac9a871b26
children f9316f7676cc
comparison
equal deleted inserted replaced
1:c0110597898e 2:8d562506f4f9
6 require(ggplot2) 6 require(ggplot2)
7 require(reshape2) 7 require(reshape2)
8 require(data.table) 8 require(data.table)
9 require(grid) 9 require(grid)
10 #require(xtable) 10 #require(xtable)
11 #dat = read.csv(inFile) 11 dat = read.csv(inFile, sep="\t")
12 dat = data.frame(fread(inFile)) #faster but with a dep 12 #dat = data.frame(fread(inFile)) #faster but with a dep
13 setwd(outDir) 13 setwd(outDir)
14 dat$V_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$V_Segment_Major_Gene), "; "), "[[", 1)))
15 dat$J_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$J_Segment_Major_Gene), "; "), "[[", 1)))
16
14 dat$Frequency = ((10^dat$Log10_Frequency)*100) 17 dat$Frequency = ((10^dat$Log10_Frequency)*100)
18
19 dat$normalized_read_count = round(dat$Clone_Molecule_Count_From_Spikes / dat$Cell_Count * 1000000 / 2)
20 dat$paste = paste(dat$Sample, dat$V_Segment_Major_Gene, dat$J_Segment_Major_Gene, dat$CDR3_Sense_Sequence)
21 dat = dat[!duplicated(dat$paste),]
15 patients = split(dat, dat$Patient, drop=T) 22 patients = split(dat, dat$Patient, drop=T)
16 intervalReads = rev(c(0,2,10,100,1000,10000)) 23 intervalReads = rev(c(0,2,10,100,1000,10000))
17 intervalFreq = rev(c(0,0.01,0.1,0.5,1,5)) 24 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") 25 V_Segments = c(".*", "IGHV", "IGHD", "IGKV", "IGKV", "IgKINTR", "TRGV", "TRDV", "TRDD" , "TRBV")
19 J_Segments = c(".*", ".*", ".*", "IGKJ", "KDE", ".*", ".*", ".*", ".*", ".*") 26 J_Segments = c(".*", ".*", ".*", "IGKJ", "KDE", ".*", ".*", ".*", ".*", ".*")
26 onShort = "reads" 33 onShort = "reads"
27 if(on == "Frequency"){ 34 if(on == "Frequency"){
28 onShort = "freq" 35 onShort = "freq"
29 } 36 }
30 splt = split(x, x$Sample, drop=T) 37 splt = split(x, x$Sample, drop=T)
38 if(length(splt) == 1){
39 print(paste(paste(x[1,which(colnames(x) == "Patient")]), "has one sample, skipping"))
40 return()
41 }
31 patient1 = splt[[1]] 42 patient1 = splt[[1]]
32 patient2 = splt[[2]] 43 patient2 = splt[[2]]
33 44
34 threshholdIndex = which(colnames(product) == "interval") 45 threshholdIndex = which(colnames(product) == "interval")
35 V_SegmentIndex = which(colnames(product) == "V_Segments") 46 V_SegmentIndex = which(colnames(product) == "V_Segments")
61 res1 = vector() 72 res1 = vector()
62 res2 = vector() 73 res2 = vector()
63 resBoth = vector() 74 resBoth = vector()
64 read1Count = vector() 75 read1Count = vector()
65 read2Count = vector() 76 read2Count = vector()
77 locussum1 = vector()
78 locussum2 = vector()
66 79
67 print(c(patient, " ", on)) 80 print(patient)
68 #for(iter in 1){ 81 #for(iter in 1){
69 for(iter in 1:length(product[,1])){ 82 for(iter in 1:length(product[,1])){
70 threshhold = product[iter,threshholdIndex] 83 threshhold = product[iter,threshholdIndex]
71 V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="") 84 V_Segment = paste(".*", as.character(product[iter,V_SegmentIndex]), ".*", sep="")
72 J_Segment = paste(".*", as.character(product[iter,J_SegmentIndex]), ".*", sep="") 85 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) 86 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)) 87 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)) 88 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)) 89 read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count) + sum(patientMerge[both,]$normalized_read_count.x))
77 read2Count = append(read2Count, sum(patient2[two,]$Clone_Molecule_Count_From_Spikes) + sum(patientMerge[both,]$Clone_Molecule_Count_From_Spikes.y)) 90 read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count) + sum(patientMerge[both,]$normalized_read_count.y))
78 res1 = append(res1, sum(one)) 91 res1 = append(res1, sum(one))
79 res2= append(res2, sum(two)) 92 res2 = append(res2, sum(two))
80 resBoth = append(resBoth, sum(both)) 93 resBoth = append(resBoth, sum(both))
94 locussum1 = append(locussum1, sum(patient1[(grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene)),]$normalized_read_count))
95 locussum2 = append(locussum2, sum(patient2[(grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene)),]$normalized_read_count))
81 #threshhold = 0 96 #threshhold = 0
82 if(threshhold != 0){ 97 if(threshhold != 0){
83 if(sum(one) > 0){ 98 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")] 99 dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence")]
85 colnames(dfOne) = c("Proximal segment", "Distal segment", "Clone_Molecule_Count_From_Spikes", "Frequency", "Sequence") 100 colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Sequence")
86 filenameOne = paste(oneSample, "_", product[iter, titleIndex], "_", threshhold, sep="") 101 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) 102 write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
89 } 103 }
90 if(sum(two) > 0){ 104 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")] 105 dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence")]
92 colnames(dfTwo) = c("Proximal segment", "Distal segment", "Clone_Molecule_Count_From_Spikes", "Frequency", "Sequence") 106 colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Sequence")
93 filenameTwo = paste(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="") 107 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) 108 write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
96 } 109 }
97 } 110 }
98 if(sum(both) > 0){ 111 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")] 112 dfBoth = patientMerge[both,c("V_Segment_Major_Gene.x", "J_Segment_Major_Gene.x", "normalized_read_count.x", "Frequency.x", "Clone_Sequence", "V_Segment_Major_Gene.y", "J_Segment_Major_Gene.y", "normalized_read_count.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)) 113 colnames(dfBoth) = c(paste("Proximal segment", oneSample), paste("Distal segment", oneSample), paste("Normalized_Read_Count", oneSample), paste("Frequency", oneSample), "Sequence", paste("Proximal segment", twoSample), paste("Distal segment", twoSample), paste("Normalized_Read_Count", twoSample), paste("Frequency", twoSample))
101 filenameBoth = paste(oneSample, "_", twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="") 114 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) 115 write.table(dfBoth, file=paste(filenameBoth, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
104 } 116 }
105 } 117 }
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)) 118 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), "Locus_sum1"=locussum1, "Locus_sum2"=locussum2)
107 if(sum(is.na(patientResult$percentage)) > 0){ 119 if(sum(is.na(patientResult$percentage)) > 0){
108 patientResult[is.na(patientResult$percentage),]$percentage = 0 120 patientResult[is.na(patientResult$percentage),]$percentage = 0
109 } 121 }
110 colnames(patientResult)[6] = oneSample 122 colnames(patientResult)[6] = oneSample
111 colnames(patientResult)[8] = twoSample 123 colnames(patientResult)[8] = twoSample
112 colnamesBak = colnames(patientResult) 124 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="")) 125 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("Normalized Read Count", oneSample), paste("Number of sequences", twoSample, sep=""), paste("Normalized Read Count", twoSample), paste("Sum number of sequences", patient), paste("Percentage of sequences ", patient, "_Both", sep=""), paste("Locus Sum", oneSample), paste("Locus Sum", twoSample))
114 write.table(patientResult, file=paste(patient, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T) 126 write.table(patientResult, file=paste(patient, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
115 colnames(patientResult) = colnamesBak 127 colnames(patientResult) = colnamesBak
116 128
117 patientResult$Locus = factor(patientResult$Locus, Titles) 129 patientResult$Locus = factor(patientResult$Locus, Titles)
118 patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep="")) 130 patientResult$cut_off_value = factor(patientResult$cut_off_value, paste(">", interval, sep=""))