Mercurial > repos > davidvanzessen > mutation_analysis

changeset 102:e6bc976760d4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author davidvanzessen
date Tue, 21 Jun 2016 03:32:50 -0400
parents 3cffb8a38bb1
children e21cbe15381f
files aa_histogram.r sequence_overview.r tmp/igat.r wrapper.sh
diffstat 4 files changed, 85 insertions(+), 30 deletions(-) [+]
line wrap: on
line diff
--- a/aa_histogram.r	Fri Jun 17 08:31:20 2016 -0400
+++ b/aa_histogram.r	Tue Jun 21 03:32:50 2016 -0400
@@ -44,8 +44,9 @@
 
 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)
-
+print("writing png") #also need this, file is haunted
 png(filename=outfile, width=1280, height=720)
 print(m)
 dev.off()
--- a/sequence_overview.r	Fri Jun 17 08:31:20 2016 -0400
+++ b/sequence_overview.r	Tue Jun 21 03:32:50 2016 -0400
@@ -17,7 +17,7 @@
 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="")
 
-before.unique = before.unique[!grepl("unmatched", before.unique$best_match),]
+#before.unique = before.unique[!grepl("unmatched", before.unique$best_match),]
 
 before.unique$seq_conc = paste(before.unique$CDR1.IMGT.seq, before.unique$FR2.IMGT.seq, before.unique$CDR2.IMGT.seq, before.unique$FR3.IMGT.seq, before.unique$CDR3.IMGT.seq)
 
@@ -46,15 +46,19 @@
 
 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)
+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)
 
 
 
 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
+unmatched=0 #all of the sequences are unmatched
+some.unmatched=0 #one or more sequences in a clone are unmatched
 matched=0 #should be the same als matched sequences
 
+sequence.id.page="by_id.html"
+
 for(i in 1:nrow(dat)){
 	
 	ca1 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^ca1", IDs$best_match),]
@@ -66,9 +70,11 @@
 	cg4 = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cg4", IDs$best_match),]
 	
 	cm = IDs[IDs$seq_conc == dat[i,c("seq_conc")] & grepl("^cm", IDs$best_match),]
-	allc = rbind(ca1, ca2, cg1, cg2, cg3, cg4, cm)
 	
-	classes = c(nrow(ca1), nrow(ca2), nrow(cg1), nrow(cg2), nrow(cg3), nrow(cg4), nrow(cm))
+	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))
 	
 	classes.sum = sum(classes)
 	
@@ -77,18 +83,27 @@
 		next
 	}
 	
+	if(nrow(un) == classes.sum){
+		unmatched = unmatched + 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
+		multiple.in.one = multiple.in.one + 1
+	} else if (nrow(un) > 0) {
+		some.unmatched = some.unmatched + 1
 	} else {
-		multiple.in.one = multiple.in.one + 1
+		in.multiple = in.multiple + 1
 	}
 	
 	id = as.numeric(dat[i,"seq_conc"])
 	
 	functionality = paste(unique(allc[,"Functionality"]), collapse=",")
-
+	
+	by.id.row = c()
+	
 	if(nrow(ca1) > 0){
 		cat(tbl(ca1), file=paste("ca1_", id, ".html", sep=""))
 	}
@@ -116,6 +131,10 @@
 	if(nrow(cm) > 0){
 		cat(tbl(cm), file=paste("cm_", id, ".html", sep=""))
 	}
+
+	if(nrow(un) > 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))
@@ -127,9 +146,18 @@
 	
 	cm.html = make.link(id, "cm", nrow(cm))
 	
-	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 = make.link(id, "un", nrow(un))
+	
+	#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)
 
 	cat(tr(rw), file=main.html, append=T)
+	
+	
+	for(i in 1:nrow(allc)){ #generate html by id
+		html = make.link(id, allc[i,"best_match"], allc[i,"Sequence.ID"])
+		cat(paste(html, "<br />"), file=sequence.id.page, append=T)
+	}
 }
 
 cat("</table>", file=main.html, append=T)
@@ -137,6 +165,7 @@
 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("Matched with unmatched:", some.unmatched))
 print(paste("Count that should match 'matched' sequences:", matched))
 
 #ACGT overview
--- a/tmp/igat.r	Fri Jun 17 08:31:20 2016 -0400
+++ b/tmp/igat.r	Tue Jun 21 03:32:50 2016 -0400
@@ -6,8 +6,6 @@
 
 merged = read.table(merged.file, header=T, sep="\t", fill=T, stringsAsFactors=F)
 
-print(head(merged$best_match))
-
 if(gene != "-"){
 	merged = merged[grepl(gene, merged$best_match),]
 }
--- a/wrapper.sh	Fri Jun 17 08:31:20 2016 -0400
+++ b/wrapper.sh	Tue Jun 21 03:32:50 2016 -0400
@@ -3,8 +3,9 @@
 dir="$(cd "$(dirname "$0")" && pwd)"
 input=$1
 method=$2
-output=$3
+log=$3 #becomes the main html page at the end
 outdir=$4
+output="$outdir/index.html" #copied to $log location at the end
 title=$5
 include_fr1=$6
 functionality=$7
@@ -17,7 +18,7 @@
 mkdir $outdir
 
 echo "---------------- read parameters ----------------"
-echo "---------------- read parameters ----------------<br />" > $output
+echo "---------------- read parameters ----------------<br />" > $log
 
 echo "unpacking IMGT file"
 
@@ -56,7 +57,7 @@
 echo "${BLASTN_DIR}"
 
 echo "---------------- identification ($method) ----------------"
-echo "---------------- identification ($method) ----------------<br />" >> $output
+echo "---------------- identification ($method) ----------------<br />" >> $log
 
 if [[ "${method}" == "custom" ]] ; then
 	python $dir/gene_identification.py --input $PWD/summary.txt --output $outdir/identified_genes.txt
@@ -73,12 +74,12 @@
 fi
 
 echo "---------------- merge_and_filter.r ----------------"
-echo "---------------- merge_and_filter.r ----------------<br />" >> $output
+echo "---------------- merge_and_filter.r ----------------<br />" >> $log
 
 Rscript $dir/merge_and_filter.r $PWD/summary.txt $PWD/sequences.txt $PWD/mutationanalysis.txt $PWD/mutationstats.txt $PWD/hotspots.txt $outdir/identified_genes.txt $outdir/merged.txt $outdir/before_unique_filter.txt $outdir/unmatched.txt $method $functionality $unique ${filter_unique} ${class_filter} 2>&1
 
 echo "---------------- creating new IMGT zip ----------------"
-echo "---------------- creating new IMGT zip ----------------<br />" >> $output
+echo "---------------- creating new IMGT zip ----------------<br />" >> $log
 
 mkdir $outdir/new_IMGT
 
@@ -132,7 +133,7 @@
 cd $tmp
 
 echo "---------------- mutation_analysis.r ----------------"
-echo "---------------- mutation_analysis.r ----------------<br />" >> $output
+echo "---------------- mutation_analysis.r ----------------<br />" >> $log
 
 classes="ca,ca1,ca2,cg,cg1,cg2,cg3,cg4,cm,unmatched"
 echo "R mutation analysis"
@@ -140,7 +141,7 @@
 
 
 echo "---------------- mutation_analysis.py ----------------"
-echo "---------------- mutation_analysis.py ----------------<br />" >> $output
+echo "---------------- mutation_analysis.py ----------------<br />" >> $log
 
 python $dir/mutation_analysis.py --input $outdir/merged.txt --genes $classes --includefr1 "${include_fr1}" --output $outdir/hotspot_analysis.txt
 
@@ -154,7 +155,7 @@
 funcs=(sum mean median)
 
 echo "---------------- sequence_overview.r ----------------"
-echo "---------------- sequence_overview.r ----------------" >> $output
+echo "---------------- sequence_overview.r ----------------<br />" >> $log
 
 mkdir $outdir/sequence_overview
 
@@ -168,7 +169,6 @@
 	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
 
-
 echo "<html><center><h1>$title</h1></center>" > $output
 
 #display the matched/unmatched for clearity
@@ -185,10 +185,12 @@
 echo "<center><h2>Percentage unmatched: ${perc_count}</h2></center>" >> $output
 
 echo "---------------- main tables ----------------"
+echo "---------------- main tables ----------------<br />" >> $log
 for func in ${funcs[@]}
 do
 
 	echo "---------------- $func table ----------------"
+	echo "---------------- $func table ----------------<br />" >> $log
 	
 	cat $outdir/mutations_${func}.txt $outdir/hotspot_analysis_${func}.txt > $outdir/data_${func}.txt
 
@@ -218,6 +220,7 @@
 done
 
 echo "---------------- download links ----------------"
+echo "---------------- download links ----------------<br />" >> $log
 
 
 echo "<a href='unmatched.txt'>unmatched</a><br />" >> $output
@@ -247,6 +250,7 @@
 
 
 echo "---------------- images ----------------"
+echo "---------------- images ----------------<br />" >> $log
 
 echo "<img src='all.png'/><br />" >> $output
 echo "<a href='all.txt'>download data</a><br />" >> $output
@@ -297,44 +301,67 @@
 echo "</html>" >> $output
 
 echo "---------------- baseline ----------------"
+echo "---------------- baseline ----------------<br />" >> $log
 tmp="$PWD"
 
 mkdir $outdir/baseline
 
 
 mkdir $outdir/baseline/ca_cg_cm
-cd $outdir/baseline/ca_cg_cm
-bash $dir/tmp/baseline/wrapper.sh 1 1 1 1 0 0 "25:26:38:55:65:104:-" $outdir/new_IMGT.txz "ca_cg_cm" "$dir/tmp/baseline/IMGT-reference-seqs-IGHV-2015-11-05.fa" "$outdir/baseline.pdf" "Sequence.ID" "$outdir/baseline.txt"
+if [[ $(wc -l < $outdir/new_IMGT/1_Summary.txt) -gt "1" ]]; then
+	cd $outdir/baseline/ca_cg_cm
+	bash $dir/tmp/baseline/wrapper.sh 1 1 1 1 0 0 "25:26:38:55:65:104:-" $outdir/new_IMGT.txz "ca_cg_cm" "$dir/tmp/baseline/IMGT-reference-seqs-IGHV-2015-11-05.fa" "$outdir/baseline.pdf" "Sequence.ID" "$outdir/baseline.txt"	
+else
+	echo "No sequences" > "$outdir/baseline.txt"	
+fi
 
 mkdir $outdir/baseline/ca
-cd $outdir/baseline/ca
-bash $dir/tmp/baseline/wrapper.sh 1 1 1 1 0 0 "25:26:38:55:65:104:-" $outdir/new_IMGT_ca.txz "ca" "$dir/tmp/baseline/IMGT-reference-seqs-IGHV-2015-11-05.fa" "$outdir/baseline_ca.pdf" "Sequence.ID" "$outdir/baseline_ca.txt"
+if [[ $(wc -l < $outdir/new_IMGT_ca/1_Summary.txt) -gt "1" ]]; then
+	cd $outdir/baseline/ca
+	bash $dir/tmp/baseline/wrapper.sh 1 1 1 1 0 0 "25:26:38:55:65:104:-" $outdir/new_IMGT_ca.txz "ca" "$dir/tmp/baseline/IMGT-reference-seqs-IGHV-2015-11-05.fa" "$outdir/baseline_ca.pdf" "Sequence.ID" "$outdir/baseline_ca.txt"
+else
+	echo "No ca sequences" > "$outdir/baseline_ca.txt"	
+fi
 
 mkdir $outdir/baseline/cg
-cd $outdir/baseline/cg
-bash $dir/tmp/baseline/wrapper.sh 1 1 1 1 0 0 "25:26:38:55:65:104:-" $outdir/new_IMGT_cg.txz "cg" "$dir/tmp/baseline/IMGT-reference-seqs-IGHV-2015-11-05.fa" "$outdir/baseline_cg.pdf" "Sequence.ID" "$outdir/baseline_cg.txt"
+if [[ $(wc -l < $outdir/new_IMGT_cg/1_Summary.txt) -gt "1" ]]; then
+	cd $outdir/baseline/cg
+	bash $dir/tmp/baseline/wrapper.sh 1 1 1 1 0 0 "25:26:38:55:65:104:-" $outdir/new_IMGT_cg.txz "cg" "$dir/tmp/baseline/IMGT-reference-seqs-IGHV-2015-11-05.fa" "$outdir/baseline_cg.pdf" "Sequence.ID" "$outdir/baseline_cg.txt"
+else
+	echo "No cg sequences" > "$outdir/baseline_cg.txt"	
+fi
 
 mkdir $outdir/baseline/cm
-cd $outdir/baseline/cm
-bash $dir/tmp/baseline/wrapper.sh 1 1 1 1 0 0 "25:26:38:55:65:104:-" $outdir/new_IMGT_cm.txz "cm" "$dir/tmp/baseline/IMGT-reference-seqs-IGHV-2015-11-05.fa" "$outdir/baseline_cm.pdf" "Sequence.ID" "$outdir/baseline_cm.txt"
+if [[ $(wc -l < $outdir/new_IMGT_cm/1_Summary.txt) -gt "1" ]]; then
+	cd $outdir/baseline/cm
+	bash $dir/tmp/baseline/wrapper.sh 1 1 1 1 0 0 "25:26:38:55:65:104:-" $outdir/new_IMGT_cm.txz "cm" "$dir/tmp/baseline/IMGT-reference-seqs-IGHV-2015-11-05.fa" "$outdir/baseline_cm.pdf" "Sequence.ID" "$outdir/baseline_cm.txt"
+else
+	echo "No cm sequences" > "$outdir/baseline_cm.txt"	
+fi
 
 cd $tmp
 
-#optional output for naive
-
 echo "---------------- naive_output.r ----------------"
+echo "---------------- naive_output.r ----------------<br />" >> $log
 
 if [[ "$naive_output" != "None" ]]
 then
 	echo "---------------- imgt_loader.r ----------------"
+	echo "---------------- imgt_loader.r ----------------<br />" >> $log
 	#python $dir/imgt_loader.py --summ $PWD/summary.txt --aa $PWD/aa.txt --junction $PWD/junction.txt --output $naive_output
 	Rscript --verbose $dir/imgt_loader.r $PWD/summary.txt $PWD/aa.txt $PWD/junction.txt $outdir/loader_output.txt 2>&1
 
 	echo "---------------- naive_output.r ----------------"
+	echo "---------------- naive_output.r ----------------<br />" >> $log
 	Rscript $dir/naive_output.r $outdir/loader_output.txt $outdir/merged.txt ${naive_output_ca} ${naive_output_cg} ${naive_output_cm} $outdir/ntoverview.txt $outdir/ntsum.txt 2>&1
 fi
 
 echo "</table>" >> $outdir/base_overview.html
 
 echo "---------------- Done! ----------------"
+echo "---------------- Done! ----------------<br />" >> $log
 
+mv $log $outdir/log.html
+
+mv $outdir/index.html $log
+