diff sequence_overview.r @ 82:564c4f6da203 draft

Uploaded
author davidvanzessen
date Tue, 24 May 2016 03:52:14 -0400
parents a778156dad3d
children 0011a0597685
line wrap: on
line diff
--- a/sequence_overview.r	Tue May 17 04:17:38 2016 -0400
+++ b/sequence_overview.r	Tue May 24 03:52:14 2016 -0400
@@ -4,7 +4,10 @@
 
 gene.matches = args[1]
 sequence.file = args[2]
-outputdir = args[3]
+merged.file = args[3]
+outputdir = args[4]
+gene.classes = unlist(strsplit(args[5], ","))
+hotspot.analysis.sum.file = args[6]
 NToverview.file = paste(outputdir, "ntoverview.txt", sep="/")
 NTsum.file = paste(outputdir, "ntsum.txt", sep="/")
 main.html = "index.html"
@@ -13,9 +16,13 @@
 
 genes = read.table(gene.matches, header=T, sep="\t", fill=T)
 sequences = read.table(sequence.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
+merged = read.table(merged.file, header=T, sep="\t", fill=T, stringsAsFactors=F, quote="")
+hotspot.analysis.sum = read.table(hotspot.analysis.sum.file, header=F, sep=",", fill=T, stringsAsFactors=F, quote="")
 
 dat = merge(sequences, genes, by="Sequence.ID")
 
+dat = dat[dat$Sequence.ID %in% merged$Sequence.ID,]
+
 dat$seq_conc = paste(dat$CDR1.IMGT, dat$CDR2.IMGT, dat$CDR3.IMGT, dat$FR2.IMGT, dat$FR3.IMGT)
 
 IDs = dat[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")]
@@ -23,15 +30,15 @@
 
 #dat = data.frame(data.table(dat)[, list(freq=.N), by=c("best_match", "seq_conc")])
 
-dat = data.frame(table(dat$seq_conc))
+dat = data.frame(table(dat$seq_conc, dat$Functionality))
 
 dat = dat[dat$Freq > 1,]
 
-names(dat) = c("seq_conc", "Freq")
+names(dat) = c("seq_conc", "Functionality", "Freq")
 
 dat$seq_conc = factor(dat$seq_conc)
 
-dat = dat[order(nchar(as.character(dat$seq_conc))),]
+dat = dat[order(as.character(dat$seq_conc)),]
 
 #writing html from R...
 td = function(val) { paste("<td>", val, "</td>", sep="") }
@@ -54,6 +61,14 @@
 	
 	cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cm",]
 	
+	classes = c(nrow(ca1), nrow(ca2), nrow(cg1), nrow(cg2), nrow(cg3), nrow(cg4), nrow(cm))
+	
+	classes.sum = sum(classes)
+	
+	if(!any(classes != 0 & classes != classes.sum)){
+		next
+	}
+	
 	id = as.numeric(dat[i,"seq_conc"])
 	
 	functionality = dat[i,"Functionality"]
@@ -113,20 +128,48 @@
 
 NToverview = merge(NToverview, sequences[,c("Sequence.ID", "seq")], by="Sequence.ID")
 
+NToverview = NToverview[NToverview$Sequence.ID %in% merged$Sequence.ID,]
+
 NToverview$A = nchar(gsub("[^Aa]", "", NToverview$seq))
 NToverview$C = nchar(gsub("[^Cc]", "", NToverview$seq))
 NToverview$G = nchar(gsub("[^Gg]", "", NToverview$seq))
 NToverview$T = nchar(gsub("[^Tt]", "", NToverview$seq))
 
-NTsum = data.frame(Sequence.ID="-", best_match="Sum", seq="-", A = sum(NToverview$A), C = sum(NToverview$C), G = sum(NToverview$G), T = sum(NToverview$T))
+#Nsum = data.frame(Sequence.ID="-", best_match="Sum", seq="-", A = sum(NToverview$A), C = sum(NToverview$C), G = sum(NToverview$G), T = sum(NToverview$T))
+
+#NToverview = rbind(NToverview, NTsum)
+
+NTresult = data.frame(nt=c("A", "C", "T", "G"))
+
+for(clazz in gene.classes){
+	NToverview.sub = NToverview[grepl(clazz, NToverview$best_match),]
+	new.col.x = c(sum(NToverview.sub$A), sum(NToverview.sub$C), sum(NToverview.sub$T), sum(NToverview.sub$G))
+	new.col.y = sum(new.col.x)
+	new.col.z = round(new.col.x / new.col.y * 100, 2)
+	
+	tmp = names(NTresult)
+	NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z))
+	names(NTresult) = c(tmp, paste(clazz, c("x", "y", "z"), sep=""))
+}
 
-print(names(NToverview))
-print(names(NTsum))
+new.col.x = c(sum(NToverview$A), sum(NToverview$C), sum(NToverview$T), sum(NToverview$G))
+new.col.y = sum(new.col.x)
+new.col.z = round(new.col.x / new.col.y * 100, 2)
+
+tmp = names(NTresult)
+NTresult = cbind(NTresult, data.frame(new.col.x, new.col.y, new.col.z))
+names(NTresult) = c(tmp, paste("all", c("x", "y", "z"), sep=""))
 
-NToverview = rbind(NToverview, NTsum)
+names(hotspot.analysis.sum) = names(NTresult)
+
+hotspot.analysis.sum = rbind(hotspot.analysis.sum, NTresult)
+
+print(hotspot.analysis.sum)
+
+write.table(hotspot.analysis.sum, hotspot.analysis.sum.file, quote=F, sep=",", row.names=F, col.names=F, na="0")
+
 
 write.table(NToverview, NToverview.file, quote=F, sep="\t", row.names=F, col.names=T)
-#write.table(NTsum, NTsum.file, quote=F, sep="\t", row.names=F, col.names=T)