--- a/report_clonality/RScript.r	Thu Feb 09 07:20:09 2017 -0500
+++ b/report_clonality/RScript.r	Mon Feb 13 09:08:46 2017 -0500
@@ -614,86 +614,86 @@
       colnames(coincidence.table) = c("Coincidence Type",  "Raw Coincidence Freq")
       write.table(coincidence.table, file=paste("lymphclon_coincidences_", sample_id, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
     }
-  } else if(clonality_method == "old") {
-    clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "clonaltype")])
-    
-    #write files for every coincidence group of >1
-    samples = unique(clonalFreq$Sample)
-    for(sample in samples){
-		clonalFreqSample = clonalFreq[clonalFreq$Sample == sample,]
-		if(max(clonalFreqSample$Type) > 1){
-			for(i in 2:max(clonalFreqSample$Type)){
-				clonalFreqSampleType = clonalFreqSample[clonalFreqSample$Type == i,]
-				clonalityFrame.sub = clonalityFrame[clonalityFrame$clonaltype %in% clonalFreqSampleType$clonaltype,]
-				clonalityFrame.sub = clonalityFrame.sub[order(clonalityFrame.sub$clonaltype),]
-				write.table(clonalityFrame.sub, file=paste("coincidences_", sample, "_", i, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
-			}
-		}
-	}
-	
-    clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")])
-    clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count
-    clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")])
-    clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample")
-    
-    ct = c('Type\tWeight\n2\t1\n3\t3\n4\t6\n5\t10\n6\t15')
-    tcct = textConnection(ct)
-    CT  = read.table(tcct, sep="\t", header=TRUE)
-    close(tcct)
-    clonalFreqCount = merge(clonalFreqCount, CT, by.x="Type", by.y="Type", all.x=T)
-    clonalFreqCount$WeightedCount = clonalFreqCount$Count * clonalFreqCount$Weight
-    
-    ReplicateReads = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "Replicate", "clonaltype")])
-    ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")])
-    clonalFreqCount$Reads = as.numeric(clonalFreqCount$Reads)
-    ReplicateReads$Reads = as.numeric(ReplicateReads$Reads)
-    ReplicateReads$squared = as.numeric(ReplicateReads$Reads * ReplicateReads$Reads)
-    
-    ReplicatePrint <- function(dat){
-      write.table(dat[-1], paste("ReplicateReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
-    }
-    
-    ReplicateSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"])
-    lapply(ReplicateSplit, FUN=ReplicatePrint)
-    
-    ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(as.numeric(Reads)), ReadsSquaredSum=sum(as.numeric(squared))), by=c("Sample")])
-    clonalFreqCount = merge(clonalFreqCount, ReplicateReads, by.x="Sample", by.y="Sample", all.x=T)
-    
-    ReplicateSumPrint <- function(dat){
-      write.table(dat[-1], paste("ReplicateSumReads_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
-    }
-    
-    ReplicateSumSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"])
-    lapply(ReplicateSumSplit, FUN=ReplicateSumPrint)
-    
-    clonalFreqCountSum = data.frame(data.table(clonalFreqCount)[, list(Numerator=sum(WeightedCount, na.rm=T)), by=c("Sample")])
-    clonalFreqCount = merge(clonalFreqCount, clonalFreqCountSum, by.x="Sample", by.y="Sample", all.x=T)
-    clonalFreqCount$ReadsSum = as.numeric(clonalFreqCount$ReadsSum) #prevent integer overflow
-    clonalFreqCount$Denominator = (((clonalFreqCount$ReadsSum * clonalFreqCount$ReadsSum) - clonalFreqCount$ReadsSquaredSum) / 2)
-    clonalFreqCount$Result = (clonalFreqCount$Numerator + 1) / (clonalFreqCount$Denominator + 1)
-    
-    ClonalityScorePrint <- function(dat){
-      write.table(dat$Result, paste("ClonalityScore_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
-    }
-    
-    clonalityScore = clonalFreqCount[c("Sample", "Result")]
-    clonalityScore = unique(clonalityScore)
-    
-    clonalityScoreSplit = split(clonalityScore, f=clonalityScore[,"Sample"])
-    lapply(clonalityScoreSplit, FUN=ClonalityScorePrint)
-    
-    clonalityOverview = clonalFreqCount[c("Sample", "Type", "Count", "Weight", "WeightedCount")]
-    
-    
-    
-    ClonalityOverviewPrint <- function(dat){
-	  dat = dat[order(dat[,2]),]
-      write.table(dat[-1], paste("ClonalityOverView_", unique(dat[1])[1,1] , ".csv", sep=""), sep=",",quote=F,na="-",row.names=F,col.names=F)
-    }
-    
-    clonalityOverviewSplit = split(clonalityOverview, f=clonalityOverview$Sample)
-    lapply(clonalityOverviewSplit, FUN=ClonalityOverviewPrint)
+  }
+  clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "clonaltype")])
+
+  #write files for every coincidence group of >1
+  samples = unique(clonalFreq$Sample)
+  for(sample in samples){
+      clonalFreqSample = clonalFreq[clonalFreq$Sample == sample,]
+      if(max(clonalFreqSample$Type) > 1){
+          for(i in 2:max(clonalFreqSample$Type)){
+              clonalFreqSampleType = clonalFreqSample[clonalFreqSample$Type == i,]
+              clonalityFrame.sub = clonalityFrame[clonalityFrame$clonaltype %in% clonalFreqSampleType$clonaltype,]
+              clonalityFrame.sub = clonalityFrame.sub[order(clonalityFrame.sub$clonaltype),]
+              write.table(clonalityFrame.sub, file=paste("coincidences_", sample, "_", i, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
+          }
+      }
+  }
+
+  clonalFreqCount = data.frame(data.table(clonalFreq)[, list(Count=.N), by=c("Sample", "Type")])
+  clonalFreqCount$realCount = clonalFreqCount$Type * clonalFreqCount$Count
+  clonalSum = data.frame(data.table(clonalFreqCount)[, list(Reads=sum(realCount)), by=c("Sample")])
+  clonalFreqCount = merge(clonalFreqCount, clonalSum, by.x="Sample", by.y="Sample")
+
+  ct = c('Type\tWeight\n2\t1\n3\t3\n4\t6\n5\t10\n6\t15')
+  tcct = textConnection(ct)
+  CT  = read.table(tcct, sep="\t", header=TRUE)
+  close(tcct)
+  clonalFreqCount = merge(clonalFreqCount, CT, by.x="Type", by.y="Type", all.x=T)
+  clonalFreqCount$WeightedCount = clonalFreqCount$Count * clonalFreqCount$Weight
+
+  ReplicateReads = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "Replicate", "clonaltype")])
+  ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")])
+  clonalFreqCount$Reads = as.numeric(clonalFreqCount$Reads)
+  ReplicateReads$Reads = as.numeric(ReplicateReads$Reads)
+  ReplicateReads$squared = as.numeric(ReplicateReads$Reads * ReplicateReads$Reads)
+
+  ReplicatePrint <- function(dat){
+      write.table(dat[-1], paste("ReplicateReads_", unique(dat[1])[1,1] , ".txt", sep=""), sep="\t",quote=F,na="-",row.names=F,col.names=F)
   }
+
+  ReplicateSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"])
+  lapply(ReplicateSplit, FUN=ReplicatePrint)
+
+  ReplicateReads = data.frame(data.table(ReplicateReads)[, list(ReadsSum=sum(as.numeric(Reads)), ReadsSquaredSum=sum(as.numeric(squared))), by=c("Sample")])
+  clonalFreqCount = merge(clonalFreqCount, ReplicateReads, by.x="Sample", by.y="Sample", all.x=T)
+
+  ReplicateSumPrint <- function(dat){
+      write.table(dat[-1], paste("ReplicateSumReads_", unique(dat[1])[1,1] , ".txt", sep=""), sep="\t",quote=F,na="-",row.names=F,col.names=F)
+  }
+
+  ReplicateSumSplit = split(ReplicateReads, f=ReplicateReads[,"Sample"])
+  lapply(ReplicateSumSplit, FUN=ReplicateSumPrint)
+
+  clonalFreqCountSum = data.frame(data.table(clonalFreqCount)[, list(Numerator=sum(WeightedCount, na.rm=T)), by=c("Sample")])
+  clonalFreqCount = merge(clonalFreqCount, clonalFreqCountSum, by.x="Sample", by.y="Sample", all.x=T)
+  clonalFreqCount$ReadsSum = as.numeric(clonalFreqCount$ReadsSum) #prevent integer overflow
+  clonalFreqCount$Denominator = (((clonalFreqCount$ReadsSum * clonalFreqCount$ReadsSum) - clonalFreqCount$ReadsSquaredSum) / 2)
+  clonalFreqCount$Result = (clonalFreqCount$Numerator + 1) / (clonalFreqCount$Denominator + 1)
+
+  ClonalityScorePrint <- function(dat){
+      write.table(dat$Result, paste("ClonalityScore_", unique(dat[1])[1,1] , ".txt", sep=""), sep="\t",quote=F,na="-",row.names=F,col.names=F)
+  }
+
+  clonalityScore = clonalFreqCount[c("Sample", "Result")]
+  clonalityScore = unique(clonalityScore)
+
+  clonalityScoreSplit = split(clonalityScore, f=clonalityScore[,"Sample"])
+  lapply(clonalityScoreSplit, FUN=ClonalityScorePrint)
+
+  clonalityOverview = clonalFreqCount[c("Sample", "Type", "Count", "Weight", "WeightedCount")]
+
+
+
+  ClonalityOverviewPrint <- function(dat){
+      dat = dat[order(dat[,2]),]
+      write.table(dat[-1], paste("ClonalityOverView_", unique(dat[1])[1,1] , ".txt", sep=""), sep="\t",quote=F,na="-",row.names=F,col.names=F)
+  }
+
+  clonalityOverviewSplit = split(clonalityOverview, f=clonalityOverview$Sample)
+  lapply(clonalityOverviewSplit, FUN=ClonalityOverviewPrint)
+  
 }
 
 bak = PRODF
@@ -724,16 +724,14 @@
   
   PRODF.with.D = PRODF[nchar(PRODF$Top.D.Gene, keepNA=F) > 2,]
   PRODF.no.D = PRODF[nchar(PRODF$Top.D.Gene, keepNA=F) < 4,]
+  write.table(PRODF.no.D, "productive_no_D.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=T)
   
   UNPROD.with.D = UNPROD[nchar(UNPROD$Top.D.Gene, keepNA=F) > 2,]
   UNPROD.no.D = UNPROD[nchar(UNPROD$Top.D.Gene, keepNA=F) < 4,]
+  write.table(UNPROD.no.D, "unproductive_no_D.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=T)
   
   num_median = function(x, na.rm=T) { as.numeric(median(x, na.rm=na.rm)) }
 
-  print("---- table prod.with.d cdr3.length ----")
-  print(table(PRODF.with.D$CDR3.Length, useNA="ifany"))
-  print(median(PRODF.with.D$CDR3.Length, na.rm=T))
-
   newData = data.frame(data.table(PRODF.with.D)[,list(unique=.N, 
                                                VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),
                                                P1=mean(.SD$P3V.nt.nb, na.rm=T),
@@ -772,10 +770,6 @@
   newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1)
   write.table(newData, "junctionAnalysisProd_median_wD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F)
   
-  print("---- table unprod.with.d cdr3.length ----")
-  print(table(UNPROD.with.D$CDR3.Length, useNA="ifany"))
-  print(median(UNPROD.with.D$CDR3.Length, na.rm=T))
-  
   newData = data.frame(data.table(UNPROD.with.D)[,list(unique=.N, 
                                                 VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),
                                                 P1=mean(.SD$P3V.nt.nb, na.rm=T),
@@ -819,11 +813,11 @@
   newData = data.frame(data.table(PRODF.no.D)[,list(unique=.N, 
                                                VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),
                                                P1=mean(.SD$P3V.nt.nb, na.rm=T),
-                                               N1=mean(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)),
+                                               N1=mean(.SD$N.REGION.nt.nb, na.rm=T),
                                                P2=mean(.SD$P5J.nt.nb, na.rm=T),
                                                DEL.JH=mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T),
                                                Total.Del=mean(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb"), with=F], na.rm=T)),
-                                               Total.N=mean(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)),
+                                               Total.N=mean(.SD$N.REGION.nt.nb, na.rm=T),
                                                Total.P=mean(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
                                                Median.CDR3.l=as.double(as.numeric(median(.SD$CDR3.Length, na.rm=T)))),
                                          by=c("Sample")])
@@ -833,39 +827,43 @@
   newData = data.frame(data.table(PRODF.no.D)[,list(unique=.N, 
                                                VH.DEL=num_median(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),
                                                P1=num_median(.SD$P3V.nt.nb, na.rm=T),
-                                               N1=num_median(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)),
+                                               N1=median(.SD$N.REGION.nt.nb, na.rm=T),
                                                P2=num_median(.SD$P5J.nt.nb, na.rm=T),
                                                DEL.JH=num_median(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T),
 											   Total.Del=num_median(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb"), with=F], na.rm=T)),
-											   Total.N=num_median(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)),
+											   Total.N=median(.SD$N.REGION.nt.nb, na.rm=T),
 											   Total.P=num_median(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
 											   Median.CDR3.l=as.double(as.numeric(median(.SD$CDR3.Length, na.rm=T)))),
                                          by=c("Sample")])
   newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1)
   write.table(newData, "junctionAnalysisProd_median_nD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F)
   
+  print(paste("mean N:", mean(UNPROD.no.D$N.REGION.nt.nb, na.rm=T)))
+  
   newData = data.frame(data.table(UNPROD.no.D)[,list(unique=.N, 
                                                 VH.DEL=mean(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),
                                                 P1=mean(.SD$P3V.nt.nb, na.rm=T),
-                                                N1=mean(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)),
+                                                N1=mean(.SD$N.REGION.nt.nb, na.rm=T),
                                                 P2=mean(.SD$P5J.nt.nb, na.rm=T),
                                                 DEL.JH=mean(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T),
                                                 Total.Del=mean(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb"), with=F], na.rm=T)),
-                                                Total.N=mean(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)),
+                                                Total.N=mean(.SD$N.REGION.nt.nb, na.rm=T),
                                                 Total.P=mean(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
                                                 Median.CDR3.l=as.double(as.numeric(median(.SD$CDR3.Length, na.rm=T)))),
                                           by=c("Sample")])
   newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1)
   write.table(newData, "junctionAnalysisUnProd_mean_nD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F)
   
+  print(paste("median N:", num_median(UNPROD.no.D$N.REGION.nt.nb, na.rm=T)))
+  
     newData = data.frame(data.table(UNPROD.no.D)[,list(unique=.N, 
                                                 VH.DEL=num_median(.SD$X3V.REGION.trimmed.nt.nb, na.rm=T),
                                                 P1=num_median(.SD$P3V.nt.nb, na.rm=T),
-                                                N1=num_median(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)),
+                                                N1=median(.SD$N.REGION.nt.nb, na.rm=T),
                                                 P2=num_median(.SD$P5J.nt.nb, na.rm=T),
                                                 DEL.JH=num_median(.SD$X5J.REGION.trimmed.nt.nb, na.rm=T),
                                                 Total.Del=num_median(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb"), with=F], na.rm=T)),
-                                                Total.N=num_median(rowSums(.SD[,c("N.REGION.nt.nb"), with=F], na.rm=T)),
+                                                Total.N=median(.SD$N.REGION.nt.nb, na.rm=T),
                                                 Total.P=num_median(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
                                                 Median.CDR3.l=as.double(as.numeric(median(.SD$CDR3.Length, na.rm=T)))),
 															by=c("Sample")])
@@ -960,20 +958,13 @@
 clonaltype.in.replicates$clonaltype = do.call(paste, c(clonaltype.in.replicates[clonaltype], sep = ":"))
 clonaltype.in.replicates = clonaltype.in.replicates[,c("clonaltype","Replicate", "ID", "Sequence", "Sample")]
 
-print(head(clonaltype.in.replicates))
-
 clonaltype.counts = data.frame(table(clonaltype.in.replicates$clonaltype))
 names(clonaltype.counts) = c("clonaltype", "coincidence")
 
-print(head(clonaltype.counts))
-
 clonaltype.counts = clonaltype.counts[clonaltype.counts$coincidence > 1,]
 
-head(clonaltype.counts)
-
 clonaltype.in.replicates = clonaltype.in.replicates[clonaltype.in.replicates$clonaltype %in% clonaltype.counts$clonaltype,]
 clonaltype.in.replicates = merge(clonaltype.in.replicates, clonaltype.counts, by="clonaltype")
-print(head(clonaltype.in.replicates))
 clonaltype.in.replicates = clonaltype.in.replicates[order(clonaltype.in.replicates$clonaltype),c("coincidence","clonaltype", "Sample", "Replicate", "ID", "Sequence")]
 
 write.table(clonaltype.in.replicates, "clonaltypes_replicates.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=T)