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, 2020 at 12:10
  • $\begingroup$ Just for a record, X is not a gap but unknown amino acid. $\endgroup$ May 6, 2020 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$ May 6, 2020 at 13:42

2 Answers 2

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

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

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