--- a/RScript.r	Mon Jan 12 04:32:40 2015 -0500
+++ b/RScript.r	Mon Jan 12 05:35:42 2015 -0500
@@ -50,11 +50,14 @@
 inputdata$Top.J.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.J.Gene)
 inputdata$clonaltype = 1:nrow(inputdata)
 PRODF = inputdata
+UNPROD = inputdata
 if(filterproductive){
   if("Functionality" %in% colnames(inputdata)) { # "Functionality" is an IMGT column
     PRODF = inputdata[inputdata$Functionality == "productive" | inputdata$Functionality == "productive (see comment)", ]
+    UNPROD = inputdata[!(inputdata$Functionality == "productive" | inputdata$Functionality == "productive (see comment)"), ]
   } else {
     PRODF = inputdata[inputdata$VDJ.Frame != "In-frame with stop codon" & inputdata$VDJ.Frame != "Out-of-frame" & inputdata$CDR3.Found.How != "NOT_FOUND" , ]
+    UNPROD = inputdata[!(inputdata$VDJ.Frame != "In-frame with stop codon" & inputdata$VDJ.Frame != "Out-of-frame" & inputdata$CDR3.Found.How != "NOT_FOUND" ), ]
   }
 }
 
@@ -62,6 +65,8 @@
 if(clonaltype != "none"){
   PRODF$clonaltype = do.call(paste, c(PRODF[unlist(strsplit(clonaltype, ","))], sep = ":"))
   PRODF = PRODF[!duplicated(PRODF$clonaltype), ]
+  UNPROD$clonaltype = do.call(paste, c(UNPROD[unlist(strsplit(clonaltype, ","))], sep = ":"))
+  UNPROD = UNPROD[!duplicated(UNPROD$clonaltype), ]
 }
 
 PRODF$freq = 1
@@ -576,5 +581,31 @@
                                                           mean(P3D.nt.nb, na.rm=T) +
                                                           mean(P5J.nt.nb, na.rm=T))),
                                         by=c("Sample")])
-  write.table(newData, "junctionAnalysis.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F)
+  write.table(newData, "junctionAnalysisProd.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F)
+  
+	newData = data.frame(data.table(UNPROD)[,list(unique=.N, 
+                                              VH.DEL=mean(X3V.REGION.trimmed.nt.nb, na.rm=T),
+                                              P1=mean(P3V.nt.nb, na.rm=T),
+                                              N1=mean(N1.REGION.nt.nb, na.rm=T),
+                                              P2=mean(P5D.nt.nb, na.rm=T),
+                                              DEL.DH=mean(X5D.REGION.trimmed.nt.nb, na.rm=T),
+                                              DH.DEL=mean(X3D.REGION.trimmed.nt.nb, na.rm=T),
+                                              P3=mean(P3D.nt.nb, na.rm=T),
+                                              N2=mean(N2.REGION.nt.nb, na.rm=T),
+                                              P4=mean(P5J.nt.nb, na.rm=T),
+                                              DEL.JH=mean(X5J.REGION.trimmed.nt.nb, na.rm=T),
+                                              Total.Del=(	mean(X3V.REGION.trimmed.nt.nb, na.rm=T) + 
+                                                            mean(X5D.REGION.trimmed.nt.nb, na.rm=T) + 
+                                                            mean(X3D.REGION.trimmed.nt.nb, na.rm=T) +
+                                                            mean(X5J.REGION.trimmed.nt.nb, na.rm=T)),
+                                              
+                                              Total.N=(	mean(N1.REGION.nt.nb, na.rm=T) +
+                                                          mean(N2.REGION.nt.nb, na.rm=T)),
+                                              
+                                              Total.P=(	mean(P3V.nt.nb, na.rm=T) +
+                                                          mean(P5D.nt.nb, na.rm=T) +
+                                                          mean(P3D.nt.nb, na.rm=T) +
+                                                          mean(P5J.nt.nb, na.rm=T))),
+                                        by=c("Sample")])
+  write.table(newData, "junctionAnalysisUnProd.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F)
 }