Added consensus rank determination

Author: gregdenay

.gitmodules

@@ -4,3 +4,6 @@
 [submodule "taxidTools"]
 	path = taxidTools
 	url = https://github.com/CVUA-RRW/taxidTools
+[submodule "scripts/taxidTools"]
+	path = scripts/taxidTools
+	url = https://github.com/CVUA-RRW/taxidTools

Snakefile

@@ -317,6 +317,9 @@ rule get_consensus_level:
         cons = "reports/consensus.tsv"
     message:
         "Determining consensus ranks"
+    params:
+        lineage = config["taxonomy"]["rankedlineage_dmp"],
+        nodes = config["taxonomy"]["nodes_dmp"],
     script:
         "../scripts/consensus_levels.py"
 
@@ -330,6 +333,7 @@ rule write_report:
         nderep = "reports/derep_number.txt",
         nNfilt = "reports/high_N.txt",
         clusterSize = "reports/cluster_size.tsv"
+        consensus = "reports/consensus.tsv"
     output:
         "reports/report.html"
     params:

scripts/consensus_levels.py

@@ -0,0 +1,36 @@
+#!/usr/bin/env python3
+# -*- coding: utf-8 -*-
+
+import pandas as pd
+from taxidTools.taxidTools import Taxdump
+
+def get_consensus_rank(txd, idlist):
+	lca = txd.lowestCommonNode([str(t) for t in idlist])
+	return(txd.getRank(lca))
+
+def main(table, outfile, lineage, nodes):
+	df = pd.read_csv(table, sep = '\t')
+	
+	df["entries"] = df["query"] + df["query_size"]
+	
+	dfout = pd.DataFrame()
+	txd = Taxdump(lineage, nodes)
+	
+	for entry in set(df["entries"]):
+		sub = df[df["entries"] == entry].reindex()
+		
+		new = sub.head(1)[["query_name", "query_taxid", "query_size", "query_relsize"]]
+		
+		for i in range(0,4):
+			taxids = list(sub[sub["distance"] <= i]["target_taxid"])
+			taxids.extend(list(sub.head(1)["query_taxid"]))
+			rank = get_consensus_rank(txd, taxids)
+			
+			new[f"Consensus rank with {i} mismatches"] = rank
+		
+		dfout = dfout.append(new, ignore_index=True)
+		
+		dfout.to_csv(outfile, sep = '\t')
+
+if __name__ == '__main__':
+	main(snakemake.input[0], snakemake.output[0], snakemake.params["lineage"], snakemake.params["nodes"])

scripts/taxidTools

@@ -40,6 +40,7 @@ derep <- snakemake@input[["derep"]]
 nderep <- snakemake@input[["nderep"]]
 nNfilt <- snakemake@input[["nNfilt"]]
 clusterSize <- snakemake@input[["clusterSize"]]
+consensus <- snakemake@input[["consensus"]]
 
 trimming <- snakemake@params[["trimming"]]
 database <- snakemake@params[["database"]]
@@ -70,6 +71,37 @@ primers <- snakemake@params[["primers"]]
 cat("* Primer file: ", primers) 
 ```
 
+## Consensus ranks
+
+Achievable Taxonomic resolution for each sequence cluster at different distance (number of mismatch) levels.
+
+```{r consensus}
+data_table <- read.csv(file = consensus, sep = '\t', check.names = FALSE)
+
+datatable(data_table, filter = list(position = 'top'), rownames = FALSE, escape = FALSE,
+		extensions= c("ColReorder", "Buttons"),
+		selection = 'none',
+		options = list(dom = 'BRrlftpi',
+						autoWidth=FALSE, 
+						scrollX = TRUE,
+						lengthMenu = list(c(50, 100, 500, -1), 
+										  c('50', '100', '500', 'All')),
+						ColReorder = TRUE,
+						buttons = list('copy',
+										'print',
+										list(
+										extend = 'collection',
+										buttons = c('csv', 'excel', 'pdf'),
+										text = 'Download'
+										),
+										I('colvis')
+										),
+						deferRender = TRUE,
+						scroller = TRUE
+						)
+		)
+```
+
 ## Pairwise distances