How to filter a fasta file?

Question

[oussama.badad]

Dear All

I am trying to remove the sequences with no functional information from a functional annotation fasta file

>Oeu043104.1|---NA---
MIESNFWDACWPHCLLRVLLSLLAESASQPLCPPLQRYNPKYLEDDYGVNQATEWLFYTPRDRENEIENIRNGVAVDGY
>Oeu043107.1|g7lir3_medtr alpha beta-hydrolase superfamily protein os=medicago truncatula gn=mtr_8g086260 pe=4 sv=1
MSTTQGTSPRGNINVKDEPDHLLVLVHGIMGSPSDWTYFEADLKRRLGKRFLIYASSCNTYTKTFTGIDGAGKRLAEEVMEIVRNTESLKKISFLAHSLG
GLFSRYAIAVLYMPNTSSDDSSVIAGSTNTSLKTSCYSNTGLIAGLEPSNFITLATPHLGVRGKKQVNPFSIILIDGPVLPFLLGLPFLEKIAAPLAPIF
TGRTGSQLFLTDGQPDRPPLLLRMASDCKDGKFVSALGAFRCRLLYANVSYDHMVGWRTSSIRRETELIKPPLQSLDGYKHVVSVEYCPPVSSEGPHFPE
EAAKAKQAAQNEPNNQNTVEYHETMEEEMIRGLQRLGWKKVDVSFHSAFWPFFAHNNINVKNEWLYNAGVGVVAHVADNIKQQENQQGSTYVAASL

I was wondering if someone can help me with a python script or a shell script

Thank you
Oussama

[Sukhi Singh]

What have you tried?

Take some thoughts from these posts:

• Filtering Fasta Sequences By Chromosomes Names From A Big Fasta File
• How To Remove Certain Sequences From A Fasta File
• Removing sequences from .fasta files based on name prefix
• http://qiime.org/scripts/filter_fasta.html

[GouthamAtla]

You could convert to single line fasta and use a simple grep to remove lines with --NA-- or use bio Python to have more control.

Answer 1

[st.ph.n]

Since you mentioned you wanted python:

Edit 'parse_fasta.py':

#!/usr/bin/env python

import sys

input_file = sys.argv[1]
output file = sys.argv[2]

with open(input_file, 'r') as f:
    headers = []
    seqs = []

    for line in f:
        if line.startswith(">"):
            headers.append(line.strip())
        else:
            seqs.append(line.strip())

myseqs = dict(zip(headers,seqs))

with open(outfile, 'w') as out: 
    for m in myseqs:
        if '---NA---' not in m:
            print >> out, m, '\n', myseqs[m]

Usage:

python parse_fasta.py input.fasta output.fasta

Answer 2

[novice]

Too easy bro:

#!/usr/bin/perl

use strict; use warnings;

my $print;

while (<>) {
    $print = m/---NA---/ ? 0: 1 if m/>/;
    print if $print;
}

Usage: perl filter.pl sequences.fasta > filtered.fasta or ./filter.pl ... if you made the script executable.

Since I've got time on my hands and like Perl, I wrote you another script. This one would "slurp" the entire file into memory, so it's probably best to avoid if your file is huge.

#!/usr/bin/perl

use strict; use warnings;

my $whole = do { local $/; <> };
my @keep = map { s/\A([^>])/>$1/; $_ }
    grep { ! m/---NA---/; } split />/, $whole;
print @keep;