<p>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:</p>
<ol>
<li><p>Extract the accessions into a new file</p>
<pre><code>awk '!/^#/{print $2}' blast.out &gt; names
</code></pre>
</li>
<li><p>Map those accessions to their taxids.</p>
<pre class="lang-bash prettyprint-override"><code>while read -r accession; do
   printf '%s\t%s\n' &quot;$accession&quot;  \
       $(curl -s &quot;https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&amp;id=${accession}&amp;rettype=summary&amp;retmode=text&quot; | 
          awk -F'[&lt;&gt;]' '/Name=&quot;TaxId&quot;/{print $3}'); 
done &lt; names &gt; taxids
</code></pre>
<p>With your example input, that produces:</p>
<pre class="lang-bash prettyprint-override"><code>$ while read -r accession; do
&gt;    printf '%s\t%s\n' &quot;$accession&quot;  \
&gt;        $(curl -s &quot;https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&amp;id=${accession}&amp;rettype=summary&amp;retmode=text&quot; | 
&gt;           awk -F'[&lt;&gt;]' '/Name=&quot;TaxId&quot;/{print $3}'); 
&gt; done &lt; names
MN464146.1   1435485
MN148428.1   2661902
MK671313.1   1239566
MK671312.1   1239566
MK671311.1   1239566
</code></pre>
<p>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:</p>
<pre class="lang-bash prettyprint-override"><code>perl -ne 'chomp; print $. % 50 == 0 ? &quot;$_\n&quot; : &quot;$_,&quot; ' names | 
 | sed 's/,$/\n/' | 
    while read accession; do 
      curl -s &quot;https://eutils.ncbi.nlm.nih.gov/entrez/eutils/esummary.fcgi?db=nuccore&amp;id=${accession}&amp;rettype=summary&amp;retmode=text&quot; | 
        awk -F'[&lt;&gt;]' '/Name=&quot;Caption&quot;/{ acc=$3 }
                   /Name=&quot;TaxId&quot;/{printf &quot;%s\t%s\n&quot;, acc, $3}'
    done  &lt; names &gt; taxids
</code></pre>
<p>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 (<code>$. % 50 == 0</code>), it will instead print it with a newline. The <code>sed</code> 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.</p>
<p>The result is a list of 50 comma-separated identifiers per line (fewer, on the last line). This is then fed to the <code>curl</code> command and then a naive awk that parses it and prints out the accession and the taxid.</p>
<p><strong>Important</strong>: this approach loses the version information. So <code>MK671311.1</code> becomes <code>MK671311</code>. 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:</p>
<pre><code>MN464146    1435485
MN148428    2661902
MK671313    1239566
MK671312    1239566
MK671311    1239566
</code></pre>
</li>
<li><p>Remove the SARS-CoV-2 sequences.</p>
<p>The taxid for SARS-CoV-2 is <a href="https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=2697049" rel="nofollow noreferrer">2697049</a>, so you can get the list of its accessions with:</p>
<pre><code>awk -F'\t' '$2==&quot;2697049&quot;{print $1}' taxids 
</code></pre>
<p>You can now filter your file to remove them using:</p>
<pre><code>awk '(NR==FNR){ 
       if($2==&quot;2697049&quot;){
         sars[$1]
       } 
       next
      }
      { 
        acc=$2; 
        sub(/\.[0-9]*$/, &quot;&quot;, acc); 
        if(acc in sars){ next } 
          print
      }' taxids blast.out &gt; blast.filtered.out
</code></pre>
<p>This awk script will read the first input file (<code>NR</code> is the current input line number while <code>FNR</code> 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 <code>sars</code>. The <code>next</code> ensures we don't do anything else while reading the first file.</p>
<p>When we move to the second file, we make a copy of the second field, the accession, and save it as <code>acc</code>, and then we remove the version information (<code> sub(/\.[0-9]*$/, &quot;&quot;, acc)</code>. If the resulting accession is found in the <code>sars</code> array, we skip the line with <code>next</code> and if not, we print.</p>
<p>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</p>
<pre><code> awk '(NR==FNR){ if($2==&quot;2697049&quot;){ sars[$1]} next }
      { if($2 in sars){next} print}' taxids blast.out &gt; blast.filtered.out
</code></pre>
</li>
</ol>
<p>That should give you a filtered blast output file that no loner has sequences from taxid 2697049. <strong>However</strong>, 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!</p>

CC BY-SA 4.0
2023-11-05T02:16:21.277Z
@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
M__
CC BY-SA 4.0
2023-11-04T12:59:15.487Z
2023-11-04T18:38:31.503Z
2023-11-04T18:38:31.503Z
terdon
<p>I would personally look at is the version number of the blast you are using.</p>
<p>Bugs around blast+ versions &gt;= 2.9.0 (April 2019) for <code>-taxids</code> and I would assume <code>-negative_taxids</code> flag should have been fixed as described in the the Blast+ release notes <a href="https://www.ncbi.nlm.nih.gov/books/NBK131777/" rel="nofollow noreferrer">here</a>.</p>
<blockquote>
<p>Fix taxID filtering combined with mask-based alias BLAST databases.</p>
</blockquote>
<p>Bug fix descriptions are usually short, hence the lack of further detail.</p>
<p>Blast <code>taxids</code> can be implemented via lists here <a href="https://www.ncbi.nlm.nih.gov/books/NBK569846/" rel="nofollow noreferrer">https://www.ncbi.nlm.nih.gov/books/NBK569846/</a></p>

CC BY-SA 4.0
2023-11-05T03:10:32.087Z
2023-11-05T14:33:18.217Z
2023-11-05T14:33:18.217Z
M__
<p>I would like to perform blast+, using the command line, here is the code:</p>
<pre><code>-task blastn -db nr -outfmt 7 -query myQuery.fas -out result.txt -remote.
</code></pre>
<p>I have to exclude SARS-CoV-2 hits from my results, so I must include <code>-negative_taxids</code> mentioning SARS-CoV-2 NCBI taxid number. However, the problem is, <code>-negative_taxids</code> is NOT compatible with <code>-remote</code> command and my command line encounters an error.
I have to use <code>-remote</code> 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?</p>
<p>My output txt file for this command line is:</p>
<pre><code># 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
</code></pre>

CC BY-SA 4.0
2023-11-03T10:43:49.380Z
Have you considered running your query against all of `nr` and then removing the SARS hits from the results instead?
terdon
CC BY-SA 4.0
2023-11-03T11:54:09.953Z
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:
Farzad Beikpour
CC BY-SA 4.0
2023-11-03T11:54:21.637Z
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.
Farzad Beikpour
CC BY-SA 4.0
2023-11-03T11:58:17.633Z
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.
terdon
CC BY-SA 4.0
2023-11-05T11:25:23.117Z
What do you want to do with these BLAST hits?
gringer
<blast><sars-cov-2><ncbi><command-line>
How can I limit blast results to exclude a specific taxonomic ID?
# How can I limit blast results to exclude a specific taxonomic ID?

(Asked by: username_0 on 2023-11-03)

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

**Tags:** blast, sars-cov-2, ncbi, command-line
> **Comments:**
> **username_1**: Have you considered running your query against all of `nr` and then removing the SARS hits from the results instead?  
> **username_0**: 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:  
> **username_0**: 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.  
> **username_1**: 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.  
> **username_2**: What do you want to do with these BLAST hits?  

---

## 2 answers

### Answer by username_1 on 2023-11-04. Score: 1

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:

Extract the accessions into a new file
awk '!/^#/{print $2}' blast.out > names


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


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!

> **Comments:**
> **username_3**: @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  

---

### Answer by username_3 on 2023-11-05. Score: 0

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](https://www.ncbi.nlm.nih.gov/books/NBK131777/).

> 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/>


---

