Mercurial > repos > davidvanzessen > argalaxy_tools

changeset 25:94765af0db1f draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Thu, 09 Feb 2017 07:20:09 -0500
parents d5d203d38c8a
children 28fbbdfd7a87
files report_clonality/RScript.r report_clonality/r_wrapper.sh
diffstat 2 files changed, 79 insertions(+), 15 deletions(-) [+]
line wrap: on
line diff
--- a/report_clonality/RScript.r	Wed Feb 01 09:48:38 2017 -0500
+++ b/report_clonality/RScript.r	Thu Feb 09 07:20:09 2017 -0500
@@ -158,7 +158,7 @@
 #write the complete dataset that is left over, will be the input if 'none' for clonaltype and 'no' for filterproductive
 write.table(PRODF, "allUnique.txt", sep="\t",quote=F,row.names=F,col.names=T)
 #write.table(PRODF, "allUnique.csv", sep=",",quote=F,row.names=F,col.names=T)
-write.table(UNPROD, "allUnproductive.csv", sep=",",quote=F,row.names=F,col.names=T)
+write.table(UNPROD, "allUnproductive.txt", sep="\t",quote=F,row.names=F,col.names=T)
 
 print("SAMPLE TABLE:")
 print(table(PRODF$Sample))
@@ -697,6 +697,7 @@
 }
 
 bak = PRODF
+bakun = UNPROD
 
 imgtcolumns = c("X3V.REGION.trimmed.nt.nb","P3V.nt.nb", "N1.REGION.nt.nb", "P5D.nt.nb", "X5D.REGION.trimmed.nt.nb", "X3D.REGION.trimmed.nt.nb", "P3D.nt.nb", "N2.REGION.nt.nb", "P5J.nt.nb", "X5J.REGION.trimmed.nt.nb", "X3V.REGION.trimmed.nt.nb", "X5D.REGION.trimmed.nt.nb", "X3D.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb", "N1.REGION.nt.nb", "N2.REGION.nt.nb", "P3V.nt.nb", "P5D.nt.nb", "P3D.nt.nb", "P5J.nt.nb")
 if(all(imgtcolumns %in% colnames(inputdata)))
@@ -724,8 +725,15 @@
   PRODF.with.D = PRODF[nchar(PRODF$Top.D.Gene, keepNA=F) > 2,]
   PRODF.no.D = PRODF[nchar(PRODF$Top.D.Gene, keepNA=F) < 4,]
   
+  UNPROD.with.D = UNPROD[nchar(UNPROD$Top.D.Gene, keepNA=F) > 2,]
+  UNPROD.no.D = UNPROD[nchar(UNPROD$Top.D.Gene, keepNA=F) < 4,]
+  
   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),
@@ -740,7 +748,7 @@
                                                Total.Del=mean(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5D.REGION.trimmed.nt.nb", "X3D.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb"), with=F], na.rm=T)),
                                                Total.N=mean(rowSums(.SD[,c("N.REGION.nt.nb", "N1.REGION.nt.nb", "N2.REGION.nt.nb", "N3.REGION.nt.nb", "N4.REGION.nt.nb"), with=F], na.rm=T)),
                                                Total.P=mean(rowSums(.SD[,c("P3V.nt.nb", "P5D.nt.nb", "P3D.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
-                                               Median.CDR3.l=as.double(median(.SD$CDR3.Length))),
+                                               Median.CDR3.l=as.double(median(as.numeric(.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_mean_wD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F)
@@ -759,12 +767,16 @@
 											   Total.Del=num_median(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5D.REGION.trimmed.nt.nb", "X3D.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", "N1.REGION.nt.nb", "N2.REGION.nt.nb", "N3.REGION.nt.nb", "N4.REGION.nt.nb"), with=F], na.rm=T)),
 											   Total.P=num_median(rowSums(.SD[,c("P3V.nt.nb", "P5D.nt.nb", "P3D.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
-											   Median.CDR3.l=as.double(median(.SD$CDR3.Length))),
+											   Median.CDR3.l=as.double(median(as.numeric(.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_wD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F)
   
-  newData = data.frame(data.table(PRODF.with.D)[,list(unique=.N, 
+  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),
                                                 N1=mean(rowSums(.SD[,c("N.REGION.nt.nb", "N1.REGION.nt.nb"), with=F], na.rm=T)),
@@ -778,12 +790,12 @@
                                                 Total.Del=mean(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5D.REGION.trimmed.nt.nb", "X3D.REGION.trimmed.nt.nb", "X5J.REGION.trimmed.nt.nb"), with=F], na.rm=T)),
                                                 Total.N=mean(rowSums(.SD[,c("N.REGION.nt.nb", "N1.REGION.nt.nb", "N2.REGION.nt.nb", "N3.REGION.nt.nb", "N4.REGION.nt.nb"), with=F], na.rm=T)),
                                                 Total.P=mean(rowSums(.SD[,c("P3V.nt.nb", "P5D.nt.nb", "P3D.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
-                                                Median.CDR3.l=as.double(median(.SD$CDR3.Length))),
+                                                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_wD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F)
   
-    newData = data.frame(data.table(PRODF.with.D)[,list(unique=.N, 
+    newData = data.frame(data.table(UNPROD.with.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", "N1.REGION.nt.nb"), with=F], na.rm=T)),
@@ -797,7 +809,7 @@
                                                 Total.Del=num_median(rowSums(.SD[,c("X3V.REGION.trimmed.nt.nb", "X5D.REGION.trimmed.nt.nb", "X3D.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", "N1.REGION.nt.nb", "N2.REGION.nt.nb", "N3.REGION.nt.nb", "N4.REGION.nt.nb"), with=F], na.rm=T)),
                                                 Total.P=num_median(rowSums(.SD[,c("P3V.nt.nb", "P5D.nt.nb", "P3D.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
-                                                Median.CDR3.l=as.double(median(.SD$CDR3.Length))),
+                                                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_median_wD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F)
@@ -813,9 +825,8 @@
                                                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.P=mean(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
-                                               Median.CDR3.l=as.double(median(.SD$CDR3.Length))),
+                                               Median.CDR3.l=as.double(as.numeric(median(.SD$CDR3.Length, na.rm=T)))),
                                          by=c("Sample")])
-  print(newData)
   newData[,sapply(newData, is.numeric)] = round(newData[,sapply(newData, is.numeric)],1)
   write.table(newData, "junctionAnalysisProd_mean_nD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F)
   
@@ -828,12 +839,12 @@
 											   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.P=num_median(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
-											   Median.CDR3.l=as.double(median(.SD$CDR3.Length))),
+											   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)
   
-  newData = data.frame(data.table(PRODF.no.D)[,list(unique=.N, 
+  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)),
@@ -842,12 +853,12 @@
                                                 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.P=mean(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
-                                                Median.CDR3.l=as.double(median(.SD$CDR3.Length))),
+                                                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)
   
-    newData = data.frame(data.table(PRODF.no.D)[,list(unique=.N, 
+    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)),
@@ -856,13 +867,14 @@
                                                 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.P=num_median(rowSums(.SD[,c("P3V.nt.nb", "P5J.nt.nb"), with=F], na.rm=T)),
-                                                Median.CDR3.l=as.double(median(.SD$CDR3.Length))),
+                                                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_median_nD.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F)
 }
 
 PRODF = bak
+UNPROD = bakun
 
 
 # ---------------------- D reading frame ----------------------
@@ -939,3 +951,54 @@
 median.aa.l = data.frame(data.table(PRODF)[, list(median=as.double(median(as.numeric(.SD$CDR3.Length, na.rm=T), na.rm=T))), by=c("Sample")])
 write.table(median.aa.l, "AAMedianBySample.txt" , sep="\t",quote=F,na="-",row.names=F,col.names=F)
 
+
+#generate the "Sequences that are present in more than one replicate" dataset
+clonaltype.in.replicates = inputdata
+clonaltype = unlist(strsplit(clonaltype, ","))
+clonaltype = clonaltype[-which(clonaltype == "Sample")]
+
+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)
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
--- a/report_clonality/r_wrapper.sh	Wed Feb 01 09:48:38 2017 -0500
+++ b/report_clonality/r_wrapper.sh	Thu Feb 09 07:20:09 2017 -0500
@@ -396,6 +396,7 @@
 
 echo "<tr><td colspan='2' style='background-color:#E0E0E0;'>Clonality</td></tr>" >> $outputFile
 echo "<tr><td>The dataset used to calculate clonality score (Unique based on clonaltype, $clonalType)</td><td><a href='clonalityComplete.txt'>Download</a></td></tr>" >> $outputFile
+echo "<tr><td>Sequences that are present in more than one replicate</td><td><a href='clonaltypes_replicates.txt'>Download</a></td></tr>" >> $outputFile
 
 echo "</table>" >> $outputFile