Mercurial > repos > davidvanzessen > mutation_analysis

changeset 104:603a10976e9c draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Wed, 22 Jun 2016 10:07:28 -0400
parents e21cbe15381f
children e4957ad476a2
files aa_histogram.r mutation_analysis.py sequence_overview.r wrapper.sh
diffstat 4 files changed, 133 insertions(+), 30 deletions(-) [+]
line wrap: on
line diff
--- a/aa_histogram.r	Tue Jun 21 05:30:16 2016 -0400
+++ b/aa_histogram.r	Wed Jun 22 10:07:28 2016 -0400
@@ -4,6 +4,7 @@
 
 input = args[1]
 outfile = args[2]
+gene = args[3]
 
 print("---------------- read input ----------------")
 
@@ -40,12 +41,12 @@
 m = m + annotate("segment", x = 38.5, y = -0.05, xend=55.5, yend=-0.05, colour="darkgreen", size=1) + annotate("text", x = 47, y = -0.1, label="FR2")
 m = m + annotate("segment", x = 55.5, y = -0.07, xend=65.5, yend=-0.07, colour="darkblue", size=1) + annotate("text", x = 60.5, y = -0.15, label="CDR2")
 m = m + annotate("segment", x = 65.5, y = -0.05, xend=104.5, yend=-0.05, colour="darkgreen", size=1) + annotate("text", x = 85, y = -0.1, label="FR3")
-m = m + expand_limits(y=c(-0.1,1)) + xlab("AA position") + ylab("Frequency") + ggtitle("AA mutation frequency")
+m = m + expand_limits(y=c(-0.1,1)) + xlab("AA position") + ylab("Frequency") + ggtitle(paste(gene, "AA mutation frequency"))
 
 print("---------------- write/print ----------------")
 
 print("writing dat_dt") #need this
-write.table(dat_dt, paste(dirname(outfile), "/aa_histogram.txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
+write.table(dat_dt, paste(dirname(outfile), "/aa_histogram_", gene, ".txt", sep=""), sep="\t",quote=F,row.names=F,col.names=T)
 print("writing png") #also need this, file is haunted
 png(filename=outfile, width=1280, height=720)
 print(m)
--- a/mutation_analysis.py	Tue Jun 21 05:30:16 2016 -0400
+++ b/mutation_analysis.py	Wed Jun 22 10:07:28 2016 -0400
@@ -81,6 +81,7 @@
 AALength = (int(max(mutationList, key=lambda i: int(i[4]) if i[4] else 0)[4]) + 1)  # [4] is the position of the AA mutation, None if silent
 
 AA_mutation = [0] * AALength
+AA_mutation_dic = {"ca": AA_mutation[:], "cg": AA_mutation[:], "cm": AA_mutation[:], "un": AA_mutation[:]}
 AA_mutation_empty = AA_mutation[:]
 
 aa_mutations_by_id_file = outfile[:outfile.rindex("/")] + "/aa_id_mutations.txt"
@@ -89,8 +90,11 @@
 		AA_mutation_for_ID = AA_mutation_empty[:]
 		for mutation in mutationListByID[ID]:
 			if mutation[4]:
-				AA_mutation[int(mutation[4])] += 1
-				AA_mutation_for_ID[int(mutation[4])] += 1
+				AA_mutation_position = int(mutation[4])
+				AA_mutation[AA_mutation_position] += 1
+				AA_mutation_for_ID[AA_mutation_position] += 1
+				clss = genedic[ID][:2]
+				AA_mutation_dic[clss][AA_mutation_position] += 1
 		o.write(ID + "\t" + "\t".join([str(x) for x in AA_mutation_for_ID[1:]]) + "\n")
 
 
@@ -149,6 +153,27 @@
 	o.write("mutations.at.position\t" + "\t".join([str(x) for x in AA_mutation[1:]]) + "\n")
 	o.write("AA.at.position\t" + "\t".join([str(x) for x in absentAA]) + "\n")
 
+
+aa_mutations_file_ca = outfile[:outfile.rindex("/")] + "/aa_mutations_ca.txt"
+with open(aa_mutations_file_ca, 'w') as o:
+	o.write("row.name\t" + "\t".join([str(x) for x in range(1, AALength-1)]) + "\n")
+	o.write("mutations.at.position\t" + "\t".join([str(x) for x in AA_mutation_dic["ca"][1:]]) + "\n")
+	o.write("AA.at.position\t" + "\t".join([str(x) for x in absentAA]) + "\n")
+
+
+aa_mutations_file_cg = outfile[:outfile.rindex("/")] + "/aa_mutations_cg.txt"
+with open(aa_mutations_file_cg, 'w') as o:
+	o.write("row.name\t" + "\t".join([str(x) for x in range(1, AALength-1)]) + "\n")
+	o.write("mutations.at.position\t" + "\t".join([str(x) for x in AA_mutation_dic["cg"][1:]]) + "\n")
+	o.write("AA.at.position\t" + "\t".join([str(x) for x in absentAA]) + "\n")
+
+
+aa_mutations_file_cm = outfile[:outfile.rindex("/")] + "/aa_mutations_cm.txt"
+with open(aa_mutations_file_cm, 'w') as o:
+	o.write("row.name\t" + "\t".join([str(x) for x in range(1, AALength-1)]) + "\n")
+	o.write("mutations.at.position\t" + "\t".join([str(x) for x in AA_mutation_dic["cg"][1:]]) + "\n")
+	o.write("AA.at.position\t" + "\t".join([str(x) for x in absentAA]) + "\n")
+
 if linecount == 0:
 	print "No data, exiting"
 	with open(outfile, 'w') as o:
--- a/sequence_overview.r	Tue Jun 21 05:30:16 2016 -0400
+++ b/sequence_overview.r	Wed Jun 22 10:07:28 2016 -0400
@@ -39,14 +39,41 @@
 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>"); }
+get.bg.color = function(val){
+	if(val %in% c("TRUE", "FALSE", "T", "F")){ #if its a logical value, give the background a green/red color
+		return(ifelse(val,"#eafaf1","#f9ebea"))
+	} else if (!is.na(as.numeric(val))) { #if its a numerical value, give it a grey tint if its >0
+		return(ifelse(val > 0,"#eaecee","white"))
+	} else {
+		return("white")
+	}
+}
+td = function(val) {
+  return(paste("<td bgcolor='", get.bg.color(val), "'>", val, "</td>", sep=""))
+}
+tr = function(val) { 
+	return(paste(c("<tr>", sapply(val, td), "</tr>"), collapse="")) 
+}
+
+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><th>un</th></tr>", file=main.html, append=T)
+cat("<tr>", file=main.html, append=T)
+cat("<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><th>un</th>", file=main.html, append=T)
+cat("<th>total CA</th><th>total CG</th><th>number of subclasses</th><th>present in both Ca and Cg</th><th>Ca1+Ca2</th>", file=main.html, append=T)
+cat("<th>Cg1+Cg2</th><th>Cg1+Cg3</th><th>Cg1+Cg4</th><th>Cg2+Cg3</th><th>Cg2+Cg4</th><th>Cg3+Cg4</th>", file=main.html, append=T)
+cat("<th>Cg1+Cg2+Cg3</th><th>Cg2+Cg3+Cg4</th><th>Cg1+Cg2+Cg4</th><th>Cg1+Cg3+Cg4</th>", file=main.html, append=T)
+cat("</tr>", file=main.html, append=T)
 
 
 
@@ -74,7 +101,19 @@
 	un = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^unmatched", IDs$best_match),]
 	allc = rbind(ca1, ca2, cg1, cg2, cg3, cg4, cm, un)
 	
-	classes = c(nrow(ca1), nrow(ca2), nrow(cg1), nrow(cg2), nrow(cg3), nrow(cg4), nrow(cm), nrow(un))
+	ca1.n = nrow(ca1)
+	ca2.n = nrow(ca2)
+	
+	cg1.n = nrow(cg1)
+	cg2.n = nrow(cg2)
+	cg3.n = nrow(cg3)
+	cg4.n = nrow(cg4)
+	
+	cm.n = nrow(cm)
+	
+	un.n = nrow(un)
+	
+	classes = c(ca1.n, ca2.n, cg1.n, cg2.n, cg3.n, cg4.n, cm.n, un.n)
 	
 	classes.sum = sum(classes)
 	
@@ -83,16 +122,18 @@
 		next
 	}
 	
-	if(nrow(un) == classes.sum){
+	if(un.n == classes.sum){
 		unmatched = unmatched + 1
 		next
 	}
 	
-	matched = matched + sum(classes > 0) #count in how many subclasses the sequence occurs.
+	in.classes = sum(classes > 0)
+	
+	matched = matched + in.classes #count in how many subclasses the sequence occurs.
 	
 	if(any(classes  == classes.sum)){
 		multiple.in.one = multiple.in.one + 1
-	} else if (nrow(un) > 0) {
+	} else if (un.n > 0) {
 		some.unmatched = some.unmatched + 1
 	} else {
 		in.multiple = in.multiple + 1
@@ -104,52 +145,79 @@
 	
 	by.id.row = c()
 	
-	if(nrow(ca1) > 0){
+	if(ca1.n > 0){
 		cat(tbl(ca1), file=paste("ca1_", id, ".html", sep=""))
 	}
 
-	if(nrow(ca2) > 0){
+	if(ca2.n > 0){
 		cat(tbl(ca2), file=paste("ca2_", id, ".html", sep=""))
 	}
 
-	if(nrow(cg1) > 0){
+	if(cg1.n > 0){
 		cat(tbl(cg1), file=paste("cg1_", id, ".html", sep=""))
 	}
 
-	if(nrow(cg2) > 0){
+	if(cg2.n > 0){
 		cat(tbl(cg2), file=paste("cg2_", id, ".html", sep=""))
 	}
 
-	if(nrow(cg3) > 0){
+	if(cg3.n > 0){
 		cat(tbl(cg3), file=paste("cg3_", id, ".html", sep=""))
 	}
 
-	if(nrow(cg4) > 0){
+	if(cg4.n > 0){
 		cat(tbl(cg4), file=paste("cg4_", id, ".html", sep=""))
 	}
 
-	if(nrow(cm) > 0){
+	if(cm.n > 0){
 		cat(tbl(cm), file=paste("cm_", id, ".html", sep=""))
 	}
 
-	if(nrow(un) > 0){
+	if(un.n > 0){
 		cat(tbl(un), file=paste("un_", id, ".html", sep=""))
 	}
 	
-	ca1.html = make.link(id, "ca1", nrow(ca1))
-	ca2.html = make.link(id, "ca2", nrow(ca2))
+	ca1.html = make.link(id, "ca1", ca1.n)
+	ca2.html = make.link(id, "ca2", ca2.n)
+	
+	cg1.html = make.link(id, "cg1", cg1.n)
+	cg2.html = make.link(id, "cg2", cg2.n)
+	cg3.html = make.link(id, "cg3", cg3.n)
+	cg4.html = make.link(id, "cg4", cg4.n)
+	
+	cm.html = make.link(id, "cm", cm.n)
+	
+	un.html = make.link(id, "un", un.n)
+	
+	#extra columns
+	ca.n = ca1.n + ca2.n
+	
+	cg.n = cg1.n + cg2.n + cg3.n + cg4.n
+	
+	#in.classes
 	
-	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))
+	in.ca.cg = (ca.n > 0 & cg.n > 0)
+	
+	in.ca1.ca2 = (ca1.n > 0 & ca2.n > 0)
 	
-	cm.html = make.link(id, "cm", nrow(cm))
+	in.cg1.cg2 = (cg1.n > 0 & cg2.n > 0)
+	in.cg1.cg3 = (cg1.n > 0 & cg3.n > 0)
+	in.cg1.cg4 = (cg1.n > 0 & cg4.n > 0)
+	in.cg2.cg3 = (cg2.n > 0 & cg3.n > 0)
+	in.cg2.cg4 = (cg2.n > 0 & cg4.n > 0)
+	in.cg3.cg4 = (cg3.n > 0 & cg4.n > 0)
 	
-	un.html = make.link(id, "un", nrow(un))
+	in.cg1.cg2.cg3 = (cg1.n > 0 & cg2.n > 0 & cg3.n > 0)
+	in.cg2.cg3.cg4 = (cg2.n > 0 & cg3.n > 0 & cg4.n > 0)
+	in.cg1.cg2.cg4 = (cg1.n > 0 & cg2.n > 0 & cg4.n > 0)
+	in.cg1.cg3.cg4 = (cg1.n > 0 & cg3.n > 0 & cg4.n > 0)
+	
+	
+	
 	
 	#rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html, un.html)
 	rw = c(as.character(dat[i,"seq_conc"]), functionality, ca1.html, ca2.html, cg1.html, cg2.html, cg3.html, cg4.html, cm.html, un.html)
+	rw = c(rw, ca.n, cg.n, in.classes, in.ca.cg, in.ca1.ca2, in.cg1.cg2, in.cg1.cg3, in.cg1.cg4, in.cg2.cg3, in.cg2.cg4, in.cg3.cg4, in.cg1.cg2.cg3, in.cg2.cg3.cg4, in.cg1.cg2.cg4, in.cg1.cg3.cg4)
 
 	cat(tr(rw), file=main.html, append=T)
 	
--- a/wrapper.sh	Tue Jun 21 05:30:16 2016 -0400
+++ b/wrapper.sh	Wed Jun 22 10:07:28 2016 -0400
@@ -148,7 +148,10 @@
 echo "---------------- aa_histogram.r ----------------"
 echo "---------------- aa_histogram.r ----------------<br />" >> $output
 
-Rscript $dir/aa_histogram.r $outdir/aa_mutations.txt $outdir/aa_histogram.png 2>&1
+Rscript $dir/aa_histogram.r $outdir/aa_mutations.txt $outdir/aa_histogram.png "" 2>&1
+Rscript $dir/aa_histogram.r $outdir/aa_mutations_ca.txt $outdir/aa_histogram_ca.png "ca" 2>&1
+Rscript $dir/aa_histogram.r $outdir/aa_mutations_cg.txt $outdir/aa_histogram_cg.png "cg" 2>&1
+Rscript $dir/aa_histogram.r $outdir/aa_mutations_cm.txt $outdir/aa_histogram_cm.png "cm" 2>&1
 
 genes=(ca ca1 ca2 cg cg1 cg2 cg3 cg4 cm)
 
@@ -279,6 +282,12 @@
 then
 	echo "<img src='aa_histogram.png'/><br />" >> $output
 	echo "<a href='aa_histogram.txt'>download data</a><br />" >> $output
+	echo "<img src='aa_histogram_ca.png'/><br />" >> $output
+	echo "<a href='aa_histogram_ca.txt'>download data</a><br />" >> $output
+	echo "<img src='aa_histogram_cg.png'/><br />" >> $output
+	echo "<a href='aa_histogram_cg.txt'>download data</a><br />" >> $output
+	echo "<img src='aa_histogram_cm.png'/><br />" >> $output
+	echo "<a href='aa_histogram_cm.txt'>download data</a><br />" >> $output
 fi
 
 for gene in ${genes[@]}