Stats

inumanag · inumanag · commit a20e09b3ea94 · 2018-01-30T01:08:00.000-05:00
diff --git a/README.md b/README.md
@@ -1,5 +1,7 @@
 # SEDEF: SEgmental Duplication Evaluation Framework
 
+# UNDER CONSTRUCTION
+
 ## Compile
 
 Simple! Use
@@ -86,7 +88,7 @@ wc -l out.*bed
 >>  1656305 out/out.bed
 >>  1558896 out/out.init.bed
 >>   231472 out/out.final.bed
->>	  71204 out.hg19.bed
+>>	  67467 out.hg19.bed
 ```
 
 Then analyse:
@@ -101,90 +103,3 @@ done
 >> 0m 58s (7.5G)
 ```
 
-
-Misses & partials:
-
-+-------------------+------------+-----+----------------+-------+
-+ Stage             + Misses     + Kbp + Partial misses + Kbp   +
-+-------------------+------------+-----+----------------+-------+
-+ Seed              + 167 (0.7%) + 652 + 1,348 (5.5%)   + 6,163 +
-+ Merge & Extend    + 133 (0.5%) + 505 +    73 (0.3%)   +   384 +
-+ Chain             + 135 (0.6%) + 515 +   215 (0.9%)   + 1,536 +
-+-------------------+------------+-----+----------------+-------+
-
-
-### Mouse
-
-```bash
-for i in `seq 1 19` X Y; do 
-	for j in `seq 1 19` X Y; do  
-		SI=`awk '$1=="chr'$i'" {print $2}' mm8.fa.fai`; 
-		SJ=`awk '$1=="chr'$j'" {print $2}' mm8.fa.fai`; 
-		if [ "$SI" -le "$SJ" ] ; then 
-			for m in y n ; do
-			echo "~/mesa ./sedef search single mm8.fa chr$i chr$j $m >mouse_out/${i}_${j}_${m}.bed 2>mouse_out/log/${i}_${j}_${m}.log"
-			done; 
-		fi
-	done
-done | time parallel --will-cite -j 80 --eta
->> 12m 17s
-grep Total mouse_out/log/*.log | wc -l
->> 462
-grep Wall mouse_out/log/*.log | tr -d '(' | awk '{s+=$4}END{print s}'
->> 46962.9s (13.05 h)
-~/mesa ./sedef align bucket mouse_out mouse_out/bins 1000
->> 
-for j in mouse_out/bins/bucket_???? ; do
-	k=$(basename $j);
-	echo "~/mesa ./sedef align generate mm8.fa $j 11 >${j}.bed 2>mouse_out/log/bins/${k}.log"
-done | time parallel --will-cite -j 80 --eta
->> 
-grep Finished mouse_out/log/bins/*.log | wc -l
->> 1000
-grep Wall mouse_out/log/bins/*.log | tr -d '(' | awk '{s+=$4}END{print s}'
->> 236095 (65.58 h)cv
-cat mouse_out/bins/bucket_???? > mouse_out.init.bed
-cat mouse_out/bins/*.bed > mouse_out.final.bed
-wc -l mouse_out.*bed
->>   975511 mouse_out.final.bed
->>  1656305 mouse_out.init.bed
-
-# Merge logs and remove progress bars and headers
-
-for S in G S X Y ; do
-	for i in ${S}_*.e*; do grep -v '%' $i | tail -n+4 ; done > ${S}.log
-	cat ${S}_*.o* > ${S}.bed
-	rm -rf ${S}*.[eo]*
-done
-
-for i in G X Y S ; do 
-	mkdir -p bins/${i}
-	sedef/sedef align bucket ${i}.bed bins/${i} 2000 
-done
-
-for i in G X Y S ; do 
-	for j in bins/${i}/* ; do
-		k=$(basename $j);
-		echo "qsub -cwd -V -b y -N \"${i}_${k}\" -l h_vmem=10G -l h_rt=24:00:00 -l h_stack=8M " \
-			"python2.7 mesa sedef/sedef align generate fasta/hg19.fa $j"
-	done
-done
-
-for S in G X Y S ; do 
-	for i in ${S}_*.e*; do grep -v '\.\.' $i | tail -n+4 ; done > ${S}.align.log
-	cat ${S}_*.o* > ${S}.align.bed
-	rm -rf ${S}_*.[eo]*
-done
-
-# Count stuff
-for S in G S X Y ; do 
-	SEA=`cat ${S}.log | grep 'Wall' | awk '{print $4}' | tr -d '[()]' | awk '{s+=$1} END{print s}'`
-	ALN=`cat ${S}.align.log | grep 'Wall' | awk '{print $4}' | tr -d '[()]' | awk '{s+=$1} END{print s}'`
-	printf "%s: %5.1f %5.1f %5.1f\n" $S $((SEA/3600)) $((ALN/3600)) $(((SEA+ALN)/3600))
-done
-
-#mkdir -p output/search
-#zmv '(sedefrun_*).o*' 'output/search/$1.bed'
-```
-
-
diff --git a/src/stats_main.cc b/src/stats_main.cc
@@ -8,6 +8,8 @@
 #include <vector>
 #include <algorithm>
 #include <chrono>
+#include <bitset>
+#include <unordered_set>
 
 #include "align.h"
 #include "common.h"
@@ -164,13 +166,102 @@ void stats(const string &ref_path, const string &bed_path)
 
 /******************************************************************************/
 
-void stats_main(int argc, char **argv) 
+#include <boost/dynamic_bitset.hpp>
+
+void get_differences()
 {
-	if (argc < 2) {
-		eprn("invalid usage");
-		return;
+	map<string, boost::dynamic_bitset<>> sedef;
+	map<string, boost::dynamic_bitset<>> wgac;
+
+	string s;
+	ifstream fin("results/out.hg19.bed", ifstream::in);
+	while (getline(fin, s)) {
+		string cigar;
+		Hit h = Hit::from_bed(s, &cigar);
+
+		auto c1 = fmt::format("{}{}", h.query->name, "+-"[h.query->is_rc]);
+		auto c2 = fmt::format("{}{}", h.ref->name, "+-"[h.ref->is_rc]);
+		if (sedef.find(c1)==sedef.end()) sedef[c1]=boost::dynamic_bitset<>(250000000);
+		if (sedef.find(c2)==sedef.end()) sedef[c2]=boost::dynamic_bitset<>(250000000);
+		for (int i = h.query_start; i < h.query_end; i++) sedef[c1].set(i);
+		for (int i = h.ref_start; i < h.ref_end; i++) sedef[c2].set(i);
+	}
+
+	eprn("sedef done");
+
+	ifstream fiw("data/GRCh37GenomicSuperDup.tab");
+	getline(fiw, s);
+	unordered_set<string> seen;
+	while (getline(fiw, s)) {
+		Hit h = Hit::from_wgac(s);
+		auto c1 = fmt::format("{}{}", h.query->name, "+-"[h.query->is_rc]);
+		auto c2 = fmt::format("{}{}", h.ref->name, "+-"[h.ref->is_rc]);		
+		if (c1.size() > 6 || c2.size() > 6)
+			continue;
+		
+		if (seen.find(h.name) == seen.end()) {
+			seen.insert(h.name);
+			if (wgac.find(c1)==wgac.end()) wgac[c1]=boost::dynamic_bitset<>(250000000);
+			if (wgac.find(c2)==wgac.end()) wgac[c2]=boost::dynamic_bitset<>(250000000);
+			for (int i = h.query_start; i < h.query_end; i++) wgac[c1].set(i);
+			for (int i = h.ref_start; i < h.ref_end; i++) wgac[c2].set(i);
+		}
+	}
+
+	eprn("wgac done");
+
+	FastaReference fr("data/hg19/hg19.fa");
+
+	int intersect = 0, wgac_only = 0, wgac_span = 0, sedef_only = 0, sedef_span = 0;
+
+	int sedef_extra_upper = 0;
+	int miss_upper = 0;
+
+	for (auto &p: sedef) {
+		auto &s = p.second;
+		auto &w = wgac[p.first];
+
+		auto seq = fr.get_sequence(p.first.substr(0, p.first.size()-1));
+
+		for (int i = 0; i < seq.size(); i++) {
+			if ((s[i] & (~w[i])) && isupper(seq[i]) && seq[i] != 'N') {
+				sedef_extra_upper++;
+			}
+			if ((w[i] & (~s[i])) && isupper(seq[i]) && seq[i] != 'N') {
+				miss_upper++;
+			}
+		}
+
+		intersect += (s & w).count();
+		wgac_only += (w & (~s)).count();
+		sedef_only += (s & (~w)).count();
+		sedef_span += s.count();
+		wgac_span += w.count();
+	}
+	
+	eprn("SEDEF: span {:12n}\n"
+		 "       only {:12n}\n"
+		 "       on/u {:12n}\n"
+		 "       miss {:12n}\n"
+		 "       mi/u {:12n}\n"
+		 "WGAC:  span {:12n}\n"
+		 "       intr {:12n}", sedef_span, sedef_only, sedef_extra_upper, wgac_only, miss_upper, wgac_span, intersect);
+}
+
+/******************************************************************************/
+
+void stats_main(int argc, char **argv)
+{
+	if (argc < 3) {
+		throw fmt::format("Not enough arguments to stats");
+	}
+
+	string command = argv[0];
+	if (command == "generate") {
+		stats(argv[1], argv[2]);
+	} else if (command == "diff") {
+		get_differences(); //(argv[1], argv[2], atoi(argv[3]));
+	} else {
+		throw fmt::format("Unknown stats command");
 	}
-	string ref_path = argv[0];
-	string bed_path = argv[1];
-	stats(ref_path, bed_path);
 }
