# HG changeset patch
# User davidvanzessen
# Date 1466494370 14400
# Node ID e6bc976760d4f6adea46a5c3f5c76539cdf1c44d
# Parent 3cffb8a38bb1116d09d0f991246f30333cc9bc98
Uploaded
diff -r 3cffb8a38bb1 -r e6bc976760d4 aa_histogram.r
--- 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()
diff -r 3cffb8a38bb1 -r e6bc976760d4 sequence_overview.r
--- 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("
", file=main.html, append=F)
cat("CDR1+FR2+CDR2+FR3+CDR3 sequences that show up more than once", file=main.html, append=T)
-cat("Sequence | Functionality | ca1 | ca2 | cg1 | cg2 | cg3 | cg4 | cm |
", file=main.html, append=T)
+cat("Sequence | Functionality | ca1 | ca2 | cg1 | cg2 | cg3 | cg4 | cm | un |
", 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, "
"), file=sequence.id.page, append=T)
+ }
}
cat("
", 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
diff -r 3cffb8a38bb1 -r e6bc976760d4 tmp/igat.r
--- 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),]
}
diff -r 3cffb8a38bb1 -r e6bc976760d4 wrapper.sh
--- 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 ----------------
" > $output
+echo "---------------- read parameters ----------------
" > $log
echo "unpacking IMGT file"
@@ -56,7 +57,7 @@
echo "${BLASTN_DIR}"
echo "---------------- identification ($method) ----------------"
-echo "---------------- identification ($method) ----------------
" >> $output
+echo "---------------- identification ($method) ----------------
" >> $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 ----------------
" >> $output
+echo "---------------- merge_and_filter.r ----------------
" >> $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 ----------------
" >> $output
+echo "---------------- creating new IMGT zip ----------------
" >> $log
mkdir $outdir/new_IMGT
@@ -132,7 +133,7 @@
cd $tmp
echo "---------------- mutation_analysis.r ----------------"
-echo "---------------- mutation_analysis.r ----------------
" >> $output
+echo "---------------- mutation_analysis.r ----------------
" >> $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 ----------------
" >> $output
+echo "---------------- mutation_analysis.py ----------------
" >> $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 ----------------
" >> $log
mkdir $outdir/sequence_overview
@@ -168,7 +169,6 @@
echo "$ID | $seq | $class | $A | $C | $G | $T |
" >> $outdir/base_overview.html
done < $outdir/sequence_overview/ntoverview.txt
-
echo "$title
" > $output
#display the matched/unmatched for clearity
@@ -185,10 +185,12 @@
echo "Percentage unmatched: ${perc_count}
" >> $output
echo "---------------- main tables ----------------"
+echo "---------------- main tables ----------------
" >> $log
for func in ${funcs[@]}
do
echo "---------------- $func table ----------------"
+ echo "---------------- $func table ----------------
" >> $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 ----------------
" >> $log
echo "unmatched
" >> $output
@@ -247,6 +250,7 @@
echo "---------------- images ----------------"
+echo "---------------- images ----------------
" >> $log
echo "
" >> $output
echo "download data
" >> $output
@@ -297,44 +301,67 @@
echo "" >> $output
echo "---------------- baseline ----------------"
+echo "---------------- baseline ----------------
" >> $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 ----------------
" >> $log
if [[ "$naive_output" != "None" ]]
then
echo "---------------- imgt_loader.r ----------------"
+ echo "---------------- imgt_loader.r ----------------
" >> $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 ----------------
" >> $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 "" >> $outdir/base_overview.html
echo "---------------- Done! ----------------"
+echo "---------------- Done! ----------------
" >> $log
+mv $log $outdir/log.html
+
+mv $outdir/index.html $log
+