Answers.Body: <p>Using the <code>FastaToTbl</code> and <code>TblToFasta</code> scripts I have <a href="https://bioinformatics.stackexchange.com/a/2651/298">posted before</a>, you can do:</p>
  
  <pre class="lang-none prettyprint-override"><code>FastaToTbl file.fa | grep -vwf ids.txt | TblToFasta &gt; file.2.fa
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-20T19:02:00.587Z
Answers.LastActivityDate: 2018-03-23T08:58:29.123Z
Answers.LastEditDate: 2018-03-23T08:58:29.123Z
Answers.UserDisplayName: terdon
Answers.Body: <p>Just today I wrote a script to do exactly this using <code>Biopython</code>. I also add a warning If any the headers I wanted to filter was not present in the original fasta. So the python script <code>filter_fasta_by_list_of_headers.py</code> is :</p>
  <pre><code>#!/usr/bin/env python3
  
  from Bio import SeqIO
  import sys
  
  ffile = SeqIO.parse(sys.argv[1], &quot;fasta&quot;)
  header_set = set(line.strip() for line in open(sys.argv[2]))
  
  for seq_record in ffile:
      try:
          header_set.remove(seq_record.name)
      except KeyError:
          print(seq_record.format(&quot;fasta&quot;))
          continue
  if len(header_set) != 0:
      print(len(header_set),'of the headers from list were not identified in the input fasta file.', file=sys.stderr)
  </code></pre>
  <p>and usage :</p>
  <pre><code>filter_fasta_by_list_of_headers.py input.fasta list_of_scf_to_filter &gt; filtered.fasta
  </code></pre>
  <p>P.S. it's quite easy to turn over the script to extract the sequences from the list (just the print line would have to move after line <code>header_set.remove(seq_record.name</code>)</p>
  <p>I made a modification to use <code>in</code> and benchmarked it on normal insect genome assembly - 1Gbp, 300,000 scaffolds and I gave it lists of 5,000 and 100,000 scaffolds to filter. Regardless of the size of the list or algorithm used it took 29-30s. In other words, if LBYL is faster than EAFP, it's on a scale that totally does not matter.</p>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-23T13:30:14.580Z
Answers.Text: @KamilSJaron Nice!
Answers.UserDisplayName: Konrad Rudolph
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-26T16:08:52.223Z
Answers.Text: Ah +1 for showing my intuitions to be wrong, one other thought, it *should* be faster to write all the records at once rather than one by one, is something like this faster? https://pastebin.com/n7RUgW5x Sorry if it seems like micro-optimization, I'm just curious
Answers.UserDisplayName: Chris_Rands
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-03-21T12:28:47.583Z
Answers.LastActivityDate: 2024-01-04T21:22:20.490Z
Answers.LastEditDate: 2024-01-04T21:22:20.490Z
Answers.UserDisplayName: Kamil S Jaron
Answers.Body: <p>Suppose you keep sequence names in <code>ids.txt</code> and sequences in <code>seq.fa</code>:</p>
  
  <pre><code>awk 'BEGIN{while((getline&lt;"ids.txt")&gt;0)l["&gt;"$1]=1}/^&gt;/{f=!l[$1]}f' seq.fa
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-21T02:08:59.863Z
Answers.LastActivityDate: 2018-03-21T02:08:59.863Z
Answers.UserDisplayName: user172818
Answers.Body: <p>If you want to learn how to do things with command-line tools, you can linearize the FASTA with <code>awk</code>, pipe to <code>grep</code> to filter for items of interest named in <code>patterns.txt</code>, then pipe to <code>tr</code> to delinearize:</p>
  
  <pre><code>$ awk '{ if ((NR&gt;1)&amp;&amp;($0~/^&gt;/)) { printf("\n%s", $0); } else if (NR==1) { printf("%s", $0); } else { printf("\t%s", $0); } }' in.fa | grep -Ff patterns.txt - | tr "\t" "\n" &gt; out.fa
  </code></pre>
  
  <p>This will work on multiline FASTA.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-16T04:52:38.637Z
Answers.Text: Excellent solution, except it will keep the sequence specified in patterns.txt, not remove them
Answers.UserDisplayName: Federico Giorgi
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-04-16T17:42:04.783Z
Answers.Text: Good point. I must have missed that part of the question. One can use `grep -v` to remove elements that match specified patterns.
Answers.UserDisplayName: Alex Reynolds
Answers.ContentLicense: CC BY-SA 3.0
Answers.CreationDate: 2018-03-20T20:03:21.150Z
Answers.LastActivityDate: 2018-03-20T21:49:46.863Z
Answers.LastEditDate: 2018-03-20T21:49:46.863Z
Answers.UserDisplayName: Alex Reynolds
Answers.Body: <p>This is what <strong><code>samtools faidx</code></strong> is intended for.</p>
  <p>It needs to be called twice.</p>
  <ol>
  <li><p>First, to create a FASTA index (<code>*.fai</code> file):</p>
  <pre class="lang-sh prettyprint-override"><code>samtools faidx input.fasta
  </code></pre>
  </li>
  <li><p>Secondly, to filter the FASTA file with the help of the index:</p>
  <pre class="lang-sh prettyprint-override"><code>samtools faidx -o output.fasta input.fasta ids…
  </code></pre>
  </li>
  </ol>
  <p><code>ids…</code> are the IDs to retain, as individual arguments separated by space. If you’ve got a file <code>blocklist.txt</code> with IDs you want to <em>discard</em> (one per line), you first need to invert this, after having created the index (using Bash syntax):<sup>1</sup></p>
  <pre class="lang-bash prettyprint-override"><code>remove_ids=($(awk '{print $1}' input.fasta.fai | grep -v -f blocklist.txt))
  </code></pre>
  <p>… and then we can invoke the command as</p>
  <pre><code>samtools faidx -o output.fasta input.fasta &quot;${remove_ids[@]}&quot;
  </code></pre>
  <p>What’s nice about <code>samtools faidx</code> is that it also works with <a href="http://www.htslib.org/doc/bgzip.html" rel="nofollow noreferrer">BGZF-compressed FASTA</a> data — i.e. <code>fa.gz</code> files (but <em>only</em> if they were compressed using <code>bgzip</code>, which ships with Samtools, and the index must be generated on the compressed file). However, it will always write FASTA. If you want gzipped FASTA output, you need to manually stream the result into <code>bgzip</code>:</p>
  <pre><code>
  samtools faidx input.fasta.gz &quot;${remove_ids[@]}&quot; | bgzip -c &gt;output.fasta.gz
  </code></pre>
  <hr />
  <p><sup>1</sup> <sub><strong>Careful</strong>, this syntax is <em>generally unsafe</em> for reading a command output into a variable; we can use it here only because we assume that we want to generate a <em>space</em>-separated list of identifiers; the generally safe way is <a href="https://unix.stackexchange.com/a/485262/3651">a lot more convoluted</a>.</sub></p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2020-09-24T15:27:05.260Z
Answers.LastActivityDate: 2020-09-24T15:27:05.260Z
Answers.UserDisplayName: Konrad Rudolph
Answers.Body: <p>Firstly, you have to install <code>biopython</code> using the cmd <code>pip install biopython</code>.</p>
  <p>Suppose you have two files <code>id.txt</code> and <code>sequence.fasta</code>. The solution file name is <code>seq_rm.py</code> which contains the following codes.</p>
  <pre><code>import argparse
  from Bio import SeqIO
  from Bio.SeqRecord import SeqRecord
  
  
  if __name__ == '__main__':
      parser = argparse.ArgumentParser()
      parser.add_argument('--id_file', required=True,
                          help='file path of removable ids')
      parser.add_argument('--input', required=True,
                          help='path of fasta sequence')
      parser.add_argument('--output', required=True,
                          help='path of selected fasta sequence')
  
      args = parser.parse_args()
  
      with open(args.id_file, mode='r', encoding='utf-8') as file:
          ids = set(file.read().splitlines())
  
      records = [SeqRecord(record.seq, id=record.id, description='', name='') for record in SeqIO.parse(args.input, 'fasta')
                 if record.id not in ids]
      SeqIO.write(records, args.output, &quot;fasta&quot;)
  
  </code></pre>
  <p>Run <code>cmd</code>:</p>
  <pre><code>python3 seq_rm.py --id_file id.txt --input sequence.fasta --output selected.fasta
  </code></pre>
  <p>It will generate an output file named <code>selected.fasta</code>.</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-13T05:00:49.723Z
Answers.LastActivityDate: 2021-09-18T08:14:24.903Z
Answers.LastEditDate: 2021-09-18T08:14:24.903Z
Answers.UserDisplayName: Rajan Raju
Answers.Body: <p>One way using <a href="/questions/tagged/awk" class="post-tag" title="show questions tagged &#39;awk&#39;" rel="tag">awk</a> without <a href="https://awk.freeshell.org/AllAboutGetline" rel="nofollow noreferrer">getline</a>:</p>
  <pre><code>awk 'FNR==NR { a[$0]; next } /^&gt;/ { f = !(substr($0, 2) in a) }f' ids.txt seqs.fa
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-09-14T12:33:55.353Z
Answers.LastActivityDate: 2021-09-14T12:33:55.353Z
Answers.UserDisplayName: Steve
Answers.Body: <p>I know you do not want a Python solution but nevertheless I would recommend it like others. Here is one way to do it using Biopython:</p>
  
  <pre class="lang-python prettyprint-override"><code>from Bio import SeqIO
  
  ids_file = 'ids.txt'
  input_f  = 'input_file.fasta'
  output_f = 'result.fasta'
  
  my_dictionary = {}
  with open(ids_file, 'r') as f:
      for line in f:
          my_dictionary[line[:-1]]='value'
  
  with open(input_f, 'r') as input_handle, open(output_f, 'w') as output_handle:
      for seq_record in SeqIO.parse(input_handle, &quot;fasta&quot;):
          if not seq_record.id in my_dictionary:
              SeqIO.write(seq_record, output_handle, &quot;fasta&quot;)
  </code></pre>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2024-01-05T01:22:45.260Z
Answers.LastActivityDate: 2024-01-05T01:22:45.260Z
Answers.UserDisplayName: Supertech
Answers.Body: <p>Shortest way...</p>
  <pre><code>awk '!/Id2|Id3/' RS=&quot;&gt;&quot; ORS=&quot;&gt;&quot; fasta.fa
  </code></pre>
  <p><code>!</code>= Not match the pattern
  <code>Note:</code> Sequence doesn't need to be in a single line!</p>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2021-12-08T06:46:41.377Z
Answers.LastActivityDate: 2021-12-08T06:46:41.377Z
Answers.UserDisplayName: Aswin S Soman
Body: <p>I have a fasta file like this:</p>
  
  <pre><code>&gt;Id1
  ATCCTT
  &gt;Id2
  ATTTTCCC
  &gt;Id3
  TTTCCCCAAAA
  &gt;Id4
  CCCTTTAAA
  </code></pre>
  
  <p>I want to delete sequences that have the following IDs.</p>
  
  <pre><code>Id2
  Id3
  </code></pre>
  
  <p>The IDs are in a .txt file, and the text file will be used to match and delete those sequences. </p>
  
  <p>My output should be a <code>fasta</code> file like this</p>
  
  <pre><code>&gt;Id1
  ATCCTT
  &gt;Id4
  CCCTTTAAA
  </code></pre>
  
  <p>But I want it with <code>awk</code> and/or <code>sed</code> and/or <code>bash</code> (no python or perl).<br>
  How can I do it?</p>
  
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-20T19:03:32.310Z
Comments.Text: why the language restrictions?
Comments.UserDisplayName: Chris_Rands
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-23T09:11:55.097Z
Comments.Text: Exactly: use the right tool for the job. Unless you can make some very strong assumptions about your input files, use a proper FASTA parser, not some hacked-together cruft.
Comments.UserDisplayName: Konrad Rudolph
Comments.ContentLicense: CC BY-SA 3.0
Comments.CreationDate: 2018-03-26T05:00:55.023Z
Comments.Text: Use whatever tool you want to use if it gets your work done before a deadline. A lot of parsers are slow as mud as compared with Unix tools, but I admit they take a lot of drudgery out of understanding how things work.
Comments.UserDisplayName: Alex Reynolds
Tags: <fasta><awk><shell><benchmarking>
Title: Remove/delete sequences by ID from multifasta
ThreadText: # Remove/delete sequences by ID from multifasta
  
  (Asked by: username_0 on 2018-03-20)
  
  I have a fasta file like this:
  
  ```
  >Id1
  ATCCTT
  >Id2
  ATTTTCCC
  >Id3
  TTTCCCCAAAA
  >Id4
  CCCTTTAAA
  
  ```
  
  I want to delete sequences that have the following IDs.
  
  ```
  Id2
  Id3
  
  ```
  
  The IDs are in a .txt file, and the text file will be used to match and delete those sequences.
  
  My output should be a `fasta` file like this
  
  ```
  >Id1
  ATCCTT
  >Id4
  CCCTTTAAA
  
  ```
  
  But I want it with `awk` and/or `sed` and/or `bash` (no python or perl).  
  How can I do it?
  
  **Tags:** fasta, awk, shell, benchmarking
  
  > **Comments:**
  > **username_1**: why the language restrictions?  
  > **username_2**: Exactly: use the right tool for the job. Unless you can make some very strong assumptions about your input files, use a proper FASTA parser, not some hacked-together cruft.  
  > **username_3**: Use whatever tool you want to use if it gets your work done before a deadline. A lot of parsers are slow as mud as compared with Unix tools, but I admit they take a lot of drudgery out of understanding how things work.  
  
  ---
  
  ## 9 answers
  
  ### Answer by username_5 on 2018-03-21. Score: 7
  
  Just today I wrote a script to do exactly this using `Biopython`. I also add a warning If any the headers I wanted to filter was not present in the original fasta. So the python script `filter_fasta_by_list_of_headers.py` is :
  
  ```
  #!/usr/bin/env python3
  
  from Bio import SeqIO
  import sys
  
  ffile = SeqIO.parse(sys.argv[1], "fasta")
  header_set = set(line.strip() for line in open(sys.argv[2]))
  
  for seq_record in ffile:
      try:
          header_set.remove(seq_record.name)
      except KeyError:
          print(seq_record.format("fasta"))
          continue
  if len(header_set) != 0:
      print(len(header_set),'of the headers from list were not identified in the input fasta file.', file=sys.stderr)
  
  ```
  
  and usage :
  
  ```
  filter_fasta_by_list_of_headers.py input.fasta list_of_scf_to_filter > filtered.fasta
  
  ```
  
  P.S. it's quite easy to turn over the script to extract the sequences from the list (just the print line would have to move after line `header_set.remove(seq_record.name`)
  
  I made a modification to use `in` and benchmarked it on normal insect genome assembly - 1Gbp, 300,000 scaffolds and I gave it lists of 5,000 and 100,000 scaffolds to filter. Regardless of the size of the list or algorithm used it took 29-30s. In other words, if LBYL is faster than EAFP, it's on a scale that totally does not matter.
  
  > **Comments:**
  > **username_2**: @username_5 Nice!  
  > **username_1**: Ah +1 for showing my intuitions to be wrong, one other thought, it *should* be faster to write all the records at once rather than one by one, is something like this faster? https://pastebin.com/n7RUgW5x Sorry if it seems like micro-optimization, I'm just curious  
  
  ---
  
  ### Answer by username_6 on 2018-03-21. Score: 7
  
  Suppose you keep sequence names in `ids.txt` and sequences in `seq.fa`:
  
  ```
  awk 'BEGIN{while((getline<"ids.txt")>0)l[">"$1]=1}/^>/{f=!l[$1]}f' seq.fa
  
  ```
  
  
  ---
  
  ### Answer by username_3 on 2018-03-20. Score: 6
  
  If you want to learn how to do things with command-line tools, you can linearize the FASTA with `awk`, pipe to `grep` to filter for items of interest named in `patterns.txt`, then pipe to `tr` to delinearize:
  
  ```
  $ awk '{ if ((NR>1)&&($0~/^>/)) { printf("\n%s", $0); } else if (NR==1) { printf("%s", $0); } else { printf("\t%s", $0); } }' in.fa | grep -Ff patterns.txt - | tr "\t" "\n" > out.fa
  
  ```
  
  This will work on multiline FASTA.
  
  > **Comments:**
  > **username_7**: Excellent solution, except it will keep the sequence specified in patterns.txt, not remove them  
  > **username_3**: Good point. I must have missed that part of the question. One can use `grep -v` to remove elements that match specified patterns.  
  
  ---
  
  ### Answer by username_4 on 2018-03-20. Score: 3
  
  Using the `FastaToTbl` and `TblToFasta` scripts I have [posted before](https://bioinformatics.stackexchange.com/a/2651/298), you can do:
  
  ```none
  FastaToTbl file.fa | grep -vwf ids.txt | TblToFasta > file.2.fa
  
  ```
  
  
  ---
  
  ### Answer by username_2 on 2020-09-24. Score: 2
  
  This is what **`samtools faidx`** is intended for.
  
  It needs to be called twice.
  
  1. First, to create a FASTA index (`*.fai` file):
  
     ```sh
     samtools faidx input.fasta
  
     ```
  2. Secondly, to filter the FASTA file with the help of the index:
  
     ```sh
     samtools faidx -o output.fasta input.fasta ids…
  
     ```
  
  `ids…` are the IDs to retain, as individual arguments separated by space. If you’ve got a file `blocklist.txt` with IDs you want to *discard* (one per line), you first need to invert this, after having created the index (using Bash syntax):1
  
  ```bash
  remove_ids=($(awk '{print $1}' input.fasta.fai | grep -v -f blocklist.txt))
  
  ```
  
  … and then we can invoke the command as
  
  ```
  samtools faidx -o output.fasta input.fasta "${remove_ids[@]}"
  
  ```
  
  What’s nice about `samtools faidx` is that it also works with [BGZF-compressed FASTA](http://www.htslib.org/doc/bgzip.html) data — i.e. `fa.gz` files (but *only* if they were compressed using `bgzip`, which ships with Samtools, and the index must be generated on the compressed file). However, it will always write FASTA. If you want gzipped FASTA output, you need to manually stream the result into `bgzip`:
  
  ```
  
  samtools faidx input.fasta.gz "${remove_ids[@]}" | bgzip -c >output.fasta.gz
  
  ```
  
  ---
  
  1 **Careful**, this syntax is *generally unsafe* for reading a command output into a variable; we can use it here only because we assume that we want to generate a *space*-separated list of identifiers; the generally safe way is [a lot more convoluted](https://unix.stackexchange.com/a/485262/3651).
  
  
  ---
  
  ### Answer by username_11 on 2021-12-08. Score: 1
  
  Shortest way...
  
  ```
  awk '!/Id2|Id3/' RS=">" ORS=">" fasta.fa
  
  ```
  
  `!`= Not match the pattern
  `Note:` Sequence doesn't need to be in a single line!
  
  
  ---
  
  ### Answer by username_8 on 2021-09-13. Score: 0
  
  Firstly, you have to install `biopython` using the cmd `pip install biopython`.
  
  Suppose you have two files `id.txt` and `sequence.fasta`. The solution file name is `seq_rm.py` which contains the following codes.
  
  ```
  import argparse
  from Bio import SeqIO
  from Bio.SeqRecord import SeqRecord
  
  if __name__ == '__main__':
      parser = argparse.ArgumentParser()
      parser.add_argument('--id_file', required=True,
                          help='file path of removable ids')
      parser.add_argument('--input', required=True,
                          help='path of fasta sequence')
      parser.add_argument('--output', required=True,
                          help='path of selected fasta sequence')
  
      args = parser.parse_args()
  
      with open(args.id_file, mode='r', encoding='utf-8') as file:
          ids = set(file.read().splitlines())
  
      records = [SeqRecord(record.seq, id=record.id, description='', name='') for record in SeqIO.parse(args.input, 'fasta')
                 if record.id not in ids]
      SeqIO.write(records, args.output, "fasta")
  
  ```
  
  Run `cmd`:
  
  ```
  python3 seq_rm.py --id_file id.txt --input sequence.fasta --output selected.fasta
  
  ```
  
  It will generate an output file named `selected.fasta`.
  
  
  ---
  
  ### Answer by username_9 on 2021-09-14. Score: 0
  
  One way using [awk](/questions/tagged/awk "show questions tagged 'awk'") without [getline](https://awk.freeshell.org/AllAboutGetline):
  
  ```
  awk 'FNR==NR { a[$0]; next } /^>/ { f = !(substr($0, 2) in a) }f' ids.txt seqs.fa
  
  ```
  
  
  ---
  
  ### Answer by username_10 on 2024-01-05. Score: 0
  
  I know you do not want a Python solution but nevertheless I would recommend it like others. Here is one way to do it using Biopython:
  
  ```python
  from Bio import SeqIO
  
  ids_file = 'ids.txt'
  input_f  = 'input_file.fasta'
  output_f = 'result.fasta'
  
  my_dictionary = {}
  with open(ids_file, 'r') as f:
      for line in f:
          my_dictionary[line[:-1]]='value'
  
  with open(input_f, 'r') as input_handle, open(output_f, 'w') as output_handle:
      for seq_record in SeqIO.parse(input_handle, "fasta"):
          if not seq_record.id in my_dictionary:
              SeqIO.write(seq_record, output_handle, "fasta")
  
  ```
  
  
  ---
  
  