Mercurial > repos > davidvanzessen > clonal_sequences_in_paired_samples

changeset 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 c0110597898e
children f9316f7676cc
files RScript.r wrapper.sh
diffstat 2 files changed, 42 insertions(+), 24 deletions(-) [+]
line wrap: on
line diff
--- a/RScript.r	Thu Jun 19 11:12:44 2014 -0400
+++ b/RScript.r	Mon Aug 25 03:38:07 2014 -0400
@@ -8,10 +8,17 @@
 require(data.table)
 require(grid)
 #require(xtable)
-#dat = read.csv(inFile)
-dat = data.frame(fread(inFile)) #faster but with a dep
+dat = read.csv(inFile, sep="\t")
+#dat = data.frame(fread(inFile)) #faster but with a dep
 setwd(outDir)
+dat$V_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$V_Segment_Major_Gene), "; "), "[[", 1)))
+dat$J_Segment_Major_Gene = as.factor(as.character(lapply(strsplit(as.character(dat$J_Segment_Major_Gene), "; "), "[[", 1)))
+
 dat$Frequency = ((10^dat$Log10_Frequency)*100)
+
+dat$normalized_read_count = round(dat$Clone_Molecule_Count_From_Spikes / dat$Cell_Count * 1000000 / 2)
+dat$paste = paste(dat$Sample, dat$V_Segment_Major_Gene, dat$J_Segment_Major_Gene, dat$CDR3_Sense_Sequence)
+dat = dat[!duplicated(dat$paste),]
 patients = split(dat, dat$Patient, drop=T)
 intervalReads = rev(c(0,2,10,100,1000,10000))
 intervalFreq = rev(c(0,0.01,0.1,0.5,1,5))
@@ -28,6 +35,10 @@
     onShort = "freq"
   }
   splt = split(x, x$Sample, drop=T)
+  if(length(splt) == 1){
+    print(paste(paste(x[1,which(colnames(x) == "Patient")]), "has one sample, skipping"))
+    return()
+  }
   patient1 = splt[[1]]
   patient2 = splt[[2]]
   
@@ -63,8 +74,10 @@
   resBoth = vector()
   read1Count = vector()
   read2Count = vector()
+  locussum1 = vector()
+  locussum2 = vector()
   
-  print(c(patient, " ", on))
+  print(patient)
   #for(iter in 1){
   for(iter in 1:length(product[,1])){
     threshhold = product[iter,threshholdIndex]
@@ -73,44 +86,43 @@
     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)
     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))
     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))
-    read1Count = append(read1Count, sum(patient1[one,]$Clone_Molecule_Count_From_Spikes) + sum(patientMerge[both,]$Clone_Molecule_Count_From_Spikes.x))
-    read2Count = append(read2Count, sum(patient2[two,]$Clone_Molecule_Count_From_Spikes) + sum(patientMerge[both,]$Clone_Molecule_Count_From_Spikes.y))
+    read1Count = append(read1Count, sum(patient1[one,]$normalized_read_count) + sum(patientMerge[both,]$normalized_read_count.x))
+    read2Count = append(read2Count, sum(patient2[two,]$normalized_read_count) + sum(patientMerge[both,]$normalized_read_count.y))
     res1 = append(res1, sum(one))
-    res2= append(res2, sum(two))
+    res2 = append(res2, sum(two))
     resBoth = append(resBoth, sum(both))
+    locussum1 = append(locussum1, sum(patient1[(grepl(V_Segment, patient1$V_Segment_Major_Gene) & grepl(J_Segment, patient1$J_Segment_Major_Gene)),]$normalized_read_count))
+    locussum2 = append(locussum2, sum(patient2[(grepl(V_Segment, patient2$V_Segment_Major_Gene) & grepl(J_Segment, patient2$J_Segment_Major_Gene)),]$normalized_read_count))
     #threshhold = 0
     if(threshhold != 0){
       if(sum(one) > 0){
-        dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "Clone_Molecule_Count_From_Spikes", "Frequency", "Clone_Sequence")]
-        colnames(dfOne) = c("Proximal segment", "Distal segment", "Clone_Molecule_Count_From_Spikes", "Frequency", "Sequence")
+        dfOne = patient1[one,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence")]
+        colnames(dfOne) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Sequence")
         filenameOne = paste(oneSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
-        #print(xtable(dfOne),type="html", file=paste(filenameOne, ".html", sep=""), include.rownames=FALSE)
         write.table(dfOne, file=paste(filenameOne, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
       }
       if(sum(two) > 0){
-        dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "Clone_Molecule_Count_From_Spikes", "Frequency", "Clone_Sequence")]
-        colnames(dfTwo) = c("Proximal segment", "Distal segment", "Clone_Molecule_Count_From_Spikes", "Frequency", "Sequence")
+        dfTwo = patient2[two,c("V_Segment_Major_Gene", "J_Segment_Major_Gene", "normalized_read_count", "Frequency", "Clone_Sequence")]
+        colnames(dfTwo) = c("Proximal segment", "Distal segment", "normalized_read_count", "Frequency", "Sequence")
         filenameTwo = paste(twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
-        #print(xtable(dfTwo),type="html", file=paste(filenameTwo, ".html", sep=""), include.rownames=FALSE)
         write.table(dfTwo, file=paste(filenameTwo, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
       }
     }
     if(sum(both) > 0){
-      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")]
-      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))
+      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")]
+      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))
       filenameBoth = paste(oneSample, "_", twoSample, "_", product[iter, titleIndex], "_", threshhold, sep="")
-      #print(xtable(dfBoth),type="html", file=paste(filenameBoth, ".html", sep=""), include.rownames=FALSE)
       write.table(dfBoth, file=paste(filenameBoth, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
     }
   }
-  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))
+  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)
   if(sum(is.na(patientResult$percentage)) > 0){
     patientResult[is.na(patientResult$percentage),]$percentage = 0
   }
   colnames(patientResult)[6] = oneSample
   colnames(patientResult)[8] = twoSample
   colnamesBak = colnames(patientResult)
-  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=""))
+  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))
   write.table(patientResult, file=paste(patient, "_", onShort, ".txt", sep=""), quote=F, sep="\t", dec=",", row.names=F, col.names=T)
   colnames(patientResult) = colnamesBak
   
--- a/wrapper.sh	Thu Jun 19 11:12:44 2014 -0400
+++ b/wrapper.sh	Mon Aug 25 03:38:07 2014 -0400
@@ -38,13 +38,15 @@
 	echo "<div class='tabbertab' title='Data frequency'>" >> $html
 	echo "<table><tr><td style='vertical-align:top;'>" >> $html
 	echo "<table border = 1 class='result_table summary_table' id='summary_table_${patient}_freq'>" >> $html
-	echo "<thead><th>Ig/TCR gene rearrangement type</th><th>Proximal gene segment</th><th>Distal gene segment</th><th>Cut off value</th><th>Number of sequences ${patient}_Both</th><th>Number of sequences_$sample1</th><th>Read Count $sample1</th><th>Number of sequences_$sample2</th><th>Read Count $sample2</th><th>Sum number of sequences $patient</th><th>Percentage of sequences ${patient}_both</th></thead>" >> $html
+	echo "<thead><th>Ig/TCR gene rearrangement type</th><th>Proximal gene segment</th><th>Distal gene segment</th><th>Cut off value</th><th>Number of sequences ${patient}_Both</th><th>Number of sequences_$sample1</th><th>Normalized Read Count $sample1</th><th>Number of sequences_$sample2</th><th>Normalized Read Count $sample2</th><th>Sum number of sequences $patient</th><th>Percentage of sequences ${patient}_both</th></thead>" >> $html
 	echo "<tbody>" >> $html
-	while read locus j_segment v_segment cut_off_value both one read_count1 two read_count2 sum percent
+	readsumtable="<table class='result_table summary_table'><thead><tr><th>Ig/TCR gene rearrangement type</th><th>Proximal gene segment</th><th>Distal gene segment</th><th>Total normalized read count for $sample1</th><th>Total normalized read count for $sample2</th></tr></thead>"
+	while read locus j_segment v_segment cut_off_value both one read_count1 two read_count2 sum percent locusreadsum1 locusreadsum2
 	do
 		if [ "$locus" != "$oldLocus" ] ; then
 			echo "<tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr>" >> $html
 			echo "<tr><td><b>$locus</b></td>" >> $html
+			readsumtable="${readsumtable}<tr><td>$locus</td><td>$v_segment</td><td>$j_segment</td><td>$locusreadsum1</td><td>$locusreadsum2</td></tr>"
 		else
 			echo "<td></td>" >> $html
 		fi
@@ -71,10 +73,11 @@
 		echo "<td>$sum</td>" >> $html
 		echo "<td>${percent}&#37;</td>" >> $html
 		echo "</tr>" >> $html
-		oldLocus="$locus"
+		oldLocus="$locus"	
 	done < tmp.txt
 	echo "</tbody></table>" >> $html
-	echo "</td><td style='vertical-align:top;'><div id='result_div_${patient}_freq'></div></td></tr></table></div>" >> $html
+	echo "</td><td style='vertical-align:top;'><div id='result_div_${patient}_freq'></div></td></tr></table>" >> $html
+	echo "${readsumtable}</table></div>" >> $html
 	echo "<div class='tabbertab' title='Graphs frequency'>" >> $html
 	echo "<a href='${patient}_freq.png'><img src='${patient}_freq.png' width='1280' height='720' /></a><br />" >> $html
 	echo "<a href='${patient}_freq_both.png'><img src='${patient}_freq_both.png' width='1280' height='720' /></a><br />" >> $html
@@ -86,11 +89,13 @@
 	echo "<table border = 1 class='result_table summary_table' id='summary_table_${patient}_reads'>" >> $html
 	echo "<thead><th>Ig/TCR gene rearrangement type</th><th>Proximal gene segment</th><th>Distal gene segment</th><th>Cut off value</th><th>Number of sequences ${patient}_Both</th><th>Number of sequences_$sample1</th><th>Read Count $sample1</th><th>Number of sequences_$sample2</th><th>Read Count $sample2</th><th>Sum number of sequences $patient</th><th>Percentage of sequences ${patient}_both</th></thead>" >> $html
 	echo "<tbody>" >> $html
-	while read locus j_segment v_segment cut_off_value both one read_count1 two read_count2 sum percent
+	readsumtable="<table class='result_table summary_table'><thead><tr><th>Ig/TCR gene rearrangement type</th><th>Proximal gene segment</th><th>Distal gene segment</th><th>Total normalized read count for $sample1</th><th>Total normalized read count for $sample2</th></tr></thead>"
+	while read locus j_segment v_segment cut_off_value both one read_count1 two read_count2 sum percent locusreadsum1 locusreadsum2
 	do
 		if [ "$locus" != "$oldLocus" ] ; then
 			echo "<tr><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td><td></td></tr><tr>" >> $html
 			echo "<tr><td><b>$locus</b></td>" >> $html
+			readsumtable="${readsumtable}<tr><td>$locus</td><td>$v_segment</td><td>$j_segment</td><td>$locusreadsum1</td><td>$locusreadsum2</td></tr>"
 		else
 			echo "<td></td>" >> $html
 		fi
@@ -120,12 +125,13 @@
 		oldLocus="$locus"
 	done < tmp.txt
 	echo "</tbody></table>" >> $html
-	echo "</td><td style='vertical-align:top;'><div id='result_div_${patient}_reads'></div></td></tr></table></div>" >> $html
+	echo "</td><td style='vertical-align:top;'><div id='result_div_${patient}_reads'></div></td></tr></table>" >> $html
+	echo "${readsumtable}</table></div>" >> $html
 	echo "<div class='tabbertab' title='Graphs reads'>" >> $html
 	echo "<a href='${patient}_reads.png'><img src='${patient}_reads.png' width='1280' height='720' /></a><br />" >> $html
 	echo "<a href='${patient}_reads_both.png'><img src='${patient}_reads_both.png' width='1280' height='720' /></a><br />" >> $html
 	echo "<a href='${patient}_percent_reads.png'><img src='${patient}_percent_reads.png' width='1280' height='720' /></a></div>" >> $html
-	echo "</div>" >> $html	
+	echo "</div>" >> $html
 	echo "</div>" >> $html
 done < patients.txt
 echo "</div></html>" >> $html