diff report_clonality/RScript.r @ 47:d97e1421aa86 draft

Uploaded

--- a/report_clonality/RScript.r	Wed Jan 27 08:39:16 2016 -0500
+++ b/report_clonality/RScript.r	Wed Jan 27 10:25:43 2016 -0500
@@ -41,6 +41,7 @@
 locus = args[6] # IGH, IGK, IGL, TRB, TRA, TRG or TRD
 filterproductive = ifelse(args[7] == "yes", T, F) #should unproductive sequences be filtered out? (yes/no)
 clonality_method = args[8]
+filter_uniques = args[9]
 
 # ---------------------- Data preperation ----------------------
 
@@ -58,6 +59,28 @@
 inputdata$Top.D.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.D.Gene)
 inputdata$Top.J.Gene = gsub("[*]([0-9]+)", "", inputdata$Top.J.Gene)
 
+#filter uniques
+inputdata.removed = inputdata[NULL,]
+
+if(filter_uniques == "yes" && c("CDR1.Seq", "CDR2.Seq", "CDR3.Seq", "FR1.IMGT", "FR2.IMGT", "FR3.IMGT") %in% names(inputdata)){
+	
+	clmns = names(inputdata)
+	
+	inputdata$unique.def = paste(inputdata$CDR1.Seq, inputdata$CDR2.Seq, inputdata$CDR3.Seq, inputdata$FR1.IMGT, inputdata$FR2.IMGT, inputdata$FR3.IMGT)
+	inputdata.filtered = inputdata[duplicated(inputdata$unique.def),]
+	fltr = inputdata$unique.def %in% inputdata.filtered$unique.def
+	
+	inputdata.removed = inputdata[!fltr,]
+	inputdata.removed$samples_replicates = paste(inputdata.removed$Sample, inputdata.removed$Replicate, sep="_")
+	
+	inputdata = inputdata[fltr,]
+	
+	inputdata = inputdata[,clmns]
+	
+	write.table(inputdata.removed, "unique_removed.csv", sep=",",quote=F,row.names=F,col.names=T)
+}
+
+
 inputdata$clonaltype = 1:nrow(inputdata)
 
 PRODF = inputdata
@@ -154,6 +177,12 @@
 sample_productive_count$perc_unprod = round(sample_productive_count$Unproductive / sample_productive_count$All * 100)
 sample_productive_count$perc_unprod_un = round(sample_productive_count$Unproductive_unique / sample_productive_count$All * 100)
 
+inputdata.removed.s = data.table(inputdata.removed)[, list(UniqueRemoved=.N), by=c("Sample")]
+
+sample_productive_count = merge(sample_productive_count, inputdata.removed.s, by="Sample")
+
+sample_productive_count$perc_rem = round(sample_productive_count$UniqueRemoved / sample_productive_count$All * 100)
+
 
 sample_replicate_productive_count = inputdata.dt[, list(All=.N, 
                                                         Productive = nrow(.SD[.SD$Functionality == "productive" | .SD$Functionality == "productive (see comment)",]), 
@@ -172,6 +201,13 @@
 sample_replicate_productive_count$perc_unprod = round(sample_replicate_productive_count$Unproductive / sample_replicate_productive_count$All * 100)
 sample_replicate_productive_count$perc_unprod_un = round(sample_replicate_productive_count$Unproductive_unique / sample_replicate_productive_count$All * 100)
 
+inputdata.removed.sr = data.table(inputdata.removed)[, list(UniqueRemoved=.N), by=c("samples_replicates")]
+
+sample_replicate_productive_count = merge(sample_replicate_productive_count, inputdata.removed.sr, by="samples_replicates")
+
+sample_replicate_productive_count$perc_rem = round(sample_replicate_productive_count$UniqueRemoved / sample_productive_count$All * 100)
+
+
 setnames(sample_replicate_productive_count, colnames(sample_productive_count))
 
 counts = rbind(sample_replicate_productive_count, sample_productive_count)