r182: added calldiff

lh3 · Oct 22, 2023 · 14015f9 · 14015f9
1 parent ab49529
commit 14015f9
Show file tree

Hide file tree

Showing 3 changed files with 48 additions and 6 deletions.
diff --git a/graph.c b/graph.c
@@ -26,6 +26,9 @@ void pg_post_process(const pg_opt_t *opt, pg_data_t *d)
 		fprintf(stderr, "[M::%s::%s] genome[%d]: %s; %d hits remain, of which %d are shadowed\n",
 				__func__, pg_timestamp(), j, g->label, tot, n_shadow);
 	}
+	for (i = j = 0; i < d->n_gene; ++i)
+		if (d->gene[i].preferred) ++j;
+	fprintf(stderr, "[M::%s] there are %d preferred genes\n", __func__, j);
 }
 
 pg_graph_t *pg_graph_init(pg_data_t *d)

diff --git a/pangene.h b/pangene.h
@@ -3,7 +3,7 @@
 
 #include <stdint.h>
 
-#define PG_VERSION "0.0-r177-dirty"
+#define PG_VERSION "0.0-r182-dirty"
 
 #define PG_F_WRITE_BED_RAW      0x1
 #define PG_F_WRITE_BED_WALK     0x2

diff --git a/pangene.js b/pangene.js
@@ -797,11 +797,16 @@ function pg_cmd_call(args) {
 }
 
 function pg_cmd_call2html(args) {
-	let endpoint = "/view";
-	for (const o of getopt(args, "", [])) {
+	let endpoint = "/view", graph = null;
+	for (const o of getopt(args, "e:g:", [])) {
+		if (o.opt == "-e") endpoint = o.arg;
+		else if (o.opt == "-g") graph = o.arg;
 	}
 	if (args.length == 0) {
-		print("Usage: pangene.js call2html <pangene-call.out>");
+		print("Usage: pangene.js call2html [options] <pangene-call.out>");
+		print("Options:");
+		print(`  -e STR     endpoint [${endpoint}]`);
+		print(`  -g STR     graph name []`);
 		return;
 	}
 	print(`<head>`);
@@ -822,7 +827,9 @@ function pg_cmd_call2html(args) {
 		const st = (t[4][0] == ">"? "&gt;" : "&lt;") + t[4].substr(1);
 		const en = (t[5][0] == ">"? "&gt;" : "&lt;") + t[5].substr(1);
 		const genes = [t[4].substr(1), t[8], t[5].substr(1)].join(",");
-		const link = `${endpoint}?gene=${genes}&step=0&ori=` + t[4].substr(1);
+		let link = `${endpoint}?`;
+		if (graph != null) link += `graph=${graph}&`;
+		link += `gene=${genes}&step=0&ori=` + t[4].substr(1);
 		const gene_space = t[8].replace(/,/g, ", ");
 		let out = `<tr><td style="text-align: right;">${t[1]}<td style="text-align: right;">${t[2]}<td style="text-align: right;">${t[6]}`;
 		out += `<td style="white-space: nowrap;"><a href="${link}" target="_blank">${st} &mdash; ${en}</a><td>${gene_space}</tr>`;
@@ -832,6 +839,36 @@ function pg_cmd_call2html(args) {
 	print(`</body>`);
 }
 
+function pg_cmd_calldiff(args) {
+	for (const o of getopt(args, "", [])) {
+	}
+	if (args.length < 2) {
+		print("Usage: pangene.js calldiff <call1.out> <call2.out>");
+		return;
+	}
+	let h = {};
+	for (const line of k8_readline(args[0])) {
+		let t = line.split("\t");
+		if (t[0] != "BB") continue;
+		const g1 = t[4].substr(1), g2 = t[5].substr(1);
+		const key = g1 < g2? `${g1}\t${g2}` : `${g2}\t${g1}`;
+		h[key] = [false, t.slice(1).join("\t")];
+	}
+	for (const line of k8_readline(args[1])) {
+		let t = line.split("\t");
+		if (t[0] != "BB") continue;
+		const g1 = t[4].substr(1), g2 = t[5].substr(1);
+		const key = g1 < g2? `${g1}\t${g2}` : `${g2}\t${g1}`;
+		if (key in h)
+			h[key][0] = true;
+		else
+			print("B2", t.slice(1).join("\t"));
+	}
+	for (const key in h)
+		if (h[key][0] == false)
+			print("B1", h[key][1]);
+}
+
 function pg_cmd_parse_ensembl(args) {
 	let species = null;
 	for (const o of getopt(args, "s:", [])) {
@@ -951,16 +988,18 @@ function main(args)
 		print("Commands:");
 		print("  call           call variants from a pangene graph");
 		print("  call2html      generate a HTML page from call output");
+		print("  calldiff       compare two call files");
 		print("  parse-ensembl  generate protein files from Ensembl annotations");
 		print("  gfa2matrix     generate gene_presence_absence.Rtab");
-		print("  flt-mmseqs     drop redundant proteins");
+		//print("  flt-mmseqs     drop redundant proteins");
 		//print("  parse-gfa     parse a GFA file (for debugging only)");
 		exit(1);
 	}
 
 	var cmd = args.shift();
 	if (cmd == 'call') pg_cmd_call(args);
 	else if (cmd == 'call2html') pg_cmd_call2html(args);
+	else if (cmd == 'calldiff') pg_cmd_calldiff(args);
 	else if (cmd == 'parse-ensembl') pg_cmd_parse_ensembl(args);
 	else if (cmd == 'gfa2matrix') pg_cmd_gfa2matrix(args);
 	else if (cmd == 'flt-mmseqs') pg_cmd_flt_mmseqs(args);