From 8c85c6c2662a69775208952c0f1e05c651a5b405 Mon Sep 17 00:00:00 2001 From: Jim Downie Date: Wed, 8 Oct 2025 10:10:24 +0100 Subject: [PATCH 01/13] feat(ncbi): replace deprecated ncbi-datasets-pylib library - use lightweight custom class reimplementing required functionality with the NCBI API - tidy up some logic with counting loops by just returning counts from the API - use a bit of pathlib to do some path stuff - re-format modified scripts with Ruff --- envs/dataset.yaml | 3 +- scripts/DetermineGenera.py | 517 +++++++++++++++---------- scripts/FetchGenomesRefSeq.py | 234 +++++++---- scripts/FetchGenomesRefSeqRelatives.py | 307 +++++++++------ scripts/NCBIApiTools.py | 119 ++++++ 5 files changed, 774 insertions(+), 406 deletions(-) create mode 100644 scripts/NCBIApiTools.py diff --git a/envs/dataset.yaml b/envs/dataset.yaml index 3943565..ba5e685 100644 --- a/envs/dataset.yaml +++ b/envs/dataset.yaml @@ -6,5 +6,4 @@ dependencies: - python=3.9 - pip=21.2.1 - kraken2 - - pip: - - ncbi-datasets-pylib==12.6.0 \ No newline at end of file + - requests=2.32 diff --git a/scripts/DetermineGenera.py b/scripts/DetermineGenera.py index 3bb1734..c86c1df 100644 --- a/scripts/DetermineGenera.py +++ b/scripts/DetermineGenera.py @@ -1,63 +1,117 @@ from __future__ import division + import argparse -import configparser -import json import os -import sys - -try: - import ncbi.datasets -except ImportError: - print('ncbi.datasets module not found. To install, run `pip install ncbi-datasets-pylib`.') -from ncbi.datasets.package import dataset +from NCBIApiTools import NcbiApi parser = argparse.ArgumentParser() -parser.add_argument("-i", type=str, action='store', dest='tax', metavar='TAX',help='define tax genus file') -parser.add_argument("-t", type=str, action='store', dest='type', metavar='TYPE',help='define ranking type') -parser.add_argument("-na", type=str, action='store', dest='namesfile', metavar='NAMES',help='NCBI names.dmp') -parser.add_argument("-no", type=str, action='store', dest='nodesfile', metavar='NODES',help='NCBI nodes.dmp') -parser.add_argument("-suf", type=str, action='store', dest='suffix', metavar='SUFFIX',help='suffix outputfile') -parser.add_argument("-od", type=str, action='store', dest='outdir', metavar='OUTDIR',help='output directory') -parser.add_argument("-g", type=str, action='store', dest='spoi', metavar='SPOI',help='species of interest') -parser.add_argument('--version', action='version', version='%(prog)s 1.0') +parser.add_argument( + "-i", + type=str, + action="store", + dest="tax", + metavar="TAX", + help="define tax genus file", +) +parser.add_argument( + "-t", + type=str, + action="store", + dest="type", + metavar="TYPE", + help="define ranking type", +) +parser.add_argument( + "-na", + type=str, + action="store", + dest="namesfile", + metavar="NAMES", + help="NCBI names.dmp", +) +parser.add_argument( + "-no", + type=str, + action="store", + dest="nodesfile", + metavar="NODES", + help="NCBI nodes.dmp", +) +parser.add_argument( + "-suf", + type=str, + action="store", + dest="suffix", + metavar="SUFFIX", + help="suffix outputfile", +) +parser.add_argument( + "-od", + type=str, + action="store", + dest="outdir", + metavar="OUTDIR", + help="output directory", +) +parser.add_argument( + "-g", + type=str, + action="store", + dest="spoi", + metavar="SPOI", + help="species of interest", +) +parser.add_argument( + "-k", + type=str, + action="store", + dest="key", + metavar="NCBI API KEY", + help="NCBI API Key", + required=False, + default=os.environ.get("NCBI_API_KEY"), +) +parser.add_argument("--version", action="version", version="%(prog)s 1.0") args = parser.parse_args() + def readNames(names_tax_file): - ''' + """ input: - name.dmp (NCBI Taxonomy) output: - dictionary of form {node: name} - dictionary of form {sci name: node} - ''' + """ tax_names = {} - tax_names_reverse= {} + tax_names_reverse = {} tax_names_sci = {} - multiple_names={} - with open(names_tax_file, 'r') as nodes_tax: + multiple_names = {} + with open(names_tax_file, "r") as nodes_tax: for line in nodes_tax: - node = [field.strip() for field in line.split('|')] - if 'scientific' in line or 'synonym' in line: - if 'synonym' in line and 'Bacteria' == node[1]: - #apparently there exists as class of walking sticks called Bacteria Latreilla (629395), that have as synonym name Bacteria - print('wrong Bacteria') + node = [field.strip() for field in line.split("|")] + if "scientific" in line or "synonym" in line: + if "synonym" in line and "Bacteria" == node[1]: + # apparently there exists as class of walking sticks called Bacteria Latreilla (629395), that have as synonym name Bacteria + print("wrong Bacteria") else: if node[1] in tax_names: - orig=tax_names[node[1]] + orig = tax_names[node[1]] if not node[1] in multiple_names: - multiple_names[node[1]]=[] + multiple_names[node[1]] = [] multiple_names[node[1]].append(orig) multiple_names[node[1]].append(node[0]) - tax_names[node[1]]=node[0] + tax_names[node[1]] = node[0] tax_names_reverse[node[0]] = node[1] else: tax_names[node[1]] = node[0] tax_names_reverse[node[0]] = node[1] - if 'scientific' in line: + if "scientific" in line: tax_names_sci[node[0]] = node[1] - return tax_names_reverse,tax_names,tax_names_sci,multiple_names + return tax_names_reverse, tax_names, tax_names_sci, multiple_names + -''' +""" def readNodes(nodes_tax_file): input: @@ -76,44 +130,47 @@ def readNodes(nodes_tax_file): tax_nodes[node[1]].append(node[0]) #couple node with parent tax_types[node[0]] = node[2] #couple node with rank return tax_nodes,tax_types -''' +""" -def readNodes(nodes_tax_file): - ''' +def readNodes(nodes_tax_file): + """ input: - nodes.dmp (NCBI Taxonomy) output: - dictionary of form {parent: node} - dictionary of form {node: type} - ''' + """ tax_nodes = {} tax_types = {} - with open(nodes_tax_file, 'r') as nodes_tax: + with open(nodes_tax_file, "r") as nodes_tax: for line in nodes_tax: - node = [field.strip() for field in line.split('|')] #make list of line - tax_nodes[node[0]] = node[1] #couple node with parent - tax_types[node[0]] = node[2] #couple node with rank - return tax_nodes,tax_types + node = [field.strip() for field in line.split("|")] # make list of line + tax_nodes[node[0]] = node[1] # couple node with parent + tax_types[node[0]] = node[2] # couple node with rank + return tax_nodes, tax_types def getTaxChildren(tax_nodes, tax_types, taxid, ranking): - - ''' + """ input: - dictionary of form {parent: node} (readNodes output) - dictionary of form {node: type} (readNodes output) - taxid output: - dictionary of form {tax_id: [descendants]} - ''' + """ tax_descendants = [] children = [] node = str(taxid) - if node not in tax_types: #check if node in nodes.dmp - print('[Warning] Could not find {} in nodes.dmp while parsing taxonomy hierarchy\n'.format(node)) + if node not in tax_types: # check if node in nodes.dmp + print( + "[Warning] Could not find {} in nodes.dmp while parsing taxonomy hierarchy\n".format( + node + ) + ) else: children.append(node) @@ -127,221 +184,271 @@ def getTaxChildren(tax_nodes, tax_types, taxid, ranking): return tax_descendants -def getTaxParent(tax_nodes, tax_types, taxid, ranking): - ''' +def getTaxParent(tax_nodes, tax_types, taxid, ranking): + """ input: - dictionary of form {parent: node} (readNodes output) - dictionary of form {node: type} (readNodes output) - taxid output: - dictionary of form {tax_id: [descendants]} - ''' + """ tax_parents = {} node = str(taxid) - if node not in tax_types: #check if node in nodes.dmp + if node not in tax_types: # check if node in nodes.dmp tax_parents[node] = None - print('[Warning] Could not find {} in nodes.dmp while parsing taxonomy hierarchy\n'.format(node)) + print( + "[Warning] Could not find {} in nodes.dmp while parsing taxonomy hierarchy\n".format( + node + ) + ) else: - parent = tax_nodes[node] #get parent for current node - tax_parents[node] = [parent] #add node to dictionary + parent = tax_nodes[node] # get parent for current node + tax_parents[node] = [parent] # add node to dictionary - while parent != tax_nodes[parent] and tax_types[parent]!= ranking: #stop when parent = parent(parent) (i.e. 1 = 1) - parent = tax_nodes[parent] #get parent of parent - tax_parents[node].append(parent) #add parent to node in dictionary + while ( + parent != tax_nodes[parent] and tax_types[parent] != ranking + ): # stop when parent = parent(parent) (i.e. 1 = 1) + parent = tax_nodes[parent] # get parent of parent + tax_parents[node].append(parent) # add parent to node in dictionary return tax_parents -api_instance = ncbi.datasets.GenomeApi(ncbi.datasets.ApiClient()) + +api_instance = NcbiApi(args.key) # determine the lineage where your tax id belongs to (lineage taken until upper level = args.type) -taxparents,taxtypes=readNodes(args.nodesfile) -taxnames,namestax,taxnames_sci,multiple_names=readNames(args.namesfile) +taxparents, taxtypes = readNodes(args.nodesfile) +taxnames, namestax, taxnames_sci, multiple_names = readNames(args.namesfile) -eukgens=[] -prokgens=[] +eukgens = [] +prokgens = [] -spoifamily="" -spoiclade="" -spoigenus=args.spoi.split()[0] -spoispecies=args.spoi +spoifamily = "" +spoiclade = "" +spoigenus = args.spoi.split()[0] +spoispecies = args.spoi if spoispecies in namestax: - lineage=getTaxParent(taxparents,taxtypes,namestax[spoispecies],args.type) - lineage2=getTaxParent(taxparents,taxtypes,namestax[spoispecies],'order') + lineage = getTaxParent(taxparents, taxtypes, namestax[spoispecies], args.type) + lineage2 = getTaxParent(taxparents, taxtypes, namestax[spoispecies], "order") if lineage[namestax[spoispecies]] != None: - spoifamily=taxnames[lineage[namestax[spoispecies]][-1]] - spoiclade=taxnames[lineage2[namestax[spoispecies]][-1]] + spoifamily = taxnames[lineage[namestax[spoispecies]][-1]] + spoiclade = taxnames[lineage2[namestax[spoispecies]][-1]] -print(spoigenus+'\t'+spoifamily+'\t'+spoiclade) +print(spoigenus + "\t" + spoifamily + "\t" + spoiclade) -k=open(args.tax,'r') +k = open(args.tax, "r") for line in k: - line=line.strip() - sciname=line.split(';')[-2] - if 'environmental' in sciname: - sciname=line.split(';')[-3] - if sciname == 'uncultured': - sciname=line.split(';')[-3] - if sciname == 'endosymbionts': - sciname=line.split(';')[-3] - if 'Hafnia-Obesumbacterium' in sciname: - sciname=sciname.split('-')[0] - if 'Escherichia-Shigella' in sciname: - sciname=sciname.split('-')[0] - if sciname not in namestax and sciname != 'Unclassified' and sciname != 'Chloroplast' and sciname != 'Mitochondrion': - sciname=line.split(';')[-3] - print('NOW: '+sciname) + line = line.strip() + sciname = line.split(";")[-2] + if "environmental" in sciname: + sciname = line.split(";")[-3] + if sciname == "uncultured": + sciname = line.split(";")[-3] + if sciname == "endosymbionts": + sciname = line.split(";")[-3] + if "Hafnia-Obesumbacterium" in sciname: + sciname = sciname.split("-")[0] + if "Escherichia-Shigella" in sciname: + sciname = sciname.split("-")[0] + if ( + sciname not in namestax + and sciname != "Unclassified" + and sciname != "Chloroplast" + and sciname != "Mitochondrion" + ): + sciname = line.split(";")[-3] + print("NOW: " + sciname) print(sciname) if sciname in namestax: - #print(sciname) - taxid_line=0 - if sciname in multiple_names : + # print(sciname) + taxid_line = 0 + if sciname in multiple_names: print(sciname) - found_true_lineage=False - besttaxid="" - bestcounter=0 + found_true_lineage = False + besttaxid = "" + bestcounter = 0 for elem in multiple_names[sciname]: - lineage=getTaxParent(taxparents,taxtypes,elem,args.type) - fulllineage=getTaxParent(taxparents,taxtypes,elem,'superkingdom') - counterhere=0 + lineage = getTaxParent(taxparents, taxtypes, elem, args.type) + fulllineage = getTaxParent(taxparents, taxtypes, elem, "superkingdom") + counterhere = 0 for x in fulllineage[elem]: - #print(x) - cnt=line.split(';').count(taxnames[x]) + # print(x) + cnt = line.split(";").count(taxnames[x]) if cnt > 0: print(taxnames[x]) - counterhere=counterhere+1 + counterhere = counterhere + 1 if int(counterhere) > int(bestcounter): - bestcounter=counterhere - besttaxid=elem - taxid_line=besttaxid - lineage=getTaxParent(taxparents,taxtypes,taxid_line,args.type) - fulllineage=getTaxParent(taxparents,taxtypes,taxid_line,'superkingdom') - cladelineage=getTaxParent(taxparents,taxtypes,taxid_line,'order') - rootlevelname=taxnames[fulllineage[taxid_line][-1]] - cladelevelname=taxnames[cladelineage[taxid_line][-1]] + bestcounter = counterhere + besttaxid = elem + taxid_line = besttaxid + lineage = getTaxParent(taxparents, taxtypes, taxid_line, args.type) + fulllineage = getTaxParent(taxparents, taxtypes, taxid_line, "superkingdom") + cladelineage = getTaxParent(taxparents, taxtypes, taxid_line, "order") + rootlevelname = taxnames[fulllineage[taxid_line][-1]] + cladelevelname = taxnames[cladelineage[taxid_line][-1]] else: - taxid_line=namestax[sciname] - lineage=getTaxParent(taxparents,taxtypes,namestax[sciname],args.type) - fulllineage=getTaxParent(taxparents,taxtypes,namestax[sciname],'superkingdom') - cladelineage=getTaxParent(taxparents,taxtypes,namestax[sciname],'order') - rootlevelname=taxnames[fulllineage[namestax[sciname]][-1]] - cladelevelname=taxnames[cladelineage[namestax[sciname]][-1]] + taxid_line = namestax[sciname] + lineage = getTaxParent(taxparents, taxtypes, namestax[sciname], args.type) + fulllineage = getTaxParent( + taxparents, taxtypes, namestax[sciname], "superkingdom" + ) + cladelineage = getTaxParent( + taxparents, taxtypes, namestax[sciname], "order" + ) + rootlevelname = taxnames[fulllineage[namestax[sciname]][-1]] + cladelevelname = taxnames[cladelineage[namestax[sciname]][-1]] if taxtypes[taxid_line] == args.type: - print('FAMILY:'+sciname+' CLADE:'+cladelevelname) - taxlevelname=sciname - if 'Eukaryota' == rootlevelname: - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(taxlevelname),page_size=5000) - #cmd=str(args.datasets)+" assembly-descriptors tax-name '"+str(taxlevelname)+"' > "+str(args.outdir)+"/log."+str(taxlevelname)+".json" + print("FAMILY:" + sciname + " CLADE:" + cladelevelname) + taxlevelname = sciname + if "Eukaryota" == rootlevelname: + n_genomes = api_instance.assembly_count_for_taxon( + taxon=str(taxlevelname) + ) else: - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(taxlevelname),filters_assembly_source='refseq',page_size=5000) - #cmd=str(args.datasets)+" assembly-descriptors --refseq tax-name '"+str(taxlevelname)+"' > "+str(args.outdir)+"/log."+str(taxlevelname)+".json" - #os.system(cmd) + n_genomes = api_instance.assembly_count_for_taxon( + taxon=str(taxlevelname), + filters_assembly_source="refseq", + ) foundlevel = False - i=0 - if genome_summary.assemblies is not None: - for assembly in map(lambda d: d.assembly, genome_summary.assemblies): - i=i+1 - if i > 0: + if n_genomes > 0: foundlevel = True - #f.close() - if foundlevel == True: - if 'Eukaryota' == rootlevelname and taxlevelname not in eukgens and taxlevelname != spoifamily and cladelevelname != spoiclade: - fulllineage_euk=taxlevelname - for elem in getTaxParent(taxparents,taxtypes,namestax[taxlevelname],'superkingdom')[namestax[taxlevelname]]: - fulllineage_euk=fulllineage_euk+','+taxnames[elem] + if foundlevel: + if ( + "Eukaryota" == rootlevelname + and taxlevelname not in eukgens + and taxlevelname != spoifamily + and cladelevelname != spoiclade + ): + fulllineage_euk = taxlevelname + for elem in getTaxParent( + taxparents, taxtypes, namestax[taxlevelname], "superkingdom" + )[namestax[taxlevelname]]: + fulllineage_euk = fulllineage_euk + "," + taxnames[elem] print(fulllineage_euk) eukgens.append(fulllineage_euk) - elif taxlevelname not in prokgens and taxlevelname != spoifamily and cladelevelname != spoiclade: + elif ( + taxlevelname not in prokgens + and taxlevelname != spoifamily + and cladelevelname != spoiclade + ): prokgens.append(taxlevelname) else: - print('No genomes in databases') - if cladelevelname != 'root': - taxlevelname=cladelevelname - if 'Eukaryota' == rootlevelname and taxlevelname not in eukgens and taxlevelname != spoifamily and cladelevelname != spoiclade: - fulllineage_euk=taxlevelname - for elem in getTaxParent(taxparents,taxtypes,namestax[taxlevelname],'superkingdom')[namestax[taxlevelname]]: - fulllineage_euk=fulllineage_euk+','+taxnames[elem] + print("No genomes in databases") + if cladelevelname != "root": + taxlevelname = cladelevelname + if ( + "Eukaryota" == rootlevelname + and taxlevelname not in eukgens + and taxlevelname != spoifamily + and cladelevelname != spoiclade + ): + fulllineage_euk = taxlevelname + for elem in getTaxParent( + taxparents, taxtypes, namestax[taxlevelname], "superkingdom" + )[namestax[taxlevelname]]: + fulllineage_euk = fulllineage_euk + "," + taxnames[elem] print(fulllineage_euk) - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(taxlevelname.split(',')[0]),page_size=5000) - i=0 - if genome_summary.assemblies is not None: - for assembly in map(lambda d: d.assembly, genome_summary.assemblies): - i=i+1 - if i > 0: - eukgens.append(fulllineage_euk) - elif taxlevelname not in prokgens and taxlevelname != spoifamily and cladelevelname != spoiclade: - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(taxlevelname),page_size=5000) - i=0 - if genome_summary.assemblies is not None: - for assembly in map(lambda d: d.assembly, genome_summary.assemblies): - i=i+1 - if i > 0: - prokgens.append(taxlevelname) + n_genomes = api_instance.assembly_count_for_taxon( + taxon=str(taxlevelname.split(",")[0]) + ) + if n_genomes > 0: + eukgens.append(fulllineage_euk) + elif ( + taxlevelname not in prokgens + and taxlevelname != spoifamily + and cladelevelname != spoiclade + ): + n_genomes = api_instance.assembly_count_for_taxon( + taxon=str(taxlevelname) + ) + if n_genomes > 0: + prokgens.append(taxlevelname) elif lineage[taxid_line] != None: - print('DIFFERENT THAN FAMILY:'+sciname+ ' CLADE:'+cladelevelname) - taxlevelname=taxnames_sci[lineage[taxid_line][-1]] + print("DIFFERENT THAN FAMILY:" + sciname + " CLADE:" + cladelevelname) + taxlevelname = taxnames_sci[lineage[taxid_line][-1]] if int(lineage[taxid_line][-1]) != 1: - if 'Eukaryota' == rootlevelname: - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(taxlevelname),page_size=5000) - #cmd=str(args.datasets)+" assembly-descriptors tax-name '"+str(taxlevelname)+"' > "+str(args.outdir)+"/log."+str(taxlevelname)+".json" + if "Eukaryota" == rootlevelname: + n_genomes = api_instance.assembly_count_for_taxon( + taxon=str(taxlevelname) + ) else: - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(taxlevelname),filters_assembly_source='refseq',page_size=5000) - #cmd=str(args.datasets)+" assembly-descriptors --refseq tax-name '"+str(taxlevelname)+"' > "+str(args.outdir)+"/log."+str(taxlevelname)+".json" - #os.system(cmd) + n_genomes = api_instance.assembly_count_for_taxon( + taxon=str(taxlevelname), + filters_assembly_source="refseq", + ) foundlevel = False - i=0 - if genome_summary.assemblies is not None: - for assembly in map(lambda d: d.assembly, genome_summary.assemblies): - i=i+1 - if i > 0: + if n_genomes > 0: foundlevel = True - #f.close() - if foundlevel == True: - #if taxtypes[namestax[sciname]] == args.type: - if 'Eukaryota' == rootlevelname and taxlevelname not in eukgens and taxlevelname != spoifamily and cladelevelname != spoiclade: - fulllineage_euk=taxlevelname - for elem in getTaxParent(taxparents,taxtypes,namestax[taxlevelname],'superkingdom')[namestax[taxlevelname]]: - fulllineage_euk=fulllineage_euk+','+taxnames[elem] + if foundlevel == True: + # if taxtypes[namestax[sciname]] == args.type: + if ( + "Eukaryota" == rootlevelname + and taxlevelname not in eukgens + and taxlevelname != spoifamily + and cladelevelname != spoiclade + ): + fulllineage_euk = taxlevelname + for elem in getTaxParent( + taxparents, taxtypes, namestax[taxlevelname], "superkingdom" + )[namestax[taxlevelname]]: + fulllineage_euk = fulllineage_euk + "," + taxnames[elem] print(fulllineage_euk) eukgens.append(fulllineage_euk) - elif taxlevelname not in prokgens and taxlevelname != spoifamily and cladelevelname != spoiclade: + elif ( + taxlevelname not in prokgens + and taxlevelname != spoifamily + and cladelevelname != spoiclade + ): prokgens.append(taxlevelname) else: - print('No genomes in databases') - if cladelevelname != 'root': - taxlevelname=cladelevelname - if 'Eukaryota' == rootlevelname and taxlevelname not in eukgens and taxlevelname != spoifamily and cladelevelname != spoiclade: - fulllineage_euk=taxlevelname - for elem in getTaxParent(taxparents,taxtypes,namestax[taxlevelname],'superkingdom')[namestax[taxlevelname]]: - fulllineage_euk=fulllineage_euk+','+taxnames[elem] + print("No genomes in databases") + if cladelevelname != "root": + taxlevelname = cladelevelname + if ( + "Eukaryota" == rootlevelname + and taxlevelname not in eukgens + and taxlevelname != spoifamily + and cladelevelname != spoiclade + ): + fulllineage_euk = taxlevelname + for elem in getTaxParent( + taxparents, + taxtypes, + namestax[taxlevelname], + "superkingdom", + )[namestax[taxlevelname]]: + fulllineage_euk = fulllineage_euk + "," + taxnames[elem] print(fulllineage_euk) - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(taxlevelname.split(',')[0]),page_size=5000) - i=0 - if genome_summary.assemblies is not None: - for assembly in map(lambda d: d.assembly, genome_summary.assemblies): - i=i+1 - if i > 0: - eukgens.append(fulllineage_euk) - elif taxlevelname not in prokgens and taxlevelname != spoifamily and cladelevelname != spoiclade: - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(taxlevelname),page_size=5000) - i=0 - if genome_summary.assemblies is not None: - for assembly in map(lambda d: d.assembly, genome_summary.assemblies): - i=i+1 - if i > 0: - prokgens.append(taxlevelname) + n_genomes = api_instance.assembly_count_for_taxon( + taxon=str(taxlevelname.split(",")[0]) + ) + if n_genomes > 0: + eukgens.append(fulllineage_euk) + elif ( + taxlevelname not in prokgens + and taxlevelname != spoifamily + and cladelevelname != spoiclade + ): + genome_summary = api_instance.assembly_count_for_taxon( + taxon=str(taxlevelname) + ) + i = 0 + if n_genomes > 0: + prokgens.append(taxlevelname) else: - print('NOT FOUND:'+sciname+'\t'+line) + print("NOT FOUND:" + sciname + "\t" + line) -file1=args.outdir+"/prok."+args.suffix -k=open(file1,'w') +file1 = args.outdir + "/prok." + args.suffix +k = open(file1, "w") for elem in prokgens: - k.write(elem+'\n') + k.write(elem + "\n") k.close() -file1=args.outdir+"/euk."+args.suffix -k=open(file1,'w') +file1 = args.outdir + "/euk." + args.suffix +k = open(file1, "w") for elem in eukgens: - k.write(elem+'\n') + k.write(elem + "\n") k.close() diff --git a/scripts/FetchGenomesRefSeq.py b/scripts/FetchGenomesRefSeq.py index 5ce92c6..30f60a1 100644 --- a/scripts/FetchGenomesRefSeq.py +++ b/scripts/FetchGenomesRefSeq.py @@ -1,116 +1,184 @@ from __future__ import division -import json import time import argparse import os -try: - import ncbi.datasets -except ImportError: - print('ncbi.datasets module not found. To install, run `pip install ncbi-datasets-pylib`.') -from ncbi.datasets.package import dataset +from pathlib import Path + +from NCBIApiTools import NcbiApi parser = argparse.ArgumentParser() -parser.add_argument('--taxname', action="store", dest="tax", type=str, help='a genus taxname to download refseq genomes for') -parser.add_argument('--dir', action="store", dest="dir", type=str, help='base directory') -parser.add_argument('--refseq', action="store", dest="refs", type=str, help='all or refseq database') +parser.add_argument( + "--taxname", + action="store", + dest="tax", + type=str, + help="a genus taxname to download refseq genomes for", +) +parser.add_argument( + "--dir", action="store", dest="dir", type=str, help="base directory" +) +parser.add_argument( + "--refseq", action="store", dest="refs", type=str, help="all or refseq database" +) +parser.add_argument( + "-k", + type=str, + action="store", + dest="key", + metavar="NCBI API KEY", + help="NCBI API Key", + required=False, + default=os.environ.get("NCBI_API_KEY"), +) args = parser.parse_args() -def Average(lst): - return sum(lst) / len(lst) -api_instance = ncbi.datasets.GenomeApi(ncbi.datasets.ApiClient()) +def Average(lst): + return sum(lst) / len(lst) + + +api_instance = NcbiApi(args.key) -taxname=str(args.tax).split("genus.")[1].split('.')[0] -taxname_orig=taxname -if '_' in taxname: - taxname=taxname.replace('_',' ') +taxname = str(args.tax).split("genus.")[1].split(".")[0] +taxname_orig = taxname +if "_" in taxname: + taxname = taxname.replace("_", " ") if not os.path.exists(args.dir): os.makedirs(args.dir) # fetch data from NCBI via 'datasets' of all species from that clade -if args.refs == 'yes': - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(taxname),filters_reference_only=True,page_size=5000) -else: - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(taxname),page_size=5000) - -SpeciesDictionary={} -i=0 -for assembly in map(lambda d: d.assembly, genome_summary.assemblies): - i=i+1 - contiguity=int(assembly.contig_n50) - acc=assembly.assembly_accession - sciname_orig=assembly.org.sci_name - strainname=assembly.org.strain - if strainname != None: - if strainname in sciname_orig and 'sp' not in sciname_orig: - sciname=sciname_orig.replace(strainname,'').strip() +is_refseq = args.refs.lower() == "yes" if args.refs else False + +assemblies = api_instance.get_assemblies_for_taxon( + taxon=str(taxname), filters_reference_only=is_refseq, page_size=1000 +) + +SpeciesDictionary = {} +for assembly in assemblies: + acc = assembly.get("accession") + date = assembly.get("release_date") + sciname_orig = assembly.get("organism").get("organism_name") + strainname = assembly.get("organism").get("infraspecific_names").get("strain") + contiguity = int(assembly.get("assembly_stats").get("contig_n50")) + size = int(assembly.get("assembly_stats").get("total_sequence_length")) + + if strainname: + if strainname in sciname_orig and "sp" not in sciname_orig: + sciname = sciname_orig.replace(strainname, "").strip() else: - sciname=sciname_orig + sciname = sciname_orig else: - sciname=sciname_orig + sciname = sciname_orig + print(sciname) + if sciname not in SpeciesDictionary: - SpeciesDictionary[sciname]={} - SpeciesDictionary[sciname]['ReleaseDate']=assembly.submission_date - SpeciesDictionary[sciname]['Identifier']=acc - SpeciesDictionary[sciname]['GenomeSize']=int(assembly.seq_length) - SpeciesDictionary[sciname]['N50']=contiguity + SpeciesDictionary[sciname] = {} + SpeciesDictionary[sciname]["ReleaseDate"] = date + SpeciesDictionary[sciname]["Identifier"] = acc + SpeciesDictionary[sciname]["GenomeSize"] = size + SpeciesDictionary[sciname]["N50"] = contiguity else: - novel_submission = time.strptime(assembly.submission_date, "%Y-%m-%d") - old_submission = time.strptime(SpeciesDictionary[sciname]['ReleaseDate'], "%Y-%m-%d") - if novel_submission > old_submission or contiguity > SpeciesDictionary[sciname]['N50']: - SpeciesDictionary[sciname]['ReleaseDate']=assembly.submission_date - SpeciesDictionary[sciname]['Identifier']=acc - SpeciesDictionary[sciname]['GenomeSize']=int(assembly.seq_length) - SpeciesDictionary[sciname]['N50']=contiguity - - -accs=[] -genomesizes=[] -j=0 + novel_submission = time.strptime(date, "%Y-%m-%d") + old_submission = time.strptime( + SpeciesDictionary[sciname]["ReleaseDate"], "%Y-%m-%d" + ) + if ( + novel_submission > old_submission + or contiguity > SpeciesDictionary[sciname]["N50"] + ): + SpeciesDictionary[sciname]["ReleaseDate"] = date + SpeciesDictionary[sciname]["Identifier"] = acc + SpeciesDictionary[sciname]["GenomeSize"] = size + SpeciesDictionary[sciname]["N50"] = contiguity + + +accs = [] +genomesizes = [] for species in SpeciesDictionary: - print(accs) - j=j+1 - #if not os.path.exists(args.dir2+"/"+SpeciesDictionary[species]['Identifier']): - accs.append(SpeciesDictionary[species]['Identifier']) - genomesizes.append(SpeciesDictionary[species]['GenomeSize']) + accs.append(SpeciesDictionary[species]["Identifier"]) + genomesizes.append(SpeciesDictionary[species]["GenomeSize"]) if len(accs) > 0: - nrgfs=int(round(float(Average(genomesizes)/1000000)*20)) - print('Genomes for '+str(j)+' species') - print("Number of necessary GFs "+str(nrgfs)) - print(f'Download a package for {accs}.') - print('Begin download of genome data package ...') - print(str(len(accs))+' genomes') - zipfile_name = str(args.dir)+"/"+"/RefSeq."+str(taxname_orig)+".zip" + nrgfs = int(round(float(Average(genomesizes) / 1000000) * 20)) + print("Genomes for " + str(len(accs)) + " species") + print("Number of necessary GFs " + str(nrgfs)) + print(f"Download a package for {accs}.") + print("Begin download of genome data package ...") + print(str(len(accs)) + " genomes") + zipfile_name = str(args.dir) + "/" + "/RefSeq." + str(taxname_orig) + ".zip" for t in range(0, len(accs), 100): - accshort=accs[t:t + 100] - try: - zipfile_name_part = str(args.dir)+"/"+"/RefSeq."+str(taxname_orig)+".part"+str(t)+".zip" - api_response = api_instance.download_assembly_package(accshort,exclude_sequence=False, hydrated='FULLY_HYDRATED',_preload_content=False,filename=zipfile_name_part) - with open(zipfile_name_part, 'wb') as f: - f.write(api_response.data) - print('Download complete part '+str(t)) - cmd="unzip -d "+str(args.dir)+"/"+str(taxname_orig)+".RefSeq.part"+str(t)+" "+str(args.dir)+"/"+"/RefSeq."+str(taxname_orig)+".part"+str(t)+".zip" - os.system(cmd) - except ncbi.datasets.openapi.ApiException as e: - print("Exception when calling GenomeApi->download_assembly_package: %s\n" % e) - cmd="mkdir "+str(args.dir)+"/"+str(taxname_orig)+".Refseq" + accshort = accs[t : t + 100] + zipfile_name_part = ( + Path(args.dir) / f"/RefSeq.{str(taxname_orig)}.part{str(t)}.zip" + ) + api_response = api_instance.download_genomes( + accshort, + outfile=zipfile_name_part, + ) + print("Download complete part " + str(t)) + cmd = "unzip -d " + str(zipfile_name_part)[-4] + " " + str(zipfile_name_part) + os.system(cmd) + + cmd = "mkdir " + str(args.dir) + "/" + str(taxname_orig) + ".Refseq" os.system(cmd) - cmd="mkdir "+str(args.dir)+"/"+str(taxname_orig)+".Refseq/ncbi_dataset" + cmd = "mkdir " + str(args.dir) + "/" + str(taxname_orig) + ".Refseq/ncbi_dataset" os.system(cmd) - cmd="mkdir "+str(args.dir)+"/"+str(taxname_orig)+".Refseq/ncbi_dataset/data" + cmd = ( + "mkdir " + str(args.dir) + "/" + str(taxname_orig) + ".Refseq/ncbi_dataset/data" + ) os.system(cmd) - cmd="cp -r "+str(args.dir)+"/"+str(taxname_orig)+".RefSeq.part*/ncbi_dataset/data/G* "+str(args.dir)+"/"+str(taxname_orig)+".Refseq/ncbi_dataset/data/" + cmd = ( + "cp -r " + + str(args.dir) + + "/" + + str(taxname_orig) + + ".RefSeq.part*/ncbi_dataset/data/G* " + + str(args.dir) + + "/" + + str(taxname_orig) + + ".Refseq/ncbi_dataset/data/" + ) os.system(cmd) - cmd="cat "+str(args.dir)+"/"+str(taxname_orig)+".RefSeq.part*/ncbi_dataset/data/assembly_data_report.jsonl >> "+str(args.dir)+"/"+str(taxname_orig)+".Refseq/ncbi_dataset/data/assembly_data_report.jsonl" + cmd = ( + "cat " + + str(args.dir) + + "/" + + str(taxname_orig) + + ".RefSeq.part*/ncbi_dataset/data/assembly_data_report.jsonl >> " + + str(args.dir) + + "/" + + str(taxname_orig) + + ".Refseq/ncbi_dataset/data/assembly_data_report.jsonl" + ) os.system(cmd) - cmd="cat "+str(args.dir)+"/"+str(taxname_orig)+".RefSeq.part*/ncbi_dataset/data/dataset_catalog.json >> "+str(args.dir)+"/"+str(taxname_orig)+".Refseq/ncbi_dataset/data/dataset_catalog.json" + cmd = ( + "cat " + + str(args.dir) + + "/" + + str(taxname_orig) + + ".RefSeq.part*/ncbi_dataset/data/dataset_catalog.json >> " + + str(args.dir) + + "/" + + str(taxname_orig) + + ".Refseq/ncbi_dataset/data/dataset_catalog.json" + ) os.system(cmd) - cmd="rm -r "+str(args.dir)+"/"+"/RefSeq."+str(taxname_orig)+".part*.zip "+str(args.dir)+"/"+str(taxname_orig)+".RefSeq.part*" + cmd = ( + "rm -r " + + str(args.dir) + + "/" + + "/RefSeq." + + str(taxname_orig) + + ".part*.zip " + + str(args.dir) + + "/" + + str(taxname_orig) + + ".RefSeq.part*" + ) os.system(cmd) else: - print('No genomes available') - cmd="touch "+str(args.dir)+"/"+str(taxname_orig)+".download.log" + print("No genomes available") + cmd = "touch " + str(args.dir) + "/" + str(taxname_orig) + ".download.log" os.system(cmd) diff --git a/scripts/FetchGenomesRefSeqRelatives.py b/scripts/FetchGenomesRefSeqRelatives.py index 68acad2..e4e04e3 100644 --- a/scripts/FetchGenomesRefSeqRelatives.py +++ b/scripts/FetchGenomesRefSeqRelatives.py @@ -1,186 +1,261 @@ from __future__ import division -import json -import time + import argparse import os -try: - import ncbi.datasets -except ImportError: - print('ncbi.datasets module not found. To install, run `pip install ncbi-datasets-pylib`.') -from ncbi.datasets.package import dataset +import time +from pathlib import Path + +from NCBIApiTools import NcbiApi parser = argparse.ArgumentParser() -parser.add_argument('--taxname', action="store", dest="tax", type=str, help='scientific name of species of interest') -parser.add_argument('--dir', action="store", dest="dir", type=str, help='base directory') -parser.add_argument("-na", type=str, action='store', dest='namesfile', metavar='NAMES',help='NCBI names.dmp') -parser.add_argument("-no", type=str, action='store', dest='nodesfile', metavar='NODES',help='NCBI nodes.dmp') -parser.add_argument('--refseq', action="store", dest="refs", type=str, help='all or refseq database') +parser.add_argument( + "--taxname", + action="store", + dest="tax", + type=str, + help="scientific name of species of interest", +) +parser.add_argument( + "--dir", action="store", dest="dir", type=str, help="base directory" +) +parser.add_argument( + "-na", + type=str, + action="store", + dest="namesfile", + metavar="NAMES", + help="NCBI names.dmp", +) +parser.add_argument( + "-no", + type=str, + action="store", + dest="nodesfile", + metavar="NODES", + help="NCBI nodes.dmp", +) +parser.add_argument( + "--refseq", action="store", dest="refs", type=str, help="all or refseq database" +) +parser.add_argument( + "-k", + type=str, + action="store", + dest="key", + metavar="NCBI API KEY", + help="NCBI API Key", + required=False, + default=os.environ.get("NCBI_API_KEY"), +) args = parser.parse_args() -def Average(lst): - return sum(lst) / len(lst) + +def Average(lst): + return sum(lst) / len(lst) + def readNames(names_tax_file): - ''' + """ input: - name.dmp (NCBI Taxonomy) output: - dictionary of form {node: name} - dictionary of form {sci name: node} - ''' + """ tax_names = {} - tax_names_reverse= {} - with open(names_tax_file, 'r') as nodes_tax: + tax_names_reverse = {} + with open(names_tax_file, "r") as nodes_tax: for line in nodes_tax: - node = [field.strip() for field in line.split('|')] - if 'scientific' in line or 'synonym' in line: + node = [field.strip() for field in line.split("|")] + if "scientific" in line or "synonym" in line: tax_names[node[1]] = node[0] tax_names_reverse[node[0]] = node[1] - return tax_names_reverse,tax_names + return tax_names_reverse, tax_names -def readNodes(nodes_tax_file): - ''' +def readNodes(nodes_tax_file): + """ input: - nodes.dmp (NCBI Taxonomy) output: - dictionary of form {parent: node} - dictionary of form {node: type} - ''' + """ tax_nodes = {} tax_types = {} - with open(nodes_tax_file, 'r') as nodes_tax: + with open(nodes_tax_file, "r") as nodes_tax: for line in nodes_tax: - node = [field.strip() for field in line.split('|')] #make list of line - tax_nodes[node[0]] = node[1] #couple node with parent - tax_types[node[0]] = node[2] #couple node with rank + node = [field.strip() for field in line.split("|")] # make list of line + tax_nodes[node[0]] = node[1] # couple node with parent + tax_types[node[0]] = node[2] # couple node with rank return tax_nodes -api_instance = ncbi.datasets.GenomeApi(ncbi.datasets.ApiClient()) +api_instance = NcbiApi(args.key) if not os.path.exists(args.dir): os.makedirs(args.dir) -taxparents=readNodes(args.nodesfile) -taxnames,namestax=readNames(args.namesfile) +taxparents = readNodes(args.nodesfile) +taxnames, namestax = readNames(args.namesfile) + +is_refseq = args.refs.lower() == "yes" if args.refs else False if args.tax in namestax: taxid = namestax[args.tax] parent = taxparents[taxid] parentname = taxnames[parent] foundlevel = False - print(str(taxid)+'\t'+str(parent)+'\t'+parentname) + print(str(taxid) + "\t" + str(parent) + "\t" + parentname) while parent != taxparents[parent]: time.sleep(1) parentname = taxnames[parent] parentname_combi = parentname - if ' ' in parentname: - parentname_combi = parentname.replace(' ','_') - print(str(parent)+'\t'+parentname+'\t'+parentname_combi+'\t'+args.tax) + if " " in parentname: + parentname_combi = parentname.replace(" ", "_") + print( + str(parent) + "\t" + parentname + "\t" + parentname_combi + "\t" + args.tax + ) # fetch data from NCBI via 'datasets' of all species from that clade - if args.refs == 'yes': - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(parent),filters_reference_only=True,page_size=5000) - else: - genome_summary = api_instance.assembly_descriptors_by_taxon(taxon=str(parent),page_size=5000) - i=0 - print(genome_summary.assemblies) - if genome_summary.assemblies is not None: - for assembly in map(lambda d: d.assembly, genome_summary.assemblies): - #print(assembly.org.sci_name) - sciname_orig=assembly.org.sci_name - strainname=assembly.org.strain - #print(sciname_orig) - if strainname != None: - if strainname in sciname_orig and 'sp' not in sciname_orig: - sciname=sciname_orig.replace(strainname,'').strip() + assemblies = api_instance.get_assemblies_for_taxon( + taxon=str(parent), filters_reference_only=is_refseq, page_size=1000 + ) + + print(assemblies) + i = 0 + if len(assemblies) > 0: + for assembly in assemblies: + # print(assembly.org.sci_name) + sciname_orig = assembly.get("organism").get("organism_name") + strainname = ( + assembly.get("organism").get("infraspecific_names").get("strain") + ) + if strainname: + if strainname in sciname_orig and "sp" not in sciname_orig: + sciname = sciname_orig.replace(strainname, "").strip() else: - sciname=sciname_orig + sciname = sciname_orig else: - sciname=sciname_orig - #print(sciname) + sciname = sciname_orig if sciname != args.tax: - i=i+1 + i = i + 1 if i > 0: break parent = taxparents[parent] - - SpeciesDictionary={} - i=0 - for assembly in map(lambda d: d.assembly, genome_summary.assemblies): - contiguity=int(assembly.contig_n50) - i=i+1 - acc=assembly.assembly_accession - sciname_orig=assembly.org.sci_name - strainname=assembly.org.strain - if strainname != None: - if strainname in sciname_orig and 'sp' not in sciname_orig: - sciname=sciname_orig.replace(strainname,'').strip() + + SpeciesDictionary = {} + for assembly in assemblies: + acc = assembly.get("accession") + date = assembly.get("release_date") + sciname_orig = assembly.get("organism").get("organism_name") + strainname = assembly.get("organism").get("infraspecific_names").get("strain") + contiguity = int(assembly.get("assembly_stats").get("contig_n50")) + size = int(assembly.get("assembly_stats").get("total_sequence_length")) + + if strainname: + if strainname in sciname_orig and "sp" not in sciname_orig: + sciname = sciname_orig.replace(strainname, "").strip() else: - sciname=sciname_orig + sciname = sciname_orig else: - sciname=sciname_orig + sciname = sciname_orig + if sciname != args.tax: - print('SECOND STEP') - #print(assembly.org.sci_name) - #print(sciname_orig) - print(sciname) + print("SECOND STEP") + print(sciname) if sciname not in SpeciesDictionary: - #print(sciname) - SpeciesDictionary[sciname]={} - SpeciesDictionary[sciname]['ReleaseDate']=assembly.submission_date - SpeciesDictionary[sciname]['Identifier']=acc - SpeciesDictionary[sciname]['GenomeSize']=int(assembly.seq_length) - SpeciesDictionary[sciname]['N50']=contiguity + SpeciesDictionary[sciname] = {} + SpeciesDictionary[sciname]["ReleaseDate"] = date + SpeciesDictionary[sciname]["Identifier"] = acc + SpeciesDictionary[sciname]["GenomeSize"] = size + SpeciesDictionary[sciname]["N50"] = contiguity else: - novel_submission = time.strptime(assembly.submission_date, "%Y-%m-%d") - old_submission = time.strptime(SpeciesDictionary[sciname]['ReleaseDate'], "%Y-%m-%d") - if novel_submission > old_submission or contiguity > SpeciesDictionary[sciname]['N50']: - SpeciesDictionary[sciname]['ReleaseDate']=assembly.submission_date - SpeciesDictionary[sciname]['Identifier']=acc - SpeciesDictionary[sciname]['GenomeSize']=int(assembly.seq_length) - SpeciesDictionary[sciname]['N50']=contiguity - - accs=[] + novel_submission = time.strptime(date, "%Y-%m-%d") + old_submission = time.strptime( + SpeciesDictionary[sciname]["ReleaseDate"], "%Y-%m-%d" + ) + if ( + novel_submission > old_submission + or contiguity > SpeciesDictionary[sciname]["N50"] + ): + SpeciesDictionary[sciname]["ReleaseDate"] = date + SpeciesDictionary[sciname]["Identifier"] = acc + SpeciesDictionary[sciname]["GenomeSize"] = size + SpeciesDictionary[sciname]["N50"] = contiguity + + accs = [] for species in SpeciesDictionary: - print(accs) - #if not os.path.exists(args.dir2+"/"+SpeciesDictionary[species]['Identifier']): - accs.append(SpeciesDictionary[species]['Identifier']) - print(f'Download a package for {accs}.') - print('Begin download of genome data package ...') - #zipfile_name = str(args.dir)+"/"+"/RefSeq.relatives.zip" + accs.append(SpeciesDictionary[species]["Identifier"]) + print(f"Download a package for {accs}.") + print("Begin download of genome data package ...") for t in range(0, len(accs), 3): - zipfile_name_part = str(args.dir)+"/"+"/RefSeq.relatives.part"+str(t)+".zip" - accshort=accs[t:t + 3] - try: - api_response = api_instance.download_assembly_package(accshort,exclude_sequence=False, hydrated='FULLY_HYDRATED',_preload_content=False,filename=zipfile_name_part) - with open(zipfile_name_part, 'wb') as f: - f.write(api_response.data) - print('Download complete part '+str(t)) - cmd="unzip -d "+str(args.dir)+"/relatives.RefSeq.part"+str(t)+" "+str(args.dir)+"/"+"/RefSeq.relatives.part"+str(t)+".zip" - os.system(cmd) - except ncbi.datasets.openapi.ApiException as e: - print("Exception when calling GenomeApi->download_assembly_package: %s\n" % e) - cmd="mkdir "+str(args.dir)+"/relatives.Refseq" + zipfile_name_part = Path(str(args.dir)) / f"RefSeq.relatives.part{str(t)}.zip" + accshort = accs[t : t + 3] + api_response = api_instance.download_genomes( + accshort, + outfile=zipfile_name_part, + ) + print("Download complete part " + str(t)) + cmd = "unzip -d " + str(zipfile_name_part)[-4] + " " + str(zipfile_name_part) + os.system(cmd) + + cmd = "mkdir " + str(args.dir) + "/relatives.Refseq" + os.system(cmd) + cmd = "mkdir " + str(args.dir) + "/relatives.Refseq/ncbi_dataset" os.system(cmd) - cmd="mkdir "+str(args.dir)+"/relatives.Refseq/ncbi_dataset" + cmd = "mkdir " + str(args.dir) + "/relatives.Refseq/ncbi_dataset/data" os.system(cmd) - cmd="mkdir "+str(args.dir)+"/relatives.Refseq/ncbi_dataset/data" + cmd = ( + "cp -r " + + str(args.dir) + + "/relatives.RefSeq.part*/ncbi_dataset/data/G* " + + str(args.dir) + + "/relatives.Refseq/ncbi_dataset/data/" + ) os.system(cmd) - cmd="cp -r "+str(args.dir)+"/relatives.RefSeq.part*/ncbi_dataset/data/G* "+str(args.dir)+"/relatives.Refseq/ncbi_dataset/data/" + cmd = ( + "cat " + + str(args.dir) + + "/relatives.RefSeq.part*/ncbi_dataset/data/assembly_data_report.jsonl >> " + + str(args.dir) + + "/relatives.Refseq/ncbi_dataset/data/assembly_data_report.jsonl" + ) os.system(cmd) - cmd="cat "+str(args.dir)+"/relatives.RefSeq.part*/ncbi_dataset/data/assembly_data_report.jsonl >> "+str(args.dir)+"/relatives.Refseq/ncbi_dataset/data/assembly_data_report.jsonl" + cmd = ( + "cat " + + str(args.dir) + + "/relatives.RefSeq.part*/ncbi_dataset/data/dataset_catalog.json >> " + + str(args.dir) + + "/relatives.Refseq/ncbi_dataset/data/dataset_catalog.json" + ) os.system(cmd) - cmd="cat "+str(args.dir)+"/relatives.RefSeq.part*/ncbi_dataset/data/dataset_catalog.json >> "+str(args.dir)+"/relatives.Refseq/ncbi_dataset/data/dataset_catalog.json" + cmd = ( + "rm -r " + + str(args.dir) + + "/" + + "/RefSeq.relatives.part*.zip " + + str(args.dir) + + "/relatives.RefSeq.part*" + ) os.system(cmd) - cmd="rm -r "+str(args.dir)+"/"+"/RefSeq.relatives.part*.zip "+str(args.dir)+"/relatives.RefSeq.part*" - os.system(cmd) - cmd="cat "+str(args.dir)+"/relatives.Refseq/ncbi_dataset/data/*/*fna > "+str(args.dir)+"/relatives.Refseq/acc.fasta" + cmd = ( + "cat " + + str(args.dir) + + "/relatives.Refseq/ncbi_dataset/data/*/*fna > " + + str(args.dir) + + "/relatives.Refseq/acc.fasta" + ) os.system(cmd) - cmd="dustmasker -in "+str(args.dir)+"/relatives.Refseq/acc.fasta -outfmt fasta | sed -e '/^>/!s/[a-z]/x/g' > "+str(args.dir)+"/relatives.Refseq/masked.fna" + cmd = ( + "dustmasker -in " + + str(args.dir) + + "/relatives.Refseq/acc.fasta -outfmt fasta | sed -e '/^>/!s/[a-z]/x/g' > " + + str(args.dir) + + "/relatives.Refseq/masked.fna" + ) os.system(cmd) - cmd="rm "+str(args.dir)+"/relatives.Refseq/acc.fasta " + cmd = "rm " + str(args.dir) + "/relatives.Refseq/acc.fasta " os.system(cmd) else: - print('No taxonomic name was not found') \ No newline at end of file + print("No taxonomic name was not found") diff --git a/scripts/NCBIApiTools.py b/scripts/NCBIApiTools.py new file mode 100644 index 0000000..42a997c --- /dev/null +++ b/scripts/NCBIApiTools.py @@ -0,0 +1,119 @@ +import sys +from pathlib import Path + +import requests + + +class NcbiApi: + def __init__(self, key): + self.ncbi_api_uri = "https://api.ncbi.nlm.nih.gov/datasets/v2" + self.ncbi_api_key = key + + def get_request(self, command: str) -> requests.Response: + """ + Returns a requests response object for a given NCBI api query (command). + + args: + command -> str: request to append to NCBI API URI + """ + if self.ncbi_api_key != "": + headers = {"Accept": "application/json", "api-key": f"{self.ncbi_api_key}"} + else: + headers = {"Accept": "application/json"} + + response = requests.get(self.ncbi_api_uri + command, headers=headers) + + if response.status_code != 200: + raise Exception( + f"Cannot connect to NCBI (status code '{str(response.status_code)}')'" + ) + + return response + + def get_assemblies_for_taxon( + self, taxon: str, filters_reference_only: bool = False, page_size: int = 1000 + ): + """ + For a given taxon name or taxid, a get a list describing all available + assemblies for a taxon. + + args: + taxon -> string: Taxon name or taxid to query + filters_reference_only -> bool: Return reference assemblies only + page_size -> int: Number of assembies to return per page (max 1000) + """ + query = f"/genome/taxon/{taxon}/dataset_report" + arg_filter_refs_only = f"?filters_reference_only={filters_reference_only}" + arg_page_size = f"&page_size={page_size}" + r = self.get_request(query + arg_filter_refs_only + arg_page_size) + + assemblies = [] + + if r.get("reports"): + assemblies = [] + r.get("reports") + + while r.get("next_page_token", "") != "": + page_token_arg = f"&page_token={r.get('next_page_token')}" + r = self.get_request( + query + arg_filter_refs_only + arg_page_size + page_token_arg + ) + if r.get("reports", None): + assemblies = assemblies + r.get("reports") + + return assemblies + + def assembly_count_for_taxon( + self, taxon: str, filters_assembly_source: str = "all" + ) -> int: + """ + For a given taxon name or taxid, return the number of assemblies available + on NCBI. + + args: + taxon -> string: Taxon name or taxid to query + filters_assembly_source -> string: "refseq", "genbank", or "all" + """ + query = f"/genome/taxon/{taxon}/dataset_report" + filter_source = f"?filters_assembly_source={filters_assembly_source}" + page_size = "&page_size=1" + taxon = self.get_request(query + filter_source + page_size) + + return taxon.json().get("total_count", 0) + + def download_genomes(self, accessions: list, outfile: Path) -> Path: + if self.ncbi_api_key != "": + headers = { + "Accept": "application/zip", + "api-key": f"{self.ncbi_api_key}", + "content-type": "application/json", + } + else: + headers = {"Accept": "application/zip", "content-type": "application/json"} + + try: + r = requests.post( + self.ncbi_api_uri + "/genome/download", + headers=headers, + json={ + "accessions": accessions, + "include_annotation_type": ["GENOME_FASTA"], + "hydrated": "FULLY_HYDRATED", + }, + stream=True, + ) + r.raise_for_status() + + downloaded = 0 + with open(outfile, "wb") as f: + for chunk in r.iter_content(chunk_size=8192): + if chunk: + f.write(chunk) + downloaded += len(chunk) + print( + f"\rDownloading: {downloaded / 1024 / 1024:.1f}MB", end="" + ) + + return outfile + + except requests.HTTPError as e: + print(f"Failed to download genomes: {e}") From fbc1ab66da2b6efe6aeec4a9d16ddedb123d9638 Mon Sep 17 00:00:00 2001 From: Jim Downie Date: Wed, 8 Oct 2025 10:28:05 +0100 Subject: [PATCH 02/13] fix: properly exit when download fails --- scripts/NCBIApiTools.py | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) diff --git a/scripts/NCBIApiTools.py b/scripts/NCBIApiTools.py index 42a997c..ea0bb9b 100644 --- a/scripts/NCBIApiTools.py +++ b/scripts/NCBIApiTools.py @@ -116,4 +116,5 @@ def download_genomes(self, accessions: list, outfile: Path) -> Path: return outfile except requests.HTTPError as e: - print(f"Failed to download genomes: {e}") + print(f"Failed to download genomes for accessions: {accessions}") + sys.exit(f"Reason: {e}") From cb7b5ba7f31927596c17e464f989ab8302916801 Mon Sep 17 00:00:00 2001 From: Jim Downie <19718667+prototaxites@users.noreply.github.com> Date: Thu, 9 Oct 2025 15:06:22 +0100 Subject: [PATCH 03/13] fix: minor change to file to trigger ci --- .github/workflows/docker-publish.yml | 1 - 1 file changed, 1 deletion(-) diff --git a/.github/workflows/docker-publish.yml b/.github/workflows/docker-publish.yml index 2f1fc70..fd056be 100644 --- a/.github/workflows/docker-publish.yml +++ b/.github/workflows/docker-publish.yml @@ -18,7 +18,6 @@ permissions: jobs: build: - runs-on: ubuntu-latest steps: From 73d55cc180ade1f3c50b5a673adfcd0ccebe8a73 Mon Sep 17 00:00:00 2001 From: Jim Downie <19718667+prototaxites@users.noreply.github.com> Date: Tue, 14 Oct 2025 10:47:01 +0100 Subject: [PATCH 04/13] fix(dockerfile): force tabulate 0.8.10 for compatibility with fixed version of snakemake --- src/docker/Dockerfile | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/docker/Dockerfile b/src/docker/Dockerfile index db1153d..cc83b60 100644 --- a/src/docker/Dockerfile +++ b/src/docker/Dockerfile @@ -2,7 +2,7 @@ FROM condaforge/miniforge3:24.9.2-0 ENV SHELL=/bin/bash -RUN mamba create -q -y -c conda-forge -c bioconda -n snakemake snakemake=6.10.0 python=3.9.7 singularity \ +RUN mamba create -q -y -c conda-forge -c bioconda -n snakemake snakemake=6.10.0 python=3.9.7 singularity tabulate=0.8.10 \ && conda clean --all -y RUN echo "source activate snakemake" > ~/.bashrc From e2db290507f5f6272890b58173da9db32aab84e8 Mon Sep 17 00:00:00 2001 From: Jim Downie <19718667+prototaxites@users.noreply.github.com> Date: Tue, 14 Oct 2025 11:37:59 +0100 Subject: [PATCH 05/13] fix(ncbi): run get on json() object, not requests response directly --- scripts/NCBIApiTools.py | 13 +++++++------ 1 file changed, 7 insertions(+), 6 deletions(-) diff --git a/scripts/NCBIApiTools.py b/scripts/NCBIApiTools.py index ea0bb9b..6957eac 100644 --- a/scripts/NCBIApiTools.py +++ b/scripts/NCBIApiTools.py @@ -46,19 +46,20 @@ def get_assemblies_for_taxon( arg_filter_refs_only = f"?filters_reference_only={filters_reference_only}" arg_page_size = f"&page_size={page_size}" r = self.get_request(query + arg_filter_refs_only + arg_page_size) + r_res = r.json() assemblies = [] - if r.get("reports"): - assemblies = [] + r.get("reports") + if r_res.get("reports"): + assemblies = [] + r_res.get("reports") - while r.get("next_page_token", "") != "": - page_token_arg = f"&page_token={r.get('next_page_token')}" + while r_res.get("next_page_token", "") != "": + page_token_arg = f"&page_token={r_res.get('next_page_token')}" r = self.get_request( query + arg_filter_refs_only + arg_page_size + page_token_arg ) - if r.get("reports", None): - assemblies = assemblies + r.get("reports") + if r_res.get("reports", None): + assemblies = assemblies + r_res.get("reports") return assemblies From 1b841a960e4e803d706088cfa8fde7abe30a87ed Mon Sep 17 00:00:00 2001 From: Jim Downie Date: Tue, 14 Oct 2025 12:42:11 +0100 Subject: [PATCH 06/13] fix(ncbi): don't fail on nested gets --- scripts/FetchGenomesRefSeq.py | 4 +++- scripts/FetchGenomesRefSeqRelatives.py | 4 +++- 2 files changed, 6 insertions(+), 2 deletions(-) diff --git a/scripts/FetchGenomesRefSeq.py b/scripts/FetchGenomesRefSeq.py index 30f60a1..50a0d83 100644 --- a/scripts/FetchGenomesRefSeq.py +++ b/scripts/FetchGenomesRefSeq.py @@ -59,7 +59,9 @@ def Average(lst): acc = assembly.get("accession") date = assembly.get("release_date") sciname_orig = assembly.get("organism").get("organism_name") - strainname = assembly.get("organism").get("infraspecific_names").get("strain") + strainname = ( + assembly.get("organism", {}).get("infraspecific_names", {}).get("strain", None) + ) contiguity = int(assembly.get("assembly_stats").get("contig_n50")) size = int(assembly.get("assembly_stats").get("total_sequence_length")) diff --git a/scripts/FetchGenomesRefSeqRelatives.py b/scripts/FetchGenomesRefSeqRelatives.py index e4e04e3..4890640 100644 --- a/scripts/FetchGenomesRefSeqRelatives.py +++ b/scripts/FetchGenomesRefSeqRelatives.py @@ -129,7 +129,9 @@ def readNodes(nodes_tax_file): # print(assembly.org.sci_name) sciname_orig = assembly.get("organism").get("organism_name") strainname = ( - assembly.get("organism").get("infraspecific_names").get("strain") + assembly.get("organism", {}) + .get("infraspecific_names", {}) + .get("strain", None) ) if strainname: if strainname in sciname_orig and "sp" not in sciname_orig: From a7429a73766bb4e139cb396dfe90134664d3eb7f Mon Sep 17 00:00:00 2001 From: Jim Downie Date: Tue, 14 Oct 2025 13:19:46 +0100 Subject: [PATCH 07/13] fix(ncbi): don't fail on nested gets --- scripts/FetchGenomesRefSeqRelatives.py | 6 +++++- 1 file changed, 5 insertions(+), 1 deletion(-) diff --git a/scripts/FetchGenomesRefSeqRelatives.py b/scripts/FetchGenomesRefSeqRelatives.py index 4890640..dfc4311 100644 --- a/scripts/FetchGenomesRefSeqRelatives.py +++ b/scripts/FetchGenomesRefSeqRelatives.py @@ -151,7 +151,11 @@ def readNodes(nodes_tax_file): acc = assembly.get("accession") date = assembly.get("release_date") sciname_orig = assembly.get("organism").get("organism_name") - strainname = assembly.get("organism").get("infraspecific_names").get("strain") + strainname = ( + assembly.get("organism", {}) + .get("infraspecific_names", {}) + .get("strain", None) + ) contiguity = int(assembly.get("assembly_stats").get("contig_n50")) size = int(assembly.get("assembly_stats").get("total_sequence_length")) From 6ea527d7f9cfd00c22c2587e3a613d1621276fb0 Mon Sep 17 00:00:00 2001 From: Jim Downie Date: Tue, 14 Oct 2025 14:27:27 +0100 Subject: [PATCH 08/13] fix(ncbi): get date correctly --- .gitignore | 1 + scripts/FetchGenomesRefSeq.py | 2 +- scripts/FetchGenomesRefSeqRelatives.py | 2 +- 3 files changed, 3 insertions(+), 2 deletions(-) create mode 100644 .gitignore diff --git a/.gitignore b/.gitignore new file mode 100644 index 0000000..65e3ba2 --- /dev/null +++ b/.gitignore @@ -0,0 +1 @@ +test/ diff --git a/scripts/FetchGenomesRefSeq.py b/scripts/FetchGenomesRefSeq.py index 50a0d83..500c78e 100644 --- a/scripts/FetchGenomesRefSeq.py +++ b/scripts/FetchGenomesRefSeq.py @@ -57,7 +57,7 @@ def Average(lst): SpeciesDictionary = {} for assembly in assemblies: acc = assembly.get("accession") - date = assembly.get("release_date") + date = assembly.get("assembly_info", {}).get("release_date") sciname_orig = assembly.get("organism").get("organism_name") strainname = ( assembly.get("organism", {}).get("infraspecific_names", {}).get("strain", None) diff --git a/scripts/FetchGenomesRefSeqRelatives.py b/scripts/FetchGenomesRefSeqRelatives.py index dfc4311..ea55898 100644 --- a/scripts/FetchGenomesRefSeqRelatives.py +++ b/scripts/FetchGenomesRefSeqRelatives.py @@ -149,7 +149,7 @@ def readNodes(nodes_tax_file): SpeciesDictionary = {} for assembly in assemblies: acc = assembly.get("accession") - date = assembly.get("release_date") + date = assembly.get("assembly_info", {}).get("release_date") sciname_orig = assembly.get("organism").get("organism_name") strainname = ( assembly.get("organism", {}) From 29d95c6683c085f5c0c0b1f83a268fa1cbb1d110 Mon Sep 17 00:00:00 2001 From: Jim Downie Date: Tue, 14 Oct 2025 15:20:48 +0100 Subject: [PATCH 09/13] fix(ncbi): correct path commands --- .gitignore | 2 ++ scripts/FetchGenomesRefSeq.py | 19 +++---------------- scripts/FetchGenomesRefSeqRelatives.py | 13 +++---------- 3 files changed, 8 insertions(+), 26 deletions(-) diff --git a/.gitignore b/.gitignore index 65e3ba2..6cc6356 100644 --- a/.gitignore +++ b/.gitignore @@ -1 +1,3 @@ test/ +**.pyc +**/__pycache__/ diff --git a/scripts/FetchGenomesRefSeq.py b/scripts/FetchGenomesRefSeq.py index 500c78e..22fdacc 100644 --- a/scripts/FetchGenomesRefSeq.py +++ b/scripts/FetchGenomesRefSeq.py @@ -112,15 +112,13 @@ def Average(lst): zipfile_name = str(args.dir) + "/" + "/RefSeq." + str(taxname_orig) + ".zip" for t in range(0, len(accs), 100): accshort = accs[t : t + 100] - zipfile_name_part = ( - Path(args.dir) / f"/RefSeq.{str(taxname_orig)}.part{str(t)}.zip" - ) + zipfile_name_part = Path(args.dir) / f"/{taxname_orig}.RefSeq.part{str(t)}.zip" api_response = api_instance.download_genomes( accshort, outfile=zipfile_name_part, ) print("Download complete part " + str(t)) - cmd = "unzip -d " + str(zipfile_name_part)[-4] + " " + str(zipfile_name_part) + cmd = f"unzip -d {zipfile_name_part.stem} {zipfile_name_part}" os.system(cmd) cmd = "mkdir " + str(args.dir) + "/" + str(taxname_orig) + ".Refseq" @@ -167,18 +165,7 @@ def Average(lst): + ".Refseq/ncbi_dataset/data/dataset_catalog.json" ) os.system(cmd) - cmd = ( - "rm -r " - + str(args.dir) - + "/" - + "/RefSeq." - + str(taxname_orig) - + ".part*.zip " - + str(args.dir) - + "/" - + str(taxname_orig) - + ".RefSeq.part*" - ) + cmd = "rm -r " + str(args.dir) + f"/{str(taxname_orig)}.RefSeq.part*" os.system(cmd) else: print("No genomes available") diff --git a/scripts/FetchGenomesRefSeqRelatives.py b/scripts/FetchGenomesRefSeqRelatives.py index ea55898..347553c 100644 --- a/scripts/FetchGenomesRefSeqRelatives.py +++ b/scripts/FetchGenomesRefSeqRelatives.py @@ -196,14 +196,14 @@ def readNodes(nodes_tax_file): print(f"Download a package for {accs}.") print("Begin download of genome data package ...") for t in range(0, len(accs), 3): - zipfile_name_part = Path(str(args.dir)) / f"RefSeq.relatives.part{str(t)}.zip" + zipfile_name_part = Path(str(args.dir)) / f"relatives.RefSeq.part{str(t)}.zip" accshort = accs[t : t + 3] api_response = api_instance.download_genomes( accshort, outfile=zipfile_name_part, ) print("Download complete part " + str(t)) - cmd = "unzip -d " + str(zipfile_name_part)[-4] + " " + str(zipfile_name_part) + cmd = f"unzip -d {zipfile_name_part.stem} {zipfile_name_part}" os.system(cmd) cmd = "mkdir " + str(args.dir) + "/relatives.Refseq" @@ -236,14 +236,7 @@ def readNodes(nodes_tax_file): + "/relatives.Refseq/ncbi_dataset/data/dataset_catalog.json" ) os.system(cmd) - cmd = ( - "rm -r " - + str(args.dir) - + "/" - + "/RefSeq.relatives.part*.zip " - + str(args.dir) - + "/relatives.RefSeq.part*" - ) + cmd = "rm -r " + str(args.dir) + "/" + "relatives.RefSeq.part*" os.system(cmd) cmd = ( "cat " From 9c5b1c7c1a0b2fe5f76bbe614943b1e57e91ac63 Mon Sep 17 00:00:00 2001 From: Jim Downie Date: Wed, 15 Oct 2025 09:43:24 +0100 Subject: [PATCH 10/13] fix(envs): bump minimap version to 2.30 --- envs/minimap.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/envs/minimap.yaml b/envs/minimap.yaml index 8d21896..ff0b3d7 100644 --- a/envs/minimap.yaml +++ b/envs/minimap.yaml @@ -2,5 +2,5 @@ name: minimap channels: - bioconda dependencies: - - minimap2=2.17 + - minimap2=2.30 - seqtk=1.3 From a9dbb45b8e63c880c7053dac4ed6a6afc74f1628 Mon Sep 17 00:00:00 2001 From: Jim Downie Date: Wed, 15 Oct 2025 13:05:56 +0100 Subject: [PATCH 11/13] fix(envs): bump seqtk version to 1.5 --- envs/minimap.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/envs/minimap.yaml b/envs/minimap.yaml index ff0b3d7..1a1f5c6 100644 --- a/envs/minimap.yaml +++ b/envs/minimap.yaml @@ -3,4 +3,4 @@ channels: - bioconda dependencies: - minimap2=2.30 - - seqtk=1.3 + - seqtk=1.5 From cd8dada6e685a6cb19c8f0a53b759bf9b9c388b0 Mon Sep 17 00:00:00 2001 From: Jim Downie Date: Wed, 15 Oct 2025 16:02:04 +0100 Subject: [PATCH 12/13] fix(envs): bump busco to 6.0.0 --- envs/busco.yaml | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/envs/busco.yaml b/envs/busco.yaml index dc12654..9d82a17 100644 --- a/envs/busco.yaml +++ b/envs/busco.yaml @@ -3,4 +3,4 @@ channels: - conda-forge - bioconda dependencies: - - busco=5.2.2 \ No newline at end of file + - busco=6.0.0 From e83f86066b04e63cde4b739599b76a8b77845923 Mon Sep 17 00:00:00 2001 From: Jim Downie Date: Wed, 15 Oct 2025 16:41:47 +0100 Subject: [PATCH 13/13] fix(busco_config): don't try and guess the paths --- scripts/BuscoConfig.py | 229 +++++++++++++++++++++++++---------------- 1 file changed, 138 insertions(+), 91 deletions(-) diff --git a/scripts/BuscoConfig.py b/scripts/BuscoConfig.py index 20245a4..58fe4d6 100644 --- a/scripts/BuscoConfig.py +++ b/scripts/BuscoConfig.py @@ -5,76 +5,121 @@ import sys parser = argparse.ArgumentParser() -parser.add_argument("-na", type=str, action='store', dest='namesfile', metavar='NAMES',help='NCBI names.dmp') -parser.add_argument("-no", type=str, action='store', dest='nodesfile', metavar='NODES',help='NCBI nodes.dmp') -parser.add_argument("-f", type=str, action='store', dest='genome', metavar='GENOME FASTA',help='fasta genome assembly file') -parser.add_argument("-d", type=str, action='store', dest='dir', metavar='WORKDIR',help='define working directory for busco') -parser.add_argument("-db", type=str, action='store', dest='db', help='define available dbs file') -parser.add_argument("-dl", type=str, action='store', dest='download',help='define directory to store busco dbs') -parser.add_argument("-c", type=int, action='store', dest='cpu',help='define cpus') -parser.add_argument("-o", type=str, action='store', dest='out', metavar='OUTFILE',help='define configfile name') -parser.add_argument('--version', action='version', version='%(prog)s 1.0') +parser.add_argument( + "-na", + type=str, + action="store", + dest="namesfile", + metavar="NAMES", + help="NCBI names.dmp", +) +parser.add_argument( + "-no", + type=str, + action="store", + dest="nodesfile", + metavar="NODES", + help="NCBI nodes.dmp", +) +parser.add_argument( + "-f", + type=str, + action="store", + dest="genome", + metavar="GENOME FASTA", + help="fasta genome assembly file", +) +parser.add_argument( + "-d", + type=str, + action="store", + dest="dir", + metavar="WORKDIR", + help="define working directory for busco", +) +parser.add_argument( + "-db", type=str, action="store", dest="db", help="define available dbs file" +) +parser.add_argument( + "-dl", + type=str, + action="store", + dest="download", + help="define directory to store busco dbs", +) +parser.add_argument("-c", type=int, action="store", dest="cpu", help="define cpus") +parser.add_argument( + "-o", + type=str, + action="store", + dest="out", + metavar="OUTFILE", + help="define configfile name", +) +parser.add_argument("--version", action="version", version="%(prog)s 1.0") args = parser.parse_args() + def readNames(names_tax_file): - ''' + """ input: - name.dmp (NCBI Taxonomy) output: - dictionary of form {node: name} - dictionary of form {sci name: node} - ''' + """ tax_names = {} - tax_names_reverse= {} - with open(names_tax_file, 'r') as nodes_tax: + tax_names_reverse = {} + with open(names_tax_file, "r") as nodes_tax: for line in nodes_tax: - node = [field.strip() for field in line.split('|')] - if 'scientific' in line: + node = [field.strip() for field in line.split("|")] + if "scientific" in line: tax_names[node[1]] = node[0] tax_names_reverse[node[0]] = node[1] - return tax_names_reverse,tax_names + return tax_names_reverse, tax_names -def readNodes(nodes_tax_file): - ''' +def readNodes(nodes_tax_file): + """ input: - nodes.dmp (NCBI Taxonomy) output: - dictionary of form {parent: node} - dictionary of form {node: type} - ''' + """ tax_nodes = {} tax_types = {} - with open(nodes_tax_file, 'r') as nodes_tax: + with open(nodes_tax_file, "r") as nodes_tax: for line in nodes_tax: - node = [field.strip() for field in line.split('|')] #make list of line - tax_nodes[node[0]] = node[1] #couple node with parent - tax_types[node[0]] = node[2] #couple node with rank + node = [field.strip() for field in line.split("|")] # make list of line + tax_nodes[node[0]] = node[1] # couple node with parent + tax_types[node[0]] = node[2] # couple node with rank return tax_nodes -taxparents=readNodes(args.nodesfile) -taxnames,namestax=readNames(args.namesfile) -genus=args.out.split('/config')[0].split('/')[-1] -if '_' in genus: - genus=genus.replace('_',' ') +taxparents = readNodes(args.nodesfile) +taxnames, namestax = readNames(args.namesfile) -busco_dbs=[] -busco_short=[] -m=open(args.db,'r') +genus = args.out.split("/config")[0].split("/")[-1] +if "_" in genus: + genus = genus.replace("_", " ") + +busco_dbs = [] +busco_short = [] +m = open(args.db, "r") for line in m: - line=line.strip() - if 'db' in line: - #if 'eukaryota' in line: + line = line.strip() + if "db" in line: + # if 'eukaryota' in line: # break - dbname=line.split(' ')[-1].split('_odb')[0] - #print(dbname) - dbname2=line.split(' ')[-1].split('_')[0] + dbname = line.split(" ")[-1].split("_odb")[0] + # print(dbname) + dbname2 = line.split(" ")[-1].split("_")[0] busco_short.append(dbname2) busco_dbs.append(dbname) -buscoset = 'Bacteria' +buscoset = "Bacteria" if genus in namestax: taxid = namestax[genus] parent = taxparents[taxid] @@ -83,61 +128,63 @@ def readNodes(nodes_tax_file): if taxnames[parent].lower() in busco_short: buscoset = taxnames[parent] print(buscoset) - if taxnames[parent].lower()+"_phylum" in busco_dbs: - buscoset = taxnames[parent].lower()+"_phylum" + if taxnames[parent].lower() + "_phylum" in busco_dbs: + buscoset = taxnames[parent].lower() + "_phylum" break parent = taxparents[parent] -print(genus+'\t'+buscoset) +print(genus + "\t" + buscoset) + +condadir = os.environ["CONDA_DEFAULT_ENV"] -condadir = os.environ['CONDA_DEFAULT_ENV'] +config_content = f"""[busco_run] +# Input file +in = {args.genome} +# Run name, used in output files and folder +out = busco +# Where to store the output directory +out_path = {args.dir} +# Path to the BUSCO dataset +lineage_dataset = {buscoset.lower()} +# Which mode to run (genome / proteins / transcriptome) +mode = genome +# How many threads to use for multithreaded steps +cpu = {args.cpu} +# Force rewrite if files already exist (True/False) +;force = False +# Local destination path for downloaded lineage datasets +download_path = {args.download} +;[tblastn] +;path = {condadir}/bin/ +;command = tblastn +;[makeblastdb] +;path = {condadir}/bin/ +;command = makeblastdb +;[augustus] +;path = {condadir}/bin/ +;command = augustus +;[etraining] +;path = {condadir}/bin/ +;command = etraining +;[gff2gbSmallDNA.pl] +;path = {condadir}/bin/ +;command = gff2gbSmallDNA.pl +;[new_species.pl] +;path = {condadir}/bin/ +;command = new_species.pl +;[optimize_augustus.pl] +;path = {condadir}/bin/ +;command = optimize_augustus.pl +;[hmmsearch] +;path = {condadir}/bin/ +;command = hmmsearch +;[sepp] +;path = {condadir}/bin/ +;command = run_sepp.py +;[prodigal] +;path = {condadir}/bin/ +;command = prodigal +""" -l=open(args.out,'w') -l.write('[busco_run]'+'\n') -l.write('# Input file'+'\n') -l.write('in = '+args.genome+'\n') -l.write('# Run name, used in output files and folder'+'\n') -l.write('out = busco'+'\n') -l.write('# Where to store the output directory'+'\n') -l.write('out_path = '+args.dir+'\n') -l.write('# Path to the BUSCO dataset'+'\n') -l.write('lineage_dataset = '+buscoset.lower()+'\n') -l.write('# Which mode to run (genome / proteins / transcriptome)'+'\n') -l.write('mode = genome'+'\n') -l.write('# How many threads to use for multithreaded steps'+'\n') -l.write('cpu = '+str(args.cpu)+'\n') -l.write('# Force rewrite if files already exist (True/False)'+'\n') -l.write(';force = False'+'\n') -l.write('# Local destination path for downloaded lineage datasets'+'\n') -l.write('download_path = '+args.download+'\n') -l.write('[tblastn]'+'\n') -l.write('path = '+condadir+'/bin/'+'\n') -l.write('command = tblastn'+'\n') -l.write('[makeblastdb]'+'\n') -l.write('path = '+condadir+'/bin/'+'\n') -l.write('command = makeblastdb'+'\n') -l.write('[augustus]'+'\n') -l.write('path = '+condadir+'/bin/'+'\n') -l.write('command = augustus'+'\n') -l.write('[etraining]'+'\n') -l.write('path = '+condadir+'/bin/'+'\n') -l.write('command = etraining'+'\n') -l.write('[gff2gbSmallDNA.pl]'+'\n') -l.write('path = '+condadir+'/bin/'+'\n') -l.write('command = gff2gbSmallDNA.pl'+'\n') -l.write('[new_species.pl]'+'\n') -l.write('path = '+condadir+'/bin/'+'\n') -l.write('command = new_species.pl'+'\n') -l.write('[optimize_augustus.pl]'+'\n') -l.write('path = '+condadir+'/bin/'+'\n') -l.write('command = optimize_augustus.pl'+'\n') -l.write('[hmmsearch]'+'\n') -l.write('path = '+condadir+'/bin/'+'\n') -l.write('command = hmmsearch'+'\n') -l.write('[sepp]'+'\n') -l.write('path = '+condadir+'/bin/'+'\n') -l.write('command = run_sepp.py'+'\n') -l.write('[prodigal]'+'\n') -l.write('path = '+condadir+'/bin/'+'\n') -l.write('command = prodigal'+'\n') -l.close() +with open(args.out, "w") as l: + l.write(config_content)