--- 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