Use taxid to get BUSCO lineage

[FernandoDuarteF]

Description of feature

If taxid column is given for a specific assembly, use this instead of params.busco_lineage, as done in sanger's blobtoolkit.

[LiaOb21]

@FernandoDuarteF I didn't really get how it works in Sanger's blobtoolkit. But here is my plan:

1. Pull taxonomy using ncbi datasets with a command like this:

datasets summary taxonomy taxon 36056 --as-json-lines | dataformat tsv taxonomy --template tax-summary | cut -f'2,3,11,13,15,19,21'
Taxid	Tax name	Domain/Realm name	Kingdom name	Phylum name	Order name	Family name
36056	Boletus edulis	Eukaryota	Fungi	Basidiomycota	Boletales	Boletaceae

or:

datasets summary taxonomy taxon --inputfile txids.txt --as-json-lines | dataformat tsv taxonomy --template tax-summary | cut -f'2,3,11,13,15,17,19,21'
Taxid	Tax name	Domain/Realm name	Kingdom name	Phylum name	Class name	Order name	Family name
36056	Boletus edulis	Eukaryota	Fungi	Basidiomycota	Agaricomycetes	Boletales	Boletaceae
3702	Arabidopsis thaliana	Eukaryota	Viridiplantae	Streptophyta	Magnoliopsida	Brassicales	Brassicaceae
4577	Zea mays	Eukaryota	Viridiplantae	Streptophyta	Magnoliopsida	Poales	Poaceae

With a bit of changes to adapt (temporarily) to my script:

datasets summary taxonomy taxon --inputfile txids.txt --as-json-lines | dataformat tsv taxonomy --template tax-summary | cut -f'2,3,11,13,15,17,19,21' | sed 's/ name//g'| sed 's/Tax/Sample/g' > samples.tsv

samples.tsv:

Sampleid	Sample	Domain/Realm	Kingdom	Phylum	Class	Order	Family
36056	Boletus edulis	Eukaryota	Fungi	Basidiomycota	Agaricomycetes	Boletales	Boletaceae
3702	Arabidopsis thaliana	Eukaryota	Viridiplantae	Streptophyta	Magnoliopsida	Brassicales	Brassicaceae
4577	Zea mays	Eukaryota	Viridiplantae	Streptophyta	Magnoliopsida	Poales	Poaceae

We'll need a module for ncbi datasets. There is already entrez in nf-core modules, but the only commands available are esearch, esummary, xtract, and I think the command to do this should be efetch.

2. use busco --list > busco_lineages.txt to obtain the list of available databases
3. extract busco databases of interest:

gawk '/eukaryota_odb12/{flag=1; indent=length($0)-length(ltrim($0)); print "eukaryota_odb12"; next}
     flag && /- [a-z_]*_odb12/ {
         current_indent=length($0)-length(ltrim($0))
         if(current_indent <= indent) flag=0
         else print gensub(/.*- ([a-z_]*_odb12).*/, "\\1", "g")
     }
     function ltrim(s) { sub(/^[ \t\r\n]+/, "", s); return s }' busco_lineages.txt > eukaryota_busco_lineages.txt

4. Use my script to map the sample to the closest dataset:

python FSP_assembly_benchmarking/workflow/scripts/map_busco_db_to_sample.py -s "Arabidopsis thaliana" -t samples.tsv -b eukaryota_busco_lineages.txt -o output.tsv -e "_odb12"
Arabidopsis thaliana: brassicales_odb12

python FSP_assembly_benchmarking/workflow/scripts/map_busco_db_to_sample.py -s "Boletus edulis" -t samples.tsv -b eukaryota_busco_lineages.txt -o output.tsv -e "_odb12"
Boletus edulis: boletales_odb12

python FSP_assembly_benchmarking/workflow/scripts/map_busco_db_to_sample.py -s "Zea mays" -t samples.tsv -b eukaryota_busco_lineages.txt -o output.tsv -e "_odb12"
Zea mays: poales_odb12

I don't know if there is an easier way to achieve this. Let me know what you think! :)

[FernandoDuarteF]

Could also use this to search for telomeres using the tidk database.

[FernandoDuarteF]

This might be more complicated than i initially thought.

Will see if I can do it for tomorrow, if not might be worth to just wait for a second release.