# How to remove sequence reads contains more than 2 X from multifasta file?

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

• 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__
May 6 '20 at 12:10
• Just for a record, X is not a gap but unknown amino acid. May 6 '20 at 13:39
• Also, >Prot1  seems to have only one X, hence I suppose it should be kept as well, is that correct? May 6 '20 at 13:42

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.

• 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. stackoverflow.com/a/52839332 May 7 '20 at 7:19
• Ha, I always wondered if there is an elegant soltion for opening it with/out gz. Cheers Chris! May 7 '20 at 12:15
• It is an elegant solution even without closing the filehandle
– M__
May 7 '20 at 13:25

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

• I guess this would not work if the aa sequence had more than one line, or am I missing something? May 7 '20 at 13:10
• 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 May 7 '20 at 13:20