first commit

Author: Sherrie Wang

bin/bind_taxonkit.py

@@ -0,0 +1,44 @@
+#!/usr/bin/env python3
+
+import argparse
+import pandas as pd 
+import numpy as np
+
+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()
+
+
+
+    conditions = [
+        (taxon_results['rank'] == "genus"),
+        (taxon_results['rank'] == "species"),
+        (taxon_results['rank'] == "strain")
+
+    ]
+
+    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])]
+
+    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]))
+    
+
+    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]
+
+    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)

bin/filter_best_bitscore.py

@@ -0,0 +1,53 @@
+#!/usr/bin/env python3
+
+import argparse
+import csv
+import json
+import sys
+
+
+def parse_blast_report(blast_report_path):
+    with open(blast_report_path, 'r') as f:
+        header_line = f.readline().strip()
+        header_fieldnames = header_line.split(',')
+
+    int_fields = [
+        'bitscore',
+    ]
+    parsed_blast_report = []
+    with open(blast_report_path, 'r') as f:
+        reader = csv.DictReader(f)
+        for row in reader:
+            for field in int_fields:
+                row[field] = int(row[field])
+            parsed_blast_report.append(row)
+
+    return header_fieldnames, parsed_blast_report
+
+
+def determine_best_bitscore(parsed_blast_report):
+    best_bitscore = 0
+    for blast_record in parsed_blast_report:
+        if blast_record['bitscore'] > best_bitscore:
+            best_bitscore = blast_record['bitscore']
+
+    return best_bitscore
+        
+
+def main(args):
+    output_fieldnames, blast_report = parse_blast_report(args.input)
+    best_bitscore = determine_best_bitscore(blast_report)
+
+    filtered_blast_report = list(filter(lambda x: x['bitscore'] == best_bitscore, blast_report))
+
+    writer = csv.DictWriter(sys.stdout, fieldnames=output_fieldnames, dialect='unix', quoting=csv.QUOTE_MINIMAL)
+    writer.writeheader()
+    for row in filtered_blast_report:
+        writer.writerow(row)
+
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-i', '--input')
+    args = parser.parse_args()
+    main(args)

bin/seq_qc.py

@@ -0,0 +1,48 @@
+#!/usr/bin/env python3
+
+import argparse
+import json
+
+def parse_fasta(fasta_path):
+    fasta = {}
+    seq_id = ''
+    seq = []
+    with open(fasta_path, 'r') as f:
+        for line in f:
+            if line.startswith('>'):
+                seq_id = line.strip().lstrip('>').split()[0]
+            else:
+                seq.append(line.strip())
+
+    fasta['id'] = seq_id
+    fasta['seq'] = ''.join(seq)
+
+    return fasta
+
+def main(args):
+    iupac_ambiguous_bases = set([
+        'M', 'R', 'W', 'S', 'Y', 'K',
+        'V', 'H', 'D', 'B',
+    ])
+    fasta = parse_fasta(args.input)
+    seq_length = len(fasta['seq'])
+    num_ambiguous_bases = 0
+    num_n_bases = 0
+
+    for base in fasta['seq']:
+        if base.upper() in iupac_ambiguous_bases:
+            num_ambiguous_bases += 1
+
+    for base in fasta['seq']:
+        if base.upper() == 'N':
+            num_n_bases += 1
+
+    print('seq_length,num_ambiguous_bases,num_n_bases')
+    print(','.join(map(str, [seq_length, num_ambiguous_bases, num_n_bases])))
+    
+
+if __name__ == '__main__':
+    parser = argparse.ArgumentParser()
+    parser.add_argument('-i', '--input')
+    args = parser.parse_args()
+    main(args)

environments/environment.yml

@@ -0,0 +1,12 @@
+name: its-nf
+channels:
+  - conda-forge
+  - bioconda
+  - defaults
+dependencies:
+  - python=3
+  - blast=2.13.0
+  - taxonkit=0.14.1
+  - pandas=1.5.3
+  - biopython=1.81
+  - kraken2=2.1.3

main.nf

@@ -0,0 +1,31 @@
+#!/usr/bin/env nextflow
+
+nextflow.enable.dsl = 2
+
+include { seq_qc }               from './modules/blast.nf'
+include { blastn }               from './modules/blast.nf'
+include { filter_best_bitscore } from './modules/blast.nf'
+
+workflow {
+
+  if (params.samplesheet_input != 'NO_FILE') {
+    ch_fasta = Channel.fromPath(params.samplesheet_input).splitCsv(header: true).map{ it -> [it['ID'], it['FILE']] }
+  } else {
+    ch_fasta = Channel.fromPath(params.fasta_search_path)
+  }
+
+  ch_db = Channel.fromPath(params.db_dir).combine(Channel.of(params.db_name))
+
+  ch_seqs = ch_fasta.splitFasta(file: true)
+
+  main:
+    seq_qc(ch_seqs)
+    ch_blast = blastn(ch_seqs.combine(ch_db))
+    ch_best_blast = filter_best_bitscore(ch_blast.taxon_results)
+
+    ch_blast.taxon_results
+        .collectFile(it -> it[1], name: "combined_blast_species_genus_results.csv", storeDir: params.outdir,keepHeader: true, skip: 1)
+  
+    ch_best_blast.blast_best_bitscore_csv
+        .collectFile(it -> it[1], name: "combined_top1_results.csv", storeDir: params.outdir,keepHeader: true, skip: 1)
+}

modules/blast.nf

@@ -0,0 +1,89 @@
+process seq_qc {
+
+    tag { sample_id }
+
+    executor 'local'
+
+    publishDir "${params.outdir}/${sample_id}", mode: 'copy', pattern: "${sample_id}_seq_qc.csv"
+
+    input:
+    tuple path(seq)
+
+    output:
+    tuple val(sample_id), path("${sample_id}_seq_qc.csv"), emit: seq_qc_csv
+
+    script:
+    sample_id = seq.getName().split('\\.')[0]
+    """
+    seq_qc.py -i ${seq} > ${sample_id}_seq_qc.csv
+    """
+}
+
+process blastn {
+    //errorStrategy 'ignore'
+    
+
+    publishDir "${params.outdir}/${sample_id}", mode: 'copy', pattern: "${sample_id}*"
+
+    input:
+    tuple path(query), path(db_dir), val(db_name)
+
+    output:
+    tuple val(sample_id), path("${sample_id}_blast.csv"), emit: blast_csv, optional:true
+    tuple val(sample_id), path("${sample_id}_seq_description"), emit: seq_description, optional:true
+    tuple val(sample_id), path("${sample_id}_blast_species_genus_results.csv"), emit: taxon_results, optional:true
+    tuple val(sample_id), path("${sample_id}_taxon_results.txt"), emit: raw_taxon_results, optional:true
+    
+    script:
+    sample_id = query.getName().split('\\.')[0]
+    """
+    export BLASTDB="${db_dir}"
+
+    echo "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" > ${sample_id}_blast.csv
+
+    blastn \
+      -db ${db_name} \
+      -num_threads ${task.cpus} \
+      -perc_identity ${params.minid} \
+      -qcov_hsp_perc ${params.mincov} \
+      -query ${query} \
+      -outfmt "6 qseqid saccver sstrand qlen qstart qend slen sstart send length pident qcovhsp mismatch gaps evalue bitscore staxids sscinames" \
+    | tr \$"\\t" "," >> ${sample_id}_blast.csv
+
+    tail -qn+2 ${sample_id}_blast.csv | cut -d',' -f2 | sort -u > seqids
+    blastdbcmd -db ${db_name} -entry_batch seqids | grep '>' > ${sample_id}_seq_description
+
+
+
+
+    if [ "${db_dir}" == "2022-11-16_nt" ] || [ "${db_dir}" == "refseq_its" ] ; then
+        tail -qn+2 ${sample_id}_blast.csv | cut -d',' -f17 | sort -u > taxids
+        taxonkit lineage -r -n  taxids > ${sample_id}_taxon_results.txt
+        bind_taxonkit.py -f ${sample_id}_taxon_results.txt -b ${sample_id}_blast.csv -o ${sample_id}_blast_species_genus_results.csv
+    fi
+
+
+    """
+}
+
+
+
+process filter_best_bitscore {
+
+    tag { sample_id }
+
+    executor 'local'
+
+    publishDir "${params.outdir}/${sample_id}", mode: 'copy', pattern: "${sample_id}_blast_best_bitscore.csv"
+
+    input:
+    tuple val(sample_id), path(full_blast_report)
+
+    output:
+    tuple val(sample_id), path("${sample_id}_blast_best_bitscore.csv"), emit: blast_best_bitscore_csv
+
+    script:
+    """
+    filter_best_bitscore.py -i ${full_blast_report} > ${sample_id}_blast_best_bitscore.csv
+    """
+}

nextflow.config

@@ -0,0 +1,61 @@
+manifest {
+  author = 'Sherrie Wang, Dan Fornika'
+  name = 'BCCDC-PHL/its-nf'
+  version = '0.1.0'
+  description = 'BCCDC-PHL Taxonomic Assignment from ITS sequences'
+  mainScript = 'main.nf'
+  nextflowVersion = '>=20.01.0'
+}
+
+params {
+  profile = false
+  cache = ''
+  outdir = 'results'
+  fasta_exts = ['*.fa', '*.fasta', '*.fna']
+  fasta_search_path = makeFastaSearchPath(fasta_exts)
+  fasta_input = 'NO_FILE'
+  samplesheet_input = 'NO_FILE'
+  db_dir = 'NO_FILE'
+  db_name = 'db'
+  db_version = 'current'
+  minid = 95.0
+  mincov = 95.0
+  pipeline_short_name = parsePipelineName(manifest.toMap().get('name'))
+  pipeline_minor_version = parseMinorVersion(manifest.toMap().get('version'))
+}
+
+def makeFastaSearchPath (fasta_exts) {
+    def fasta_search_path = []
+    for (ext in fasta_exts) {
+        fasta_search_path.add(params.fasta_input.toString() + '/' + ext.toString())
+        fasta_search_path.add(params.fasta_input.toString() + '/**/' + ext.toString())
+    }
+    return fasta_search_path
+}
+
+def parseMinorVersion(version) {
+    minor_version = version.split('\\.')[0..1].join('.')
+    return minor_version
+}
+
+def parsePipelineName(name) {
+    short_name = name.split('/')[1]
+    return short_name
+}
+
+profiles {
+  conda {
+     process.conda = "$baseDir/environments/environment.yml"
+     if (params.cache){
+       conda.cacheDir = params.cache
+     }
+  }
+}
+
+process {
+
+  withName: blastn {
+    shell = ['/bin/bash', '-uo','pipefail' ]
+  }
+
+}