--- a/RScript.r	Mon Nov 25 09:21:55 2013 -0500
+++ b/RScript.r	Tue Dec 10 05:53:08 2013 -0500
@@ -1,4 +1,4 @@
-options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
+#options( show.error.messages=F, error = function () { cat( geterrmessage(), file=stderr() ); q( "no", 1, F ) } )
 
 args <- commandArgs(trailingOnly = TRUE)
 
@@ -30,11 +30,6 @@
 
 test = test[test$Sample != "",]
 
-if("Replicate" %in% colnames(test))
-{
-	test$SRID = do.call(paste, c(test[c("Sample", "Replicate")], sep = "-"))
-}
-
 test$Top.V.Gene = gsub("[*]([0-9]+)", "", test$Top.V.Gene)
 test$Top.D.Gene = gsub("[*]([0-9]+)", "", test$Top.D.Gene)
 test$Top.J.Gene = gsub("[*]([0-9]+)", "", test$Top.J.Gene)
@@ -43,6 +38,9 @@
 test$VDJCDR3 = do.call(paste, c(test[unlist(strsplit(clonalType, ","))], sep = ":"))
 
 PROD = test[test$VDJ.Frame != "In-frame with stop codon" & test$VDJ.Frame != "Out-of-frame" & test$CDR3.Found.How != "NOT_FOUND" , ]
+if("Functionality" %in% colnames(test)) {
+	PROD = test[test$Functionality == "productive" | test$Functionality == "productive (see comment)", ]
+}
 
 NONPROD = test[test$VDJ.Frame == "In-frame with stop codon" | test$VDJ.Frame == "Out-of-frame" | test$CDR3.Found.How == "NOT_FOUND" , ]
 
@@ -74,7 +72,7 @@
 PRODFJ = merge(PRODFJ, Total, by.x='Sample', by.y='Sample', all.x=TRUE)
 PRODFJ = ddply(PRODFJ, c("Sample", "Top.J.Gene"), summarise, relFreq= (Length*100 / Total))
 
-V = c("v.name\tchr.orderV\nIGHV7-81\t1\nIGHV3-74\t2\nIGHV3-73\t3\nIGHV3-72\t4\nIGHV3-71\t5\nIGHV2-70\t6\nIGHV1-69\t7\nIGHV3-66\t8\nIGHV3-64\t9\nIGHV4-61\t10\nIGHV4-59\t11\nIGHV1-58\t12\nIGHV3-53\t13\nIGHV3-52\t14\nIGHV5-a\t15\nIGHV5-51\t16\nIGHV3-49\t17\nIGHV3-48\t18\nIGHV3-47\t19\nIGHV1-46\t20\nIGHV1-45\t21\nIGHV3-43\t22\nIGHV4-39\t23\nIGHV3-35\t24\nIGHV4-34\t25\nIGHV3-33\t26\nIGHV4-31\t27\nIGHV4-30-4\t28\nIGHV4-30-2\t29\nIGHV3-30-3\t30\nIGHV3-30\t31\nIGHV4-28\t32\nIGHV2-26\t33\nIGHV1-24\t34\nIGHV3-23\t35\nIGHV3-22\t36\nIGHV3-21\t37\nIGHV3-20\t38\nIGHV3-19\t39\nIGHV1-18\t40\nIGHV3-15\t41\nIGHV3-13\t42\nIGHV3-11\t43\nIGHV3-9\t44\nIGHV1-8\t45\nIGHV3-7\t46\nIGHV2-5\t47\nIGHV7-4-1\t48\nIGHV4-4\t49\nIGHV4-b\t50\nIGHV1-3\t51\nIGHV1-2\t52\nIGHV6-1\t53")
+V = c("v.name\tchr.orderV\nIGHV7-81\t1\nIGHV3-74\t2\nIGHV3-73\t3\nIGHV3-72\t4\nIGHV2-70\t6\nIGHV1-69\t7\nIGHV3-66\t8\nIGHV3-64\t9\nIGHV4-61\t10\nIGHV4-59\t11\nIGHV1-58\t12\nIGHV3-53\t13\nIGHV5-a\t15\nIGHV5-51\t16\nIGHV3-49\t17\nIGHV3-48\t18\nIGHV1-46\t20\nIGHV1-45\t21\nIGHV3-43\t22\nIGHV4-39\t23\nIGHV3-35\t24\nIGHV4-34\t25\nIGHV3-33\t26\nIGHV4-31\t27\nIGHV4-30-4\t28\nIGHV4-30-2\t29\nIGHV3-30-3\t30\nIGHV3-30\t31\nIGHV4-28\t32\nIGHV2-26\t33\nIGHV1-24\t34\nIGHV3-23\t35\nIGHV3-21\t37\nIGHV3-20\t38\nIGHV1-18\t40\nIGHV3-15\t41\nIGHV3-13\t42\nIGHV3-11\t43\nIGHV3-9\t44\nIGHV1-8\t45\nIGHV3-7\t46\nIGHV2-5\t47\nIGHV7-4-1\t48\nIGHV4-4\t49\nIGHV4-b\t50\nIGHV1-3\t51\nIGHV1-2\t52\nIGHV6-1\t53")
 tcV = textConnection(V)
 Vchain = read.table(tcV, sep="\t", header=TRUE)
 PRODFV = merge(PRODFV, Vchain, by.x='Top.V.Gene', by.y='v.name', all.x=TRUE)
@@ -95,6 +93,8 @@
 
 setwd(outDir)
 
+write.table(PRODF, "allUnique.tsv", sep="\t",quote=F,row.names=T,col.names=T)
+
 pV = ggplot(PRODFV)
 pV = pV + geom_bar( aes( x=factor(reorder(Top.V.Gene, chr.orderV)), y=relFreq, fill=Sample), stat='identity', position="dodge") + theme(axis.text.x = element_text(angle = 90, hjust = 1))
 pV = pV + xlab("Summary of V gene") + ylab("Frequency") + ggtitle("Relative frequency of V gene usage")
@@ -131,7 +131,7 @@
 	img = ggplot() + 
 	geom_tile(data=dat, aes(x=factor(reorder(Top.D.Gene, chr.orderD)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) + 
 	theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
-	scale_fill_gradient(low="gold", high="blue", na.value="white", limits=c(0,1)) + 
+	scale_fill_gradient(low="gold", high="blue", na.value="white") + 
 	ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 
 	xlab("D genes") + 
 	ylab("V Genes")
@@ -166,7 +166,7 @@
 	img = ggplot() + 
 	geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.V.Gene, chr.orderV)), fill=relLength)) + 
 	theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
-	scale_fill_gradient(low="gold", high="blue", na.value="white", limits=c(0,1)) + 
+	scale_fill_gradient(low="gold", high="blue", na.value="white") + 
 	ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 
 	xlab("J genes") + 
 	ylab("V Genes")
@@ -198,7 +198,7 @@
 	img = ggplot() + 
 	geom_tile(data=dat, aes(x=factor(reorder(Top.J.Gene, chr.orderJ)), y=factor(reorder(Top.D.Gene, chr.orderD)), fill=relLength)) + 
 	theme(axis.text.x = element_text(angle = 90, hjust = 1)) + 
-	scale_fill_gradient(low="gold", high="blue", na.value="white", limits=c(0,1)) + 
+	scale_fill_gradient(low="gold", high="blue", na.value="white") + 
 	ggtitle(paste(unique(dat$Sample), " (N=" , sum(dat$Length, na.rm=T) ,")", sep="")) + 
 	xlab("J genes") + 
 	ylab("D Genes")
@@ -229,3 +229,88 @@
 un = paste(un, sep="\n")
 writeLines(un, sampleFile)
 close(sampleFile)
+
+
+if("Replicate" %in% colnames(test))
+{
+	clonalityFrame = PROD
+	clonalityFrame$ReplicateConcat = do.call(paste, c(clonalityFrame[c("VDJCDR3", "Sample", "Replicate")], sep = ":"))
+	clonalityFrame = clonalityFrame[!duplicated(clonalityFrame$ReplicateConcat), ]
+	write.table(clonalityFrame, "clonalityComplete.tsv", sep="\t",quote=F,row.names=T,col.names=T)
+
+	ClonalitySampleReplicatePrint <- function(dat){
+	    write.table(dat, paste("clonality_", unique(dat$Sample) , "_", unique(dat$Replicate), ".tsv", sep=""), sep="\t",quote=F,row.names=T,col.names=T)
+	}
+
+    clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,c("Sample", "Replicate")])
+    lapply(clonalityFrameSplit, FUN=ClonalitySampleReplicatePrint)
+
+    ClonalitySamplePrint <- function(dat){
+	    write.table(dat, paste("clonality_", unique(dat$Sample) , ".tsv", sep=""), sep="\t",quote=F,row.names=T,col.names=T)
+	}
+
+    clonalityFrameSplit = split(clonalityFrame, f=clonalityFrame[,"Sample"])
+    lapply(clonalityFrameSplit, FUN=ClonalitySamplePrint)
+
+	clonalFreq = data.frame(data.table(clonalityFrame)[, list(Type=.N), by=c("Sample", "VDJCDR3")])
+	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", "VDJCDR3")])
+	ReplicateReads = data.frame(data.table(ReplicateReads)[, list(Reads=.N), by=c("Sample", "Replicate")])
+	ReplicateReads$squared = 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(Reads), ReadsSquaredSum=sum(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$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){
+		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)
+}