Mercurial > repos > davidvanzessen > mutation_analysis

changeset 92:b869a126e2c4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Mon, 06 Jun 2016 08:26:54 -0400
parents 5e237c243088
children 53fb2948726e
files gene_identification.py sequence_overview.r sequence_overview.r.old wrapper.sh
diffstat 4 files changed, 5 insertions(+), 212 deletions(-) [+]
line wrap: on
line diff
--- a/gene_identification.py	Fri Jun 03 10:35:52 2016 -0400
+++ b/gene_identification.py	Mon Jun 06 08:26:54 2016 -0400
@@ -160,9 +160,9 @@
 chunksInCM = len(compiledregex["cm"])
 requiredChunkPercentage = 0.7
 varsInCA = float(len(ca1.keys()) * 2)
-varsInCG = float(len(cg1.keys()) * 2) - 1 # -1 because the sliding window doesn't hit the first nt twice
+varsInCG = float(len(cg1.keys()) * 2) - 2 # -1 because the sliding window doesn't hit the first nt twice
 varsInCM = 0
-requiredVarPercentage = 0.7
+
 
 
 first = True
--- a/sequence_overview.r	Fri Jun 03 10:35:52 2016 -0400
+++ b/sequence_overview.r	Mon Jun 06 08:26:54 2016 -0400
@@ -66,7 +66,7 @@
 	
 	id = as.numeric(dat[i,"seq_conc"])
 	
-	functionality = paste(unique(allc[,"Functionality"], sep=","))
+	functionality = paste(unique(allc[,"Functionality"]), collapse=",")
 
 	if(nrow(ca1) > 0){
 		cat(tbl(ca1), file=paste("ca1_", id, ".html", sep=""))
@@ -106,8 +106,7 @@
 	
 	cm.html = make.link(id, "cm", nrow(cm))
 	
-	rw = c(as.character(dat[i,"seq_conc"]), as.character(functionality), ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html)
-	#print(rw)
+	rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html)
 
 	cat(tr(rw), file=main.html, append=T)
 }
@@ -117,8 +116,6 @@
 
 #ACGT overview
 
-
-
 NToverview = merged
 NToverview$seq = paste(NToverview$CDR1.IMGT.seq, NToverview$FR2.IMGT.seq, NToverview$CDR2.IMGT.seq, NToverview$FR3.IMGT.seq, sep="_")
 
--- a/sequence_overview.r.old	Fri Jun 03 10:35:52 2016 -0400
+++ /dev/null	Thu Jan 01 00:00:00 1970 +0000
@@ -1,203 +0,0 @@
-library(reshape2)
-
-args <- commandArgs(trailingOnly = TRUE)
-
-gene.matches = args[1]
-sequence.file = args[2]
-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"
-
-setwd(outputdir)
-
-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$FR2.IMGT, dat$CDR2.IMGT, dat$FR3.IMGT, dat$CDR3.IMGT)
-#dat$seq_conc = paste(dat$CDR1.IMGT, dat$CDR2.IMGT, dat$CDR3.IMGT)
-
-IDs = dat[,c("Sequence.ID", "seq_conc", "best_match", "Functionality")]
-IDs$best_match = as.character(IDs$best_match)
-
-#dat = data.frame(data.table(dat)[, list(freq=.N), by=c("best_match", "seq_conc")])
-
-dat = data.frame(table(dat$seq_conc, dat$Functionality))
-
-#dat = dat[dat$Freq > 1,]
-
-names(dat) = c("seq_conc", "Functionality", "Freq")
-
-dat$seq_conc = factor(dat$seq_conc)
-
-dat = dat[order(as.character(dat$seq_conc)),]
-
-#writing html from R...
-td = function(val) { paste("<td>", val, "</td>", sep="") }
-tr = function(val) { capture.output(cat("<tr>", td(val), "</tr>", sep="")) }
-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>"); }
-
-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)
-
-for(i in 1:nrow(dat)){
-	ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "ca1",]
-	ca2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "ca2",]
-	
-	cg1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg1",]
-	cg2 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg2",]
-	cg3 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg3",]
-	cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & IDs$best_match == "cg4",]
-	
-	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(classes.sum == 1){
-		next
-	}
-	
-	id = as.numeric(dat[i,"seq_conc"])
-	
-	functionality = dat[i,"Functionality"]
-
-	if(nrow(ca1) > 0){
-		cat(tbl(ca1), file=paste("ca1_", id, ".html", sep=""))
-	}
-
-	if(nrow(ca2) > 0){
-		cat(tbl(ca2), file=paste("ca2_", id, ".html", sep=""))
-	}
-
-	if(nrow(cg1) > 0){
-		cat(tbl(cg1), file=paste("cg1_", id, ".html", sep=""))
-	}
-
-	if(nrow(cg2) > 0){
-		cat(tbl(cg2), file=paste("cg2_", id, ".html", sep=""))
-	}
-
-	if(nrow(cg3) > 0){
-		cat(tbl(cg3), file=paste("cg3_", id, ".html", sep=""))
-	}
-
-	if(nrow(cg4) > 0){
-		cat(tbl(cg4), file=paste("cg4_", id, ".html", sep=""))
-	}
-
-	if(nrow(cm) > 0){
-		cat(tbl(cm), file=paste("cm_", id, ".html", sep=""))
-	}
-	
-	ca1.html = make.link(id, "ca1", nrow(ca1))
-	ca2.html = make.link(id, "ca2", nrow(ca2))
-	
-	cg1.html = make.link(id, "cg1", nrow(cg1))
-	cg2.html = make.link(id, "cg2", nrow(cg2))
-	cg3.html = make.link(id, "cg3", nrow(cg3))
-	cg4.html = make.link(id, "cg4", nrow(cg4))
-	
-	cm.html = make.link(id, "cm", nrow(cm))
-	
-	rw = c(as.character(dat[i,"seq_conc"]), as.character(functionality), ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html)
-
-	cat(tr(rw), file=main.html, append=T)
-}
-
-cat("</table>", file=main.html, append=T)
-
-
-#ACGT overview
-
-
-
-NToverview = genes[,c("Sequence.ID", "best_match")]
-sequences$seq = paste(sequences$CDR1.IMGT, sequences$FR2.IMGT, sequences$CDR2.IMGT, sequences$FR3.IMGT, sep="_")
-
-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))
-
-#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, paste("^", NToverview$best_match, sep="")),]
-	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=""))
-}
-
-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=""))
-
-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)
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
-
--- a/wrapper.sh	Fri Jun 03 10:35:52 2016 -0400
+++ b/wrapper.sh	Mon Jun 06 08:26:54 2016 -0400
@@ -55,7 +55,6 @@
 
 echo "${BLASTN_DIR}"
 
-echo "identification ($method)"
 echo "---------------- identification ($method) ----------------"
 echo "---------------- identification ($method) ----------------<br />" >> $output
 
@@ -109,7 +108,7 @@
 
 echo "<table border='1'>" > $outdir/base_overview.html
 
-while read ID class seq A C G T
+while IFS=$'\t' read ID class seq A C G T
 do
 	echo "<tr><td>$ID</td><td>$seq</td><td>$class</td><td>$A</td><td>$C</td><td>$G</td><td>$T</td></tr>" >> $outdir/base_overview.html
 done < $outdir/sequence_overview/ntoverview.txt