<p>A Biopytnon solution that works for <code>.gz</code> files too :-). The script make a script <code>filter_faa.py</code>:</p>

<pre><code>#!/usr/bin/env python3

import sys
from Bio import SeqIO
from mimetypes import guess_type
import gzip

filename = sys.argv[1]

encoding = guess_type(filename)[1]
_open = partial(gzip.open, mode='rt') if encoding == 'gzip' else open

with _open(filename) as ffile:
    for seq_record in SeqIO.parse(ffile, 'fasta'):
        if seq_record.seq.count('X') &lt; 3:
            print(seq_record.format("fasta"))
</code></pre>

<p>Then usage is</p>

<pre><code>python3 filter_faa.py my.faa &gt; my_filtered.faa 
</code></pre>

<p>The <code>my_filtered.faa</code> then will contain</p>

<pre><code>&gt;Prot1
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF

&gt;Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF
</code></pre>

<p>The filtering condition is hard-coded. You could potentially parametrise it if you feel like later you might want to filter by a different number of unknown aa.</p>

CC BY-SA 4.0
2020-05-07T07:19:27.987Z
it's generally better to explicitly close the open file handle, which won't happen for your gzip.open() case, could be modified, e.g. https://stackoverflow.com/a/52839332
Chris_Rands
CC BY-SA 4.0
2020-05-07T12:15:26.887Z
Ha, I always wondered if there is an elegant soltion for opening it with/out gz. Cheers Chris!
Kamil S Jaron
CC BY-SA 4.0
2020-05-07T13:25:53.970Z
It is an elegant solution even without closing the filehandle
M__
CC BY-SA 4.0
2020-05-06T13:38:07.007Z
2020-05-07T12:26:49.433Z
2020-05-07T12:26:49.433Z
Kamil S Jaron
<p><strong>Updated</strong>: as pointed out, this will not work for a multifasta as asked for in your original question. However it would work fine for a single line fasta</p>

<hr>

<p>You could accomplish this with a single line of bash code if you would like (this assumes your fasta is in gzip format, if not just use cat):</p>

<pre class="lang-sh prettyprint-override"><code>zcat prot.fasta.gz | grep -E -v -B1 '(^&gt;.*|^(.*X.*){3}$)' | grep -E -v '^--$'  - | gzip &gt; prot.filtered.fasta.gzip
</code></pre>

<p>The <code>-v</code> option of grep inverts the regex, and the <code>-B1</code> will print out the line occuring before the current context. The <code>-B1</code> introduces <code>--</code> between the matches, which is why the second grep is there. Then you can immediately pipe it back into a zipped file </p>

CC BY-SA 4.0
2020-05-07T13:10:07.143Z
I guess this would not work if the aa sequence had more than one line, or am I missing something?
Kamil S Jaron
CC BY-SA 4.0
2020-05-07T13:20:49.537Z
You would be correct, In which case this would not work with OP's original question.  Grep would still be a useful tool for the filtering. You could probably filter out any sequences which havae More then 3 x's then filter out any Protein s that no longer have sequences
Patrick Magee
CC BY-SA 4.0
2020-05-07T12:58:36.777Z
2020-05-08T12:03:30.973Z
2020-05-08T12:03:30.973Z
Patrick Magee
<p>I have 5000 protein sequences in one multifasta file. I found more reads have gaps as X in their reads. So, want to eliminate those reads completely (Whole protein seq) from the file. I am keeping filter criteria as if a read contains morethan 2 X ( continuesly or anywhere in the read) should be removed. Thanks in advance for your help. </p>

<p>The input sequence looks like this </p>

<pre><code>&gt;Prot1 
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF
&gt;Prot2
ANSTVKKKKLLLYYYSSSEERXXXXXXXXXXXFGHYFGHYFGHFYVHFGFYVHCEDYHF
&gt;Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF
</code></pre>

<p>I want output Like this </p>

<pre><code>&gt;Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF
</code></pre>

CC BY-SA 4.0
2020-05-06T12:10:05.847Z
Easy use Biopython "SequenceClean".You could do it as a Perl one liner or awk. Which of these solutions are you more familiar with?
M__
CC BY-SA 4.0
2020-05-06T13:39:13.473Z
Just for a record, X is not a gap but unknown amino acid.
Kamil S Jaron
CC BY-SA 4.0
2020-05-06T13:42:38.110Z
Also, `>Prot1 ` seems to have only one X, hence I suppose it should be kept as well, is that correct?
Kamil S Jaron
<phylogenetics><ngs><parsing>
How to remove sequence reads contains more than 2 X from multifasta file?
# How to remove sequence reads contains more than 2 X from multifasta file?

(Asked by: username_0 on 2020-05-06)

I have 5000 protein sequences in one multifasta file. I found more reads have gaps as X in their reads. So, want to eliminate those reads completely (Whole protein seq) from the file. I am keeping filter criteria as if a read contains morethan 2 X ( continuesly or anywhere in the read) should be removed. Thanks in advance for your help. 
The input sequence looks like this 
>Prot1 
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF
>Prot2
ANSTVKKKKLLLYYYSSSEERXXXXXXXXXXXFGHYFGHYFGHFYVHFGFYVHCEDYHF
>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

I want output Like this 
>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

**Tags:** phylogenetics, ngs, parsing
> **Comments:**
> **username_1**: Easy use Biopython "SequenceClean".You could do it as a Perl one liner or awk. Which of these solutions are you more familiar with?  
> **username_2**: Just for a record, X is not a gap but unknown amino acid.  
> **username_2**: Also, `>Prot1 ` seems to have only one X, hence I suppose it should be kept as well, is that correct?  

---

## 2 answers

### Answer by username_2 on 2020-05-06. Score: 1

A Biopytnon solution that works for .gz files too :-). The script make a script filter_faa.py:
#!/usr/bin/env python3

import sys
from Bio import SeqIO
from mimetypes import guess_type
import gzip

filename = sys.argv[1]

encoding = guess_type(filename)[1]
_open = partial(gzip.open, mode='rt') if encoding == 'gzip' else open

with _open(filename) as ffile:
    for seq_record in SeqIO.parse(ffile, 'fasta'):
        if seq_record.seq.count('X') < 3:
            print(seq_record.format("fasta"))

Then usage is
python3 filter_faa.py my.faa > my_filtered.faa 

The my_filtered.faa then will contain
>Prot1
ANSTVKKKKLLLYYYSSSEERXFGHYFGHYFGHFYVHFGFYVHCEDYHF

>Prot3
ANSTVKKKKLLLYYYSSSEERFGHYFGHYFGHFYVHFGFYVHCEDYHF

The filtering condition is hard-coded. You could potentially parametrise it if you feel like later you might want to filter by a different number of unknown aa.

> **Comments:**
> **username_3**: it's generally better to explicitly close the open file handle, which won't happen for your gzip.open() case, could be modified, e.g. https://stackoverflow.com/a/52839332  
> **username_2**: Ha, I always wondered if there is an elegant soltion for opening it with/out gz. Cheers Chris!  
> **username_1**: It is an elegant solution even without closing the filehandle  

---

### Answer by username_4 on 2020-05-07. Score: 0

Updated: as pointed out, this will not work for a multifasta as asked for in your original question. However it would work fine for a single line fasta

You could accomplish this with a single line of bash code if you would like (this assumes your fasta is in gzip format, if not just use cat):
zcat prot.fasta.gz | grep -E -v -B1 '(^>.*|^(.*X.*){3}$)' | grep -E -v '^--$'  - | gzip > prot.filtered.fasta.gzip

The -v option of grep inverts the regex, and the -B1 will print out the line occuring before the current context. The -B1 introduces -- between the matches, which is why the second grep is there. Then you can immediately pipe it back into a zipped file

> **Comments:**
> **username_2**: I guess this would not work if the aa sequence had more than one line, or am I missing something?  
> **username_4**: You would be correct, In which case this would not work with OP's original question.  Grep would still be a useful tool for the filtering. You could probably filter out any sequences which havae More then 3 x's then filter out any Protein s that no longer have sequences  

---

