2
$\begingroup$

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
$\endgroup$
3
  • 1
    $\begingroup$ Easy use Biopython "SequenceClean".You could do it as a Perl one liner or awk. Which of these solutions are you more familiar with? $\endgroup$
    – M__
    May 6 '20 at 12:10
  • $\begingroup$ Just for a record, X is not a gap but unknown amino acid. $\endgroup$
    – Kamil S Jaron
    May 6 '20 at 13:39
  • $\begingroup$ Also, >Prot1 seems to have only one X, hence I suppose it should be kept as well, is that correct? $\endgroup$
    – Kamil S Jaron
    May 6 '20 at 13:42
1
$\begingroup$

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.

$\endgroup$
3
  • $\begingroup$ 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 $\endgroup$ May 7 '20 at 7:19
  • 1
    $\begingroup$ Ha, I always wondered if there is an elegant soltion for opening it with/out gz. Cheers Chris! $\endgroup$
    – Kamil S Jaron
    May 7 '20 at 12:15
  • $\begingroup$ It is an elegant solution even without closing the filehandle $\endgroup$
    – M__
    May 7 '20 at 13:25
0
$\begingroup$

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

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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