changeset 5:fb219dfdcc15 draft

Uploaded
author davidvanzessen
date Thu, 21 Aug 2014 10:44:40 -0400
parents 0748ef4d42d3
children f4ff7450ef16
files RScript.r
diffstat 1 files changed, 22 insertions(+), 6 deletions(-) [+]
line wrap: on
line diff
--- a/RScript.r	Thu May 15 10:00:20 2014 -0400
+++ b/RScript.r	Thu Aug 21 10:44:40 2014 -0400
@@ -17,11 +17,11 @@
 if (!("ggplot2" %in% rownames(installed.packages()))) {
 	install.packages("ggplot2", repos="http://cran.xl-mirror.nl/") 
 }
-require(ggplot2)
+library(ggplot2)
 if (!("plyr" %in% rownames(installed.packages()))) {
 	install.packages("plyr", repos="http://cran.xl-mirror.nl/") 
 }			
-require(plyr)
+library(plyr)
 
 if (!("data.table" %in% rownames(installed.packages()))) {
 	install.packages("data.table", repos="http://cran.xl-mirror.nl/") 
@@ -38,7 +38,7 @@
 
 test = test[test$Sample != "",]
 
-
+print("test1\n")
 
 test$Top.V.Gene = gsub("[*]([0-9]+)", "", test$Top.V.Gene)
 test$Top.D.Gene = gsub("[*]([0-9]+)", "", test$Top.D.Gene)
@@ -57,13 +57,18 @@
 #PRODF = PROD[ -1]
 
 PRODF = PROD
-
+print("test2\n")
 #PRODF = unique(PRODF)
-if(any(grepl(pattern="_", x=PRODF$ID))){
-
+if(any(grepl(pattern="_", x=PRODF$ID))){ #dumb and way to simple
 	PRODF$freq = gsub("^[0-9]+_", "", PRODF$ID)
 	PRODF$freq = gsub("_.*", "", PRODF$freq)
 	PRODF$freq = as.numeric(PRODF$freq)
+	if(any(is.na(PRODF$freq))){ #fix the dumbness if it fails
+		PRODF$freq = 1
+		if(selection == "unique"){
+			PRODF = PRODF[!duplicated(PRODF$VDJCDR3), ]
+		}
+	}
 } else {
 	PRODF$freq = 1
 	if(selection == "unique"){
@@ -96,6 +101,8 @@
 D = c("v.name\tchr.orderD\n")	
 J = c("v.name\tchr.orderJ\n")
 
+print("test3\n")
+
 if(species == "human"){
 	if(locus == "trb"){		
 		V = c("v.name\tchr.orderV\nTRBV2\t1\nTRBV3-1\t2\nTRBV4-1\t3\nTRBV5-1\t4\nTRBV6-1\t5\nTRBV4-2\t6\nTRBV6-2\t7\nTRBV4-3\t8\nTRBV6-3\t9\nTRBV7-2\t10\nTRBV6-4\t11\nTRBV7-3\t12\nTRBV9\t13\nTRBV10-1\t14\nTRBV11-1\t15\nTRBV10-2\t16\nTRBV11-2\t17\nTRBV6-5\t18\nTRBV7-4\t19\nTRBV5-4\t20\nTRBV6-6\t21\nTRBV5-5\t22\nTRBV7-6\t23\nTRBV5-6\t24\nTRBV6-8\t25\nTRBV7-7\t26\nTRBV6-9\t27\nTRBV7-8\t28\nTRBV5-8\t29\nTRBV7-9\t30\nTRBV13\t31\nTRBV10-3\t32\nTRBV11-3\t33\nTRBV12-3\t34\nTRBV12-4\t35\nTRBV12-5\t36\nTRBV14\t37\nTRBV15\t38\nTRBV16\t39\nTRBV18\t40\nTRBV19\t41\nTRBV20-1\t42\nTRBV24-1\t43\nTRBV25-1\t44\nTRBV27\t45\nTRBV28\t46\nTRBV29-1\t47\nTRBV30\t48")
@@ -131,6 +138,7 @@
 	cat("No D Genes in this species/locus")
 }
 
+print("test4\n")
 
 tcV = textConnection(V)
 Vchain = read.table(tcV, sep="\t", header=TRUE)
@@ -181,6 +189,8 @@
 pJ
 dev.off();
 
+print("test5\n")
+
 VGenes = PRODF[,c("Sample", "Top.V.Gene", "freq")]
 VGenes$Top.V.Gene = gsub("-.*", "", VGenes$Top.V.Gene)
 VGenes = data.frame(data.table(VGenes)[, list(Count=sum(freq)), by=c("Sample", "Top.V.Gene")])
@@ -245,6 +255,7 @@
 revVchain$chr.orderV = rev(revVchain$chr.orderV)
 revDchain$chr.orderD = rev(revDchain$chr.orderD)
 
+print("test6\n")
 
 plotVD <- function(dat){
 	if(length(dat[,1]) == 0){
@@ -353,6 +364,7 @@
 writeLines(un, sampleFile)
 close(sampleFile)
 
+print("test7\n")
 
 if("Replicate" %in% colnames(test))
 {
@@ -439,6 +451,8 @@
 	lapply(clonalityOverviewSplit, FUN=ClonalityOverviewPrint)
 }
 
+print("test8\n")
+
 if("Functionality" %in% colnames(test))
 {
 	newData = data.frame(data.table(PROD)[,list(unique=.N, 
@@ -467,3 +481,5 @@
 				by=c("Sample")])
 	write.table(newData, "junctionAnalysis.csv" , sep=",",quote=F,na="-",row.names=F,col.names=F)
 }
+
+print("test9\n")