diff sequence_overview.r @ 99:86206431cbb0 draft

Uploaded
author davidvanzessen
date Thu, 16 Jun 2016 10:01:54 -0400
parents 5ffbf40cdd4b
children ff5be711382b
line wrap: on
line diff
--- a/sequence_overview.r	Thu Jun 16 05:05:47 2016 -0400
+++ b/sequence_overview.r	Thu Jun 16 10:01:54 2016 -0400
@@ -3,10 +3,9 @@
 args <- commandArgs(trailingOnly = TRUE)
 
 input.file = args[1]
-input.file = args[2]
-outputdir = args[3]
-gene.classes = unlist(strsplit(args[4], ","))
-hotspot.analysis.sum.file = args[5]
+outputdir = args[2]
+gene.classes = unlist(strsplit(args[3], ","))
+hotspot.analysis.sum.file = args[4]
 NToverview.file = paste(outputdir, "ntoverview.txt", sep="/")
 NTsum.file = paste(outputdir, "ntsum.txt", sep="/")
 main.html = "index.html"
@@ -41,12 +40,17 @@
 make.link = function(id, clss, val) { paste("<a href='", clss, "_", id, ".html'>", val, "</a>", sep="") }
 tbl = function(df) { res = "<table border='1'>"; for(i in 1:nrow(df)){ res = paste(res, tr(df[i,]), sep=""); }; res = paste(res, "</table>"); }
 
-print(paste("Number of unique sequences to be written to the sequence overview page", nrow(dat)))
-
 cat("<table border='1'>", file=main.html, append=F)
 cat("<caption>CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once</caption>", file=main.html, append=T)
 cat("<tr><th>Sequence</th><th>Functionality</th><th>ca1</th><th>ca2</th><th>cg1</th><th>cg2</th><th>cg3</th><th>cg4</th><th>cm</th></tr>", file=main.html, append=T)
 
+
+
+single.sequences=0 #sequence only found once, skipped
+in.multiple=0 #same sequence across multiple subclasses
+multiple.in.one=0 #same sequence multiple times in one subclass
+matched=0 #should be the same als matched sequences
+
 for(i in 1:nrow(dat)){
 	
 	ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^ca1", IDs$best_match),]
@@ -65,10 +69,16 @@
 	classes.sum = sum(classes)
 	
 	if(classes.sum == 1){
-		print(paste("next", i, classes.sum))
+		single.sequences = single.sequences + 1
 		next
+	}
+	
+	matched = matched + sum(classes > 0) #count in how many subclasses the sequence occurs.
+	
+	if(any(classes  == classes.sum)){
+		in.multiple = in.multiple + 1
 	} else {
-		print(i)
+		multiple.in.one = multiple.in.one + 1
 	}
 	
 	id = as.numeric(dat[i,"seq_conc"])
@@ -120,6 +130,10 @@
 
 cat("</table>", file=main.html, append=T)
 
+print(paste("Single sequences:", single.sequences))
+print(paste("Sequences in multiple subclasses:", in.multiple))
+print(paste("Multiple sequences in one subclass:", multiple.in.one))
+print(paste("Count that should match 'matched' sequences:", matched))
 
 #ACGT overview
 
@@ -131,8 +145,6 @@
 NToverview$G = nchar(gsub("[^Gg]", "", NToverview$seq))
 NToverview$T = nchar(gsub("[^Tt]", "", NToverview$seq))
 
-print(sum(colSums(NToverview[,c("A", "C", "T", "G")])))
-
 #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)
@@ -166,8 +178,6 @@
 
 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")