ASVs missing from output

[owen-hughes42]

Hello

I'm running taxonomizr to assign taxonomy of some BLASTed sequences, but the output is missing some ASVs, as in they don't even come up as NA for assignment.

I'm running based on this method: [https://github.com/ewan-harney/hpc_blast2taxonomizr?tab=readme-ov-file] but adapted for my university's HPC.

Here's the code taxonomizr code:
`.libPaths(c( .libPaths(), "/mnt/iusers01/fatpou01/bmh01/ug-bioinf-2025-2026/d93291oh/R/x86_64-pc-linux-gnu-library/4.5"))

Load libraries

library(taxonomizr)
library(dplyr)
library(optparse)

parse command line options

option_list = list(
make_option(c("-E", "--email"), type="character", default=NULL,
help="Provide an email address to receive an email notification when the job has finished.", metavar="character"),
make_option(c("-T", "--topperc"), type="numeric", default=NULL,
help="Top percent value for taxonomizr", metavar="character"))

Parse arguments

opt_parser = OptionParser(option_list=option_list)
opt = parse_args(opt_parser)

Get the directory containing the raw files

path<-getwd()

Read in the filtered_blast.out.tab file

blastResults<-read.table(file = paste(path, "/filtered_blast.out.tab", sep=""),header=FALSE,stringsAsFactors=FALSE)

Filter file by (user-supplied) top percent value of bit score

Note that the top percent value is converted to a fraction and subtracted from 1:

thus a top percent value of 2 becomes 0.98

topPercent <- data.frame(blastResults %>% group_by(V1) %>%
arrange(V1, desc(V12)) %>%
filter(V12 > max(V12*(1 - (opt$topperc/100)))))

Split the top percent values in to asvs and accessions

asvs<-topPercent[,1]
accessions<-topPercent[,2]

Get the taxonomic information for top percent accessions from accessionTaxa.sql

taxaId<-accessionToTaxa(accessions,"./../accessionTaxa.sql")
taxa<-getTaxonomy(taxaId,'./../accessionTaxa.sql', desiredTaxa=c("kingdom", "phylum", "class", "order", "family", "genus", "species"))

LCA algorithm and merge back to top percent asvs

taxon_path <- condenseTaxa(taxa,asvs)

Write output

write.table(taxon_path, "taxonomizr_taxon_path.tsv", sep="\t", quote = FALSE)

Email user

email_plot_command <- paste("echo "Taxonomizr is complete." | mail -s "Taxonomizr table" -a taxonomizr_taxon_path.tsv", opt$email, sep=" ")
system(email_plot_command)`

My parameters are as follows: -P 97 -L 70 -T 5

My blast_out file contains 1235 ASVs, but the taxonomizr_taxon_path.tsv contains 950

Any help would be massively appreciated!

Thanks
Owen

[sherrillmix]

Hi Owen. Sorry just saw this. Hope everything worked out. Let me know if any additional issues.

[owen-hughes42]

Hey! No worries.

I'm having another related issue (background is all the same, except T is now 95.

The taxonomizr output seems to assign lots of ASVs "NA" identity even though the blast output has taxa with high percent identity values etc. I've tried changing the parameters around lots but still can't find one that fixes things.

Here's the blast output for ASV20:

It's also still there in the filtered blast_out.

But the taxonomizr output gives it NA for all levels of assignment.

Any ideas as to what I'm doing wrong would be massively appreciated!

Thanks
Owen :)

[sherrillmix]

Can you specify an accession number and/or taxa ID that resulted in NA that you think is inappropriate?

[owen-hughes42]

Did the screenshot in the previous message upload? That ASV has lots of blast hits with high pident and bit scores (e.g the top one with ID 4545) but got NA for all taxonomic levels, even ones that are shared across all the hits like Eukaryota.
Sorry if I’m not explaining properly I’m quite inexperienced.

[sherrillmix]

Yeah I kind of guessed the experience part. :) A random google hit on how to file a bug report:
https://www.chiark.greenend.org.uk/~sgtatham/bugs.html
That might be helpful.

I see a screenshot but it is of a table not produced by taxonomizr so not particularly helpful and I don't see any NAs. So I don't really have anyway to replicate your problem or debug anything. In general, screenshots are often suboptimal since you're asking someone to spend their time trying to figure out what you're asking about or making them manually enter text from the picture.

[owen-hughes42]

Haha I'll try my best to be more clear/specific :)

Here's a a tabular file with the accession numbers etc of the BLAST hits for ASV20 that is used for the input for taxonomizr (blastResults in above code):

asv20.xlsx

Here's the output from the condenseTaxa function (taxon path in above code) for ASV20:

output.xlsx

I hope that's more clear :)

[sherrillmix]

None of these are taxonomizr output so you should probably be checking with whoever wrote the pipeline that is generating them. Especially since I don't know what any of the columns are.

But trying to be helpful and read between the lines, is there an intermediate file between those two steps? If it shows taxonomies then you might want to check if you're getting hits to two different kingdoms e.g. cyanobacteria and plants. Then an all NA output would be reasonable since the common ancestor is upstream of kingdom. I guess your screenshot example also apparently shows kingdom (although I'm not sure if that's taxonomizr derived or something else) so could scroll through there.

[sherrillmix]

Alternatively, it could be that you're getting accession numbers where taxa aren't assigned and the pipeline isn't excluding them. For example using an old accession database that doesn't have a taxa ID for recent upload PX975919:

> taxId=taxonomizr::accessionToTaxa(c('GQ324408.1','PX975919.1'),'OLD_DATABASE.sqlite')
> taxId
[1] 764749     NA

> taxa=taxonomizr::getTaxonomy(taxId,'db/taxo/nameNode.sqlite',desiredTaxa=c('kingdom',"phylum", "class", "order", "family", "genus", "species"))
> taxa
          kingdom         phylum         class           order    family    genus species        
764749 "Viridiplantae" "Streptophyta" "Magnoliopsida" "Poales" "Poaceae" "Poa" "Poa clelandii"
    NA NA              NA             NA              NA       NA        NA    NA       
> taxonomizr::condenseTaxa(taxa)
  kingdom phylum class order family genus species
1 NA      NA     NA    NA    NA     NA    NA

Whoever wrote your pipeline might want to exclude anything with an NA for kingdom. They would have to be careful though since I believe many bacteria don't have kingdom annotated or they could pull an additional domain column and filter out anything with NAs there. Or just exclude any NA taxa IDs if that's the root of the problem.

[owen-hughes42]

thanks so much I think that's the solution!
is there a way to exclude NAs with the taxonomizr functions ?

[sherrillmix]

Just using R's selection functions should do the job, so something like taxId[!is.na(taxId)]. So the above code would be:

> taxa=taxonomizr::getTaxonomy(taxId[!is.na(taxId)],'db/taxo/nameNode.sqlite',desiredTaxa=c('kingdom',"phylum", "class", "order", "family", "genus", "species"))
> taxa
       kingdom         phylum         class           order    family    genus species        
764749 "Viridiplantae" "Streptophyta" "Magnoliopsida" "Poales" "Poaceae" "Poa" "Poa clelandii"

or in your original post:

taxaId<-accessionToTaxa(accessions,"./../accessionTaxa.sql")
taxa<-getTaxonomy(taxaId[!is.na(taxaId)],'./../accessionTaxa.sql', desiredTaxa=c("kingdom", "phylum", "class", "order", "family", "genus", "species"))

Also some of these NAs might be related to an old taxonomy database not containing info for newly deposited sequences so redownloading might help.

[owen-hughes42]

amazing thank you so much!

do you know of a similar way to filter out anything assigned NA at kingdom level instead of just the missing IDs? I've tried lots of things and nothing seems to work

no worries if not you've been so helpful/patient so far!

[sherrillmix]

Square brackets select in the format matrix[rows,columns] where rows/columns can be specified by row/column names, numbers or a TRUE/FALSE vector. is.na() returns TRUE where it's argument is NA and FALSE otherwise. ! is a "not" reversing trues and falses. You want the rows in the taxa matrix (taxa[XXX,]) where the kingdom column (taxa[,'kingdom']) of taxa is not NA (!is.na()) so:

taxa=taxa[!is.na(taxa[,'kingdom']),]

should do it.

[owen-hughes42]

got it to work thank you so much!

I had to merge the taxa with the asvs to filter the asvs based on the assignment as well but its worked now :)

really appreciate all the help
owen :)