merge dev into main (#1)

* populate the bin folder

* updated conda environment

* populate modules folder

* remove need for DBVERSION

* added README

* fixed package versioning and bumped to latest

* improved filter_best_bitscore.py

* updated filter_best_bitscore process

* updated extra info parsing

Author: jpalmer37

.gitignore

@@ -0,0 +1 @@
+bin/__pycache__

README.md

@@ -0,0 +1,154 @@
+# its-nf
+
+Prepare a report for taxonomic assignment based on [ITS](https://en.wikipedia.org/wiki/Internal_transcribed_spacer) sequences, using [BLAST](https://blast.ncbi.nlm.nih.gov/Blast.cgi).
+
+## Usage
+
+The pipeline requires a list of BLAST databases to run against. It should follow the following format:
+
+```csv
+ID,DBNAME,PATH
+ncbi,its_ncbi,/path/to/ncbi/2024-05-16_its_ncbi
+unite,its_unite,/path/to/unite/2024-05-13_its_unite
+```
+
+...where we expect to find the actual database files at:
+
+```
+/path/to/ncbi/2024-05-16_its_ncbi/its_ncbi.ndb
+/path/to/ncbi/2024-05-16_its_ncbi/its_ncbi.nhr
+/path/to/ncbi/2024-05-16_its_ncbi/its_ncbi.nin
+...etc
+/path/to/unite/2024-05-13_its_unite/its_unite.ndb
+/path/to/unite/2024-05-13_its_unite/its_unite.nhr
+/path/to/unite/2024-05-13_its_unite/its_unite.nin
+...etc
+```
+
+The pipeline also assumes that there is a `metadata.json` file alongside the database files
+
+```
+/path/to/ncbi/2024-05-16_its_ncbi/metadata.json
+/path/to/unite/2024-05-13_its_unite/metadata.json
+```
+
+The contents of the metadata file may vary by database, but we assume that:
+
+- The file contains a single top-level object (not an array or atomic value).
+- The top-level object includes these fields:
+
+```
+version
+date
+```
+
+The values associated with those fields will be incorporated into the blast results. All other fields in
+the `metadata.json` file are ignored.
+
+```
+nextflow run BCCDC-PHL/its-nf \
+  --databases </path/to/blast/databases.csv> \
+  --taxonkit_db </path/to/taxonkit/database/> \
+  --fasta_input </path/to/fasta_dir> \
+  --outdir </path/to/output_dir>
+```
+
+By default, minimum identity and coverage thresholds of 95% will be applied to the blast results.
+Alternate thresholds can be applied using the `--minid` and `--mincov` flags.
+
+```
+nextflow run BCCDC-PHL/its-nf \
+  --databases </path/to/blast/databases.csv> \
+  --taxonkit_db </path/to/taxonkit/database/> \
+  --fasta_input </path/to/fasta_dir> \
+  --minid 99.0 \
+  --mincov 97.5 \
+  --outdir </path/to/output_dir>
+```
+
+Collecting database metadata from the `metadata.json` file can be skipped using the `--no_db_metadata` flag.
+
+```
+nextflow run BCCDC-PHL/its-nf \
+  --databases </path/to/blast/databases.csv> \
+  --taxonkit_db </path/to/taxonkit/database/> \
+  --no_db_metadata \
+  --fasta_input </path/to/fasta_dir> \
+  --outdir </path/to/output_dir>
+```
+
+
+## Outputs
+
+Each sequence will have a separate output directory, named using the seq ID parsed from
+the fasta header. That directory will contain:
+
+```
+<seq_id>_<db_id>_blast.csv
+<seq_id>_<db_id>_blast_best_bitscore.csv
+<seq_id>_<db_id>_blast_filtered.csv
+<seq_id>_<db_id>_lineages.tsv
+<seq_id>_<db_id>_seq_qc.csv
+```
+
+The `_blast.csv`, `_blast_filtered.csv` and `blast_best_bitscore.csv` files have the following headers:
+
+```
+query_seq_id
+subject_accession
+subject_strand
+query_length
+query_start
+query_end
+subject_length
+subject_start
+subject_end
+alignment_length
+percent_identity
+percent_coverage
+num_mismatch
+num_gaps
+e_value
+bitscore
+subject_taxids
+subject_names
+genus
+species
+database_name
+database_version
+database_date
+```
+
+...though if the `--no_db_metadata` flag is used when running the pipeline, the last three fields will be omitted.
+
+The `_blast_best_bitscore.csv` file will only include one entry per species per database if there are multiple matches from the same
+species with equally-good bitscores.
+
+The `_lineages.tsv` file is generated by `taxonkit`, and has the following headers:
+
+```
+query_taxid
+lineage
+lineage_taxids
+query_taxon_name
+lineage_ranks
+```
+
+...where the `lineage`, `lineage_taxids`, and `lineage_ranks` are themselves semicolon-separated lists.
+
+The `seq_qc.csv` file has the following headers:
+
+```
+seq_length
+num_ambiguous_bases
+num_n_bases
+```
+
+There will also be collected ouputs in the top-level of the `--outdir` directory, named:
+
+```
+collected_blast.csv
+collected_blast_best_bitscore.csv
+```
+
+...which will include results from all sequences.

bin/add_db_metadata.py

@@ -0,0 +1,61 @@
+#!/usr/bin/env python3
+
+import argparse
+import csv
+import json
+import sys
+
+
+def parse_metadata(metadata_path):
+    with open(metadata_path, 'r') as f:
+        metadata = json.load(f)
+
+    return metadata
+
+
+def parse_blast_results(blast_results_path):
+    header_fieldnames = []
+    blast_results = []
+
+    with open(blast_results_path, 'r') as f:
+        header_line = f.readline().strip()
+        header_fieldnames = header_line.split(',')
+
+    with open(blast_results_path, 'r') as f:
+        reader = csv.DictReader(f)
+        for row in reader:
+            blast_results.append(row)
+
+    return header_fieldnames, blast_results
+
+
+def main(args):
+    metadata = parse_metadata(args.metadata)
+    output_fieldnames, blast_results = parse_blast_results(args.blastresult)
+    for record in blast_results:
+        record['database_name'] = args.database_name
+        record['database_version'] = metadata.get('version', None)
+        record['database_date'] = metadata.get('date', None)
+
+
+    output_fieldnames += [
+        'database_name',
+        'database_version',
+        'database_date',
+    ]
+
+    writer = csv.DictWriter(sys.stdout, fieldnames=output_fieldnames, delimiter=',', lineterminator='\n', quoting=csv.QUOTE_MINIMAL, extrasaction='ignore')
+
+    writer.writeheader()
+    for record in blast_results:
+        writer.writerow(record)
+
+
+
+if __name__ == "__main__":
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-m','--metadata')
+    parser.add_argument('-d','--database_name')
+    parser.add_argument('-b','--blastresult')
+    args = parser.parse_args()
+    main(args)

bin/bind_taxonkit.py

@@ -1,44 +1,100 @@
 #!/usr/bin/env python3
 
 import argparse
-import pandas as pd 
-import numpy as np
+import csv
+import json
+import sys
 
-def main(args):
-    cols = 'subject_taxids,lineage,name,rank'.split(',')
-    #taxon_results = pd.read_csv('F1910235_taxon_results.txt',sep = '\t', names=cols)
-    taxon_results = pd.read_csv(args.taxonresult,sep = '\t', names=cols)
 
-    blast_results = pd.read_csv(args.blastresult,sep = ',')
-    taxon_results = taxon_results.dropna()
+def parse_taxonkit_lineage(taxonkit_path):
+    """
+    Parse taxonkit lineage outputs
 
+    :param taxonkit_path: path to taxonkit lineage output
+    :type taxonkit_path: str
+    :return: taxonkit lineages, by taxid. dict of dicts. keys of inner dicts: ['taxid', 'lineage', 'name', 'rank']
+    :rtype: dict[str, dict[str, str]]
+    """
+    taxonkit_lineage_by_taxid = {}
+    with open(taxonkit_path, 'r') as f:
+        reader = csv.DictReader(f, delimiter='\t')
+        for row in reader:
+            taxonkit_lineage_record = {}
+            query_taxid = row['query_taxid']
+            lineage = row['lineage']
+            lineage_split = lineage.split(';')
+            taxids = row['lineage_taxids']
+            taxids_split = taxids.split(';')
+            name = row['query_taxon_name']
+            ranks = row['lineage_ranks']
+            ranks_split = ranks.split(';')
 
+            taxonkit_lineage_record['query_taxid'] = query_taxid
+            for idx, rank in enumerate(ranks_split):
+                if rank == 'species':
+                    taxonkit_lineage_record['species_taxid'] = taxids_split[idx]
+                    taxonkit_lineage_record['species_name'] = lineage_split[idx]
+                elif rank == 'genus':
+                    taxonkit_lineage_record['genus_taxid'] = taxids_split[idx]
+                    taxonkit_lineage_record['genus_name'] = lineage_split[idx]
 
-    conditions = [
-        (taxon_results['rank'] == "genus"),
-        (taxon_results['rank'] == "species"),
-        (taxon_results['rank'] == "strain")
+            taxonkit_lineage_by_taxid[query_taxid] = taxonkit_lineage_record
+            
+    return taxonkit_lineage_by_taxid
 
-    ]
 
-    choices_species = [None, taxon_results['lineage'].apply(lambda x: x.split(';')[-1]), taxon_results['lineage'].apply(lambda x: x.split(';')[-2])]
-    choices_genus = [taxon_results['lineage'].apply(lambda x: x.split(';')[-1]), taxon_results['lineage'].apply(lambda x: x.split(';')[-2]),taxon_results['lineage'].apply(lambda x: x.split(';')[-3])]
+def parse_blast_results(blast_results_path):
+    """
+    Parse blast results
+
+    :param blast_results_path: path to blast results
+    :type blast_results_path: str
+    
+    """
+    header_fieldnames = []
+    blast_results = []
+    with open(blast_results_path, 'r') as f:
+        header_line = f.readline().strip()
+        header_fieldnames = header_line.split(',')
+
+    with open(blast_results_path, 'r') as f:
+        reader = csv.DictReader(f)
+        for row in reader:
+            blast_results.append(row)
+
+    return header_fieldnames, blast_results
 
-    taxon_results['species'] = np.select(conditions,choices_species, default = taxon_results['lineage'].apply(lambda x: x.split(';')[-1]))
-    taxon_results['genus'] = np.select(conditions,choices_genus, default = taxon_results['lineage'].apply(lambda x: x.split(';')[-2]))
+
+def main(args):
 
+    taxonkit_lineage_by_taxid = parse_taxonkit_lineage(args.taxonresult)
+
+    output_fieldnames, blast_results = parse_blast_results(args.blastresult)
 
-    merged = pd.merge(blast_results,taxon_results[['subject_taxids','species','genus']],on='subject_taxids', how='left')
-    fil = merged['species'].str.contains('uncultured')
-    filtered_merged = merged[~fil]
+    for blast_result in blast_results:
+        subject_taxid = blast_result['subject_taxids']
+        if subject_taxid in taxonkit_lineage_by_taxid:
+            blast_result['species'] = taxonkit_lineage_by_taxid[subject_taxid].get('species_name', None)
+            blast_result['genus'] = taxonkit_lineage_by_taxid[subject_taxid].get('genus_name', None)
+        else:
+            blast_result['species'] = None
+            blast_result['genus'] = None
+    
+    output_fieldnames += [
+        'genus',
+        'species',
+    ]
+    writer = csv.DictWriter(sys.stdout, fieldnames=output_fieldnames, delimiter=',', lineterminator='\n', quoting=csv.QUOTE_MINIMAL, extrasaction='ignore')
+    writer.writeheader()
+    for blast_result in blast_results:
+        writer.writerow(blast_result)
+    
 
-    filtered_merged.to_csv(args.outfile)
 
 if __name__ == "__main__":
     parser = argparse.ArgumentParser()
 
     parser.add_argument('-f','--taxonresult')
     parser.add_argument('-b','--blastresult')
-    parser.add_argument('-o','--outfile')
     args = parser.parse_args()
-    main(args)
+    main(args)