1
$\begingroup$

I would like to perform blast+, using the command line, here is the code:

-task blastn -db nr -outfmt 7 -query myQuery.fas -out result.txt -remote.

I have to exclude SARS-CoV-2 hits from my results, so I must include -negative_taxids mentioning SARS-CoV-2 NCBI taxid number. However, the problem is, -negative_taxids is NOT compatible with -remote command and my command line encounters an error. I have to use -remote as I cannot download entire +400 GB NR database to run it locally. Could you please kindly help me out with this contradiction in my work?

My output txt file for this command line is:

# Fields: query acc.ver, subject acc.ver, % identity, alignment length, mismatches, gap opens, q. start, q. end, s. start, s. end, evalue, bit score
# 600483 hits found
AstV-VA1    MN464146.1  100.000 32  0   0   1   32  6362    6393    7.14e-06    59.0
AstV-VA1    MN148428.1  100.000 32  0   0   1   32  6685    6716    7.14e-06    59.0
AstV-VA1    MK671313.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0
AstV-VA1    MK671312.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0
AstV-VA1    MK671311.1  100.000 32  0   0   1   32  6682    6713    7.14e-06    59.0
$\endgroup$
5
  • $\begingroup$ Have you considered running your query against all of nr and then removing the SARS hits from the results instead? $\endgroup$
    – terdon
    Nov 3, 2023 at 10:43
  • $\begingroup$ Thanks terdon for your constant and positive contribution to my questions. Actually, I have already done it and I got +600,000 hits but the point is the txt result file contains only the following information: "query", "acc.ver", "% identity", "alignment length", "mismatches", "gap opens", "q. start", "q. end", "s. start", "s. end", "Evalue", "bit score" so, there are no organism and taxid columns in the result. One possible solution can be using R to download the annotation to add the results automatically, but it has two disadvantages: $\endgroup$ Nov 3, 2023 at 11:54
  • $\begingroup$ 1. downloading annotation of +600,000 hits will take around 200 hours (more than a week) 2. not all hits automatically yield annotations. Therefore, I think filtering the results from the beginning can be a better option. $\endgroup$ Nov 3, 2023 at 11:54
  • 1
    $\begingroup$ No, no, there is no reason to download entire annotations! All you need to do is extract the IDs of the target sequences and then get the taxids for those sequences. If you edit your question and show us a few examples of your output, we should be able to help, but it will just be a simple entrez query whose input is the list of sequence identifiers and output will be the associated taxids. $\endgroup$
    – terdon
    Nov 3, 2023 at 11:58
  • $\begingroup$ What do you want to do with these BLAST hits? $\endgroup$
    – gringer
    Nov 5, 2023 at 11:25

2 Answers 2

1
$\begingroup$

This isn't particularly elegant, but it should work. First, extract the accessions into a file. Next, get the taxid for each of the accessions. Once you have that, you can identify the accessions that are for SARS-CoV-2 sequences, and remove them from your output. Something like this:

  1. Extract the accessions into a new file

    awk '!/^#/{print $2}' blast.out > names
    
  2. Map those accessions to their taxids.

    while read -r accession; do
       printf '%s\t%s\n' "$accession"  \
           $(curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" | 
              awk -F'[<>]' '/Name="TaxId"/{print $3}'); 
    done < names > taxids
    

    With your example input, that produces:

    $ while read -r accession; do
    >    printf '%s\t%s\n' "$accession"  \
    >        $(curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" | 
    >           awk -F'[<>]' '/Name="TaxId"/{print $3}'); 
    > done < names
    MN464146.1   1435485
    MN148428.1   2661902
    MK671313.1   1239566
    MK671312.1   1239566
    MK671311.1   1239566
    

    Now, this will probably be way too slow, doing them one at a time. So if you have many IDs, try and combine multiple queries into one. Like this:

    perl -ne 'chomp; print $. % 50 == 0 ? "$_\n" : "$_," ' names | 
     | sed 's/,$/\n/' | 
        while read accession; do 
          curl -s "https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&id=${accession}&rettype=summary&retmode=text" | 
            awk -F'[<>]' '/Name="Caption"/{ acc=$3 }
                       /Name="TaxId"/{printf "%s\t%s\n", acc, $3}'
        done  < names > taxids
    

    Here, the first Perl command will read through the list of accessions and print the accession followed by a comma and, once every 50 lines ($. % 50 == 0), it will instead print it with a newline. The sed command will replace a trailing comma with a newline for the very last line of output where, unless the number of accessions is an exact multiple of 50, you will have a line ending with a comma instead of a newline.

    The result is a list of 50 comma-separated identifiers per line (fewer, on the last line). This is then fed to the curl command and then a naive awk that parses it and prints out the accession and the taxid.

    Important: this approach loses the version information. So MK671311.1 becomes MK671311. This shouldn't be a problem, and we will take it into account later, but it should be understood. This is what the output looks like:

    MN464146    1435485
    MN148428    2661902
    MK671313    1239566
    MK671312    1239566
    MK671311    1239566
    
  3. Remove the SARS-CoV-2 sequences.

    The taxid for SARS-CoV-2 is 2697049, so you can get the list of its accessions with:

    awk -F'\t' '$2=="2697049"{print $1}' taxids 
    

    You can now filter your file to remove them using:

    awk '(NR==FNR){ 
           if($2=="2697049"){
             sars[$1]
           } 
           next
          }
          { 
            acc=$2; 
            sub(/\.[0-9]*$/, "", acc); 
            if(acc in sars){ next } 
              print
          }' taxids blast.out > blast.filtered.out
    

    This awk script will read the first input file (NR is the current input line number while FNR the line number of the file being read; the two are only equal while reading the first file) and if its second field is the desired taxid, it saves the first field, the accession, in the array sars. The next ensures we don't do anything else while reading the first file.

    When we move to the second file, we make a copy of the second field, the accession, and save it as acc, and then we remove the version information ( sub(/\.[0-9]*$/, "", acc). If the resulting accession is found in the sars array, we skip the line with next and if not, we print.

    This assumes you are using the faster, batch retrieval I gave above, the one that removes the versioning information. If you instead use the very slow first one, you can filter with

     awk '(NR==FNR){ if($2=="2697049"){ sars[$1]} next }
          { if($2 in sars){next} print}' taxids blast.out > blast.filtered.out
    

That should give you a filtered blast output file that no loner has sequences from taxid 2697049. However, please make sure to test this, count the number of input and output lines and make sure everything makes sense. I don't have access to your actual data and while this does work with the example you provided, I cannot guarantee it will work out of the box for your real data. As always, you need to be careful!

$\endgroup$
1
  • $\begingroup$ @FarzadBeikpour kindly consider upvoting. This is a good answer and a lot of work. If you can't get -negative_taxids working I'd use this approach $\endgroup$
    – M__
    Nov 5, 2023 at 2:16
0
$\begingroup$

I would personally look at is the version number of the blast you are using.

Bugs around blast+ versions >= 2.9.0 (April 2019) for -taxids and I would assume -negative_taxids flag should have been fixed as described in the the Blast+ release notes here.

Fix taxID filtering combined with mask-based alias BLAST databases.

Bug fix descriptions are usually short, hence the lack of further detail.

Blast taxids can be implemented via lists here https://www.ncbi.nlm.nih.gov/books/NBK569846/

$\endgroup$

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.