Extract fasta sequences from a large file using a list of names

Question

[fhsantanna]

I have a large fasta file of 16S sequences and I want to retrieve sequences using a list of organism names. Do you know a script capable of doing it?

EDIT:

Headers look like that:

>S000000859 Bacillus sp. USC14; AF346495
sequence
>S000001027 Paenibacillus borealis; KN25; AJ011325
sequence

And I have a list like the following:

Paenibacilus borealis
Paenibacillus sp. 1-18
Paenibacillus sp. 1-49
Paenibacillus sp. A9
Paenibacillus sp. Aloe-11

I want to retrieve those sequences that match with names present in the list.

[Dan D]

Can you show us some example headers from your FASTA file? Right now I'm thinking to put all the organism names in a text file and simply use grep -f

[fhsantanna]

I believe I was not clear. My 16S fasta file has sequences of hundreds of species. I have to retrieve just dozens of them using a list of organisms of interest.

[marina.v.yurieva]

If the fasta sequence is in one line it's":

grep -A1 'organism name'

If it's not it would be easier to format it into one line and then do grep.

[fhsantanna]

A good way to start. Very simple. Thanks. The problem is when the sequences have different sizes.

Answer 1

[tlawrence3]

You can do this easily by installing the FAST: Fast Analysis of Sequences Toolbox (publication)(github).

You can install it by using this command (only use sudo if you need to):

(sudo) perl -MCPAN -e 'install FAST'

Once it is installed here is a small bash script to do what you need:

#!/bin/bash

while read line ;do
    cat original_fasta.fa | fasgrep -di "$line"  >> reduced_fasta_file.fa
done < species.txt

[venu]

Asking just for clarification. In above bash script which one is the file containing the fasta sequence and which one is the file containing list of organism names?

[tlawrence3]

Updated the example to make it more clear.

[fhsantanna]

Very good. I was able to extract the sequences running the program one by one.

But when I try the bash script, the following error message appears: "fasgrep: expects one regex argument when input is on STDIN."

[venu]

Try ./faSomeRecords very simple one as shown

[fhsantanna]

tlawrence, do you know why I had trouble with your bash script (error message appears: fasgrep: expects one regex argument when input is on STDIN.)?

It seems that fasgrep is not recognizing the $line variable.

[tlawrence3]

It is something in the bash script. I don't know bash scripting well enough to figure it out.

This works on the command line:

fasgrep -di "Paenibacilus borealis" fasta_file.fa

However, it doesn't work in the script. Maybe someone here can help.

[tlawrence3]

I updated the script. It should work now.

[fhsantanna]

It did not work for me.

[tlawrence3]

What didn't work? Did you get an empty fasta file? The wrong sequences? An error? Does doing the command manually give the expected results?

[fhsantanna]

As stated before, fasgrep works fine manually (I accept your answer, as you can see).

However, the bash script did not work (I had to write one without using a loop, but writing a line for every organism of interest). Your bash script did not filtered the organisms of interest, all of the original sequences were included in the output and more...

[tlawrence3]

I guess my bash scripting skills need some improvement. Since I am one of the authors of the toolkit I want to make sure it is working for as many people as possible.

[fhsantanna]

Never mind. Great toolkit. Thanks again.

Just one thing, maybe I found a bug.

Considering a fasta file like the following (note the absence of a underscore between "S000001027" and "Paenibacillus_borealis"):

>S000001027 Paenibacillus_borealis;_KN25;_AJ011325
AAAAAAAAAATTTTTTT

When I type:

fasgrep -di "Paenibacillus_borealis" input.fas

Fasgrep is able to find the sequence of interest.

However, if the fasta sequence has a underscore between S000001027 and Paenibacillus_borealis (>S000001027_Paenibacillus_borealis;_KN25;_AJ011325), fasgrep cannot find the sequence.

[tlawrence3]

It is not a bug. fastA format is loosely defined as:

>identifier<space>description<newline>
<sequence data>

So when you don't have a space in the header line all you have is an identifier which fasgrep searches by default. So instead of:

fasgrep -di

which searches the description field you would use:

fasgrep -i

[Matt Shirley]

+1 for properly handling the fasta defline info field.

Answer 2

[h.mon]

The following script should work:

#!/usr/bin/env perl
my $list_file = $ARGV[0];
my $fasta_in = $ARGV[1];
my $fasta_out = $ARGV[2];
open(LIST_FILE, "<", $list_file) or die "could not open '$list_file' : $! \n";
open(FASTA_IN, "<", $fasta_in) or die "could not open '$fasta_in' : $! \n";
open(FASTA_OUT, ">", $fasta_out) or die "could not open $fasta_out : $! \n";
my @headers = ();
while(<LIST_FILE>) {
    chomp;
    next if ( /^\s*$/ );
    push(@headers, $_.";");
}
my $pat = join '|', map quotemeta, @headers;
$/ = ">";
while(<FASTA_IN>) {
    chomp;
    if ( /$pat/ ) { print FASTA_OUT ">$_"; }
}
close(LIST_FILE);
close(FASTA_IN);
close(FASTA_OUT);

Call it with:

./getFastas.pl list.txt sequences.fasta sequences.out

[fhsantanna]

It did not work. Blank output.

[h.mon]

I tested with the following mock files:

list.txt

Paenibacillus borealis
Paenibacillus sp. 1-18
Paenibacillus sp. 1-49
Paenibacillus sp. A9
Paenibacillus sp. Aloe-11

sequences.fasta:

>S000000859 Bacillus sp. USC14; AF346495
ATGCATGCATGCATGCATGCATGC
>S000001027 Paenibacillus borealis; KN25; AJ011325
GCGCGCATATATGCGCGCATATAT
>S000001026 Paenibacillus borealisensis; KN26; AJ011326
GCGCGCATATATGCGCGCATATAT
>Simaginary Streptomyces borealis; KN25; AJ011325
AAAAAAAAAAAAAAAAAAAAAAAA
>S000000000 Paenibacillus sp. 1-49; nvhjdasnviuas
GGGGGGGGGGGGGGGGGGGGGGGG
>S000001027_Paenibacillus borealis; KN25; AJ011325
AAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAAA

worked all right for me, retrieving three sequences ( "S000001027 Paenibacillus borealis; KN25; AJ011325", "S000001027_Paenibacillus borealis; KN25; AJ011325" and "S000000000 Paenibacillus sp. 1-49; nvhjdasnviuas" ).

[fhsantanna]

For the mock sequences.fasta, it works fine. However, it does not work if the sequences.fasta has the following sequence:

>AB042938.1.1462 Paenibacillus borealis
GACGAACGCUGGCGGCGUGCCUAAUACAUGCAAGUCGAGCGGACUUUGCCUUCGGGUAAAGUUAGCGGCGGACGGGUGAG
UAACACGUGGGUAACCUGCCCAUAAGACUGGGAUAACAUUCGGAAACGAAUGCUAAUACCGGAUACGCGAAUCGGUCGCA
UGAUCGAAUCGGGAAAGGCGGAGCAAUCUGCCACUUAUGGAUGGACCCGCGGCGCAUUAGCUAGUUGGUGGGGUAACGGC
UCACCAAGGCGACGAUGCGUAGCCGACCUGAGAGGGUGAUCGGCCACACUGGGACUGAGACACGGCCCAGACUCCUACGG
GAGGCAGCAGUAGGGAAUCUUCCGCAAUGGACGAAAGUCUGACGGAGCAACGCCGCGUGAGUGAUGAAGGUUUUCGGAUC
GUAAAGCUCUGUUGCCAGGGAAGAACGCUUGCGAGAGUAACUGCUCGCAAGGUGACGGUACCUGAGAAGAAAGCCCCGGC
UAACUACGUGCCAGCAGCCGCGGUAAUACGUAGGGGGCAAGCGUUGUCCCGGAAUUAUUGGGCGUAAAGCGCGCGCAGGC
GGCCUUGUAAGUCUGUCGUUUAAACUCGGAGCUCAACUUCGAGUCGCGAUGGAAACUGCAAAGCUUGAGUGCAGAAGAGG
AAAGUGGAAUUCCACGUGUAGCGGUGAAAUGCGUAGAGAUGUGGAGGAACACCAGUGGCGAAGGCGACUUUCUGGGCUGU
AACUGACGCUGAGGCGCGAAAGCGUGGGGAGCAAACAGGAUUAGAUACCCUGGUAGUCCACGCCGUAAACGAUGAAUGCU
AGGUGUUAGGGGUUUCGAUACCCUUGGUGCCGAAGUUAACACAUUAAGCAUUCCGCCUGGGGAGUACGGUCGCAAGACUG
AAACUCAAAGGAAUUGACGGGGACCCGCACAAGCAGUGGAGUAUGUGGUUUAAUUCGAAGCAACGCGAAGAACCUUACCA
GGUCUUGACAUCCCUCUGACCGUCCUAGAGAUAGGGCUUUCCUUCGGGACAGAGGAGACAGGUGGUGCAUGGUUGUCGUC
AGCUCGUGUCGUGAGAUGUUGGGUUAAGUCCCGCAACGAGCGCAACCCUUGAUCUUAGUUGCCAGCACUUUGGGUGGGCA
CUCUAGGAUGACUGCCGGUGACAAACCGGAGGAAGGUGGGGAUGACGUCAAAUCAUCAUGCCCCUUAUGACCUGGGCUAC
ACACGUACUACAAUGGCCGAUACAACGGGAAGCGAAACCGCGAGGUGGAGCCAAUCCUAUCAAAGUCGGUCUCAGUUCGG
AUUGCAGGCUGCAACUCGCCUGCAUGAAGUCGGAAUUGCUAGUAAUCGCGGAUCAGCAUGCCGCGGUGAAUACGUUCCCG
GGUCUUGUACACACCGCCCGUCACACCACGAGAGUUUACAACACCCGAAGCCGGUGGGGUAACCGCAAGGAGCCAGCCGU
CGAAGGUGGGGUAGAUGAUUGG

[h.mon]

No it doesn't, because it lacks the semicolon ( ; ) at the end of the species / strain name in the fasta header. I've used the semicolon because all the fasta headers you provided had a semicolon, and it was a good way to, e.g., select "Paenibacillus borealis" (which is included in your list) and avoid "Paenibacillus borealisensis" (which I made up and is not included in your list). That is, it was just a way of avoiding longer names which could contain shorter ones, but were not in the list.

Anyways, if you change:

push(@headers, $_.";");

to:

push(@headers, $_);

it should work again.

Answer 3

[Juke34]

BASH

IFS=$'\n'
for i in $(cat headersFILE)
do
  line=$(grep -nr "$i" multiFastaFILE)
  if [[ ! -z $line ]]
  then
    for j in $line
    do
      lineNumber=$(echo $j | cut -d':' -f1)
      sed -n "$lineNumber p" multiFastaFILE
      awk -v nb=$lineNumber 'NR > nb {if ($0 ~ ">") exit; else print $0 }' multiFastaFILE
    done
  fi
done

Where headersFILE is the file containing one header per line ; multiFastaFILE is obviously your multifasta file.

Every time a fasta sequence header of your multiFastaFILE contains a name described in your headersFILE, the header and its sequence will be displayed.

[fhsantanna]

It worked very well! Thanks!

Answer 4

[shelly.deforte]

You could use Biopython's fasta parser. I didn't test this, but this would be the basic idea:

from Bio import SeqIO
handle = open("example.fasta", "rU")
organism_name = "Paenibacilus borealis"
for record in SeqIO.parse(handle, "fasta"):
    if organism_name in record.id:
        # Do whatever you want with these
        print record.id
        print record.seq
handle.close()

There's more info here: http://biopython.org/wiki/SeqIO

[adeel.maliks20]

I have a similar problem, I have many species name in one fasta file and I want them to copy those files which contain the desired species sequence.

The file is like this

>CanaanFir115301    PQPIVPPAASTSSEATTSTDVPQPIVHPAASTSSEATTSPVAFPARDNVAGTSTTSTDVPQPIVPPAASTSSEATTSPVA
FPARDNVAGTSSSSSVSKASPLIHDIVAAARAVVTDADLDDLIGNIFSLLEHIKAKKTETKSS
>FraserFir31480
MGKTTLADAVFARVYVEGCKYSMVRLFDDVTSTPSTSHIVKLQKYILEDLKMGTPEETTHE
>CanaanFir58880
SSTDVLPESNEIAPVDHFSVTSTAAPNATDHGPQPIVPPDASSSSTDVLPESNEIAPVDHFSVTSTAAPNATDDALKHPQ
HGPQPILPPDASSSSTDVLPDSNEPQPIVHPDLSSSSPLAFPAHDEVYQVSEAASSSNAATDNHKVADDQLQGTPSSSTV
SNAPHLINNILTAANALDT
>BalsamFir4479
MEKKSCRCRLFSCFQCCFSAQRPDPSVTSTAATDDALHHRPQPIVPSAASTSSEPTTSTDVPQPIVHPAASASSKFSCFQ
CCFSAQRPDPSVTSTAATDDALHHRPQPIVPSAASTSSEPTTSADVPQPIVHPAASTSSEATT

I want to copy the files which contain the sequence (Balsam, Fraser and Canaan) species. The rest will be discarded.

How can I do that?

Answer 5

[venu]

As your list is containing 'Paenibaci' in common, I have followed this approach

grep '^>' input.fasta > headers.txt
sed -n '/Paenibaci/p' headers.txt > filter_headers.txt
$ ./faSomeRecords input.faa filter_headers.txt output.faa

perl code to extract sequences from multi-line fasta works on all test files but not on research files

[fhsantanna]

Not worked. Output file is empty.

[venu]

Can you provide some more organism names from the organisms file? You have only provided one name, post 5-10 line from organism file

[fhsantanna]

Edited.

[venu]

try the edited one..it might work

[fhsantanna]

I think you have not understood. I already have a "header" file (organism name), which was not derived from a fasta file (so, grep and sed commands are unnecessary). I want to collect those sequences that contain specif organism names ("headers") from a FASTA file.

Answer 6

[Matt Shirley]

pip install --user pyfaidx 
xargs faidx 16s.fasta < organisms.txt

[shelly.deforte]

I wasn't familiar with that package and it looks really nice. Thanks!

[fhsantanna]

Firstly, the pip install did not work. I had to install using setup.py.

When I try to run the script, I get the following error:

Traceback (most recent call last):
  File "/usr/local/bin/faidx", line 9, in <module>
    load_entry_point('pyfaidx==0.4.0', 'console_scripts', 'faidx')()
  File "build/bdist.linux-x86_64/egg/pyfaidx/cli.py", line 149, in main
  File "build/bdist.linux-x86_64/egg/pyfaidx/cli.py", line 40, in write_sequence
  File "build/bdist.linux-x86_64/egg/pyfaidx/cli.py", line 50, in fetch_sequence
KeyError: 'Paenibacillus_alginolyticus_DSM_5050'

My command was:

xargs faidx rdp_download_2882seqs.fasta < organisms.txt

where organisms.txt is the list of organism names.

[Matt Shirley]

Regarding the pip install issue, you probably don't have pip installed, but it seems that you figured that out. The KeyError is telling you that the organism name Paenibacillus_alginolyticus_DSM_5050 does not exist as a FASTA header line in your file. I should probably catch this and warn the user instead of having the script bail.

[fhsantanna]

Yes, I had pip installed, and it installed pyfaidx. However faidx was not executable it in terminal, only in python.

In my list of species, each line has a organism name. Paenibacillus_alginolyticus_DSM_5050 was the first one. So, if the first name is not located in my fasta file, faidx does not try to locate other names?

[Matt Shirley]

That is correct, but I just submitted a change that fixes that behavior. The current github master branch will just warn the user and move to the next query without exiting with an error code.

[Matt Shirley]

I've just pushed a new version of pyfaidx (v0.4.1) that will produce a warning for each sequence that is not found in your file, but will continue without exiting. You can redirect stderr to a log file to figure out which organisms are missing:

xargs faidx rdp_download_2882seqs.fasta < organisms.txt 2> missing.log

[fhsantanna]

Not worked for me. Here is the log:

..
warning: Paenibacillus not found in file
warning: sp. not found in file
warning: 1-49 not found in file
..

Of course the strain Paenibacillus sp. 1-49 is in rdp_download_2882seqs.fasta and organisms.txt files.

[Matt Shirley]

Ah, yes. There can be no spaces in the sequence names. "xargs" passes each line in organisms.txt as the arguments for "faidx". If there are spaces, they are interpreted as separate sequence names. Generally it is advisable to not include spaces in the FASTA sequence names, as most programs will truncate anything after the first space.

[fhsantanna]

Still not working. Here is the log:

warning: Paenibacillus_sp._1-49 not found in file

In my rdp_download_2882seqs.fasta I have the following sequence:

>ASRY01000111.1.1112_Paenibacillus_sp._1-49
GAGUAACACGUAGGCAACCUGCCCAUGAGACUGGGAUAACUACCGGAAACGGUAGCUAAUACCGGAUACAUCCUUUCCCU
GCAUGGGGAGGGGAGGAAAGACGGAGCAAUCUGUCACUGAUGGAUGGGCCUGCGGCGCAUUAGCUAGUUGGUGGGGUAAU

[Matt Shirley]

The issue here is that you are expecting to find an inexact match. You can do this using:

xargs faidx rdp_download_2882seqs.fasta --regex  < organisms.txt 2> missing.log

The --regex option specifies that the sequence names are regular expressions patterns to be matched.

[fhsantanna]

Traceback (most recent call last):
  File "/usr/local/bin/faidx", line 9, in <module>
    load_entry_point('pyfaidx==0.4.1', 'console_scripts', 'faidx')()
  File "build/bdist.linux-x86_64/egg/pyfaidx/cli.py", line 153, in main
  File "build/bdist.linux-x86_64/egg/pyfaidx/cli.py", line 16, in write_sequence
  File "build/bdist.linux-x86_64/egg/pyfaidx/__init__.py", line 577, in __init__
  File "build/bdist.linux-x86_64/egg/pyfaidx/__init__.py", line 214, in __init__
IOError: [Errno 2] No such file or directory: 'Paenibacillus_alginolyticus_DSM_5050'

[Matt Shirley]

See the updated command above. Note that order of arguments matters.

[fhsantanna]

faidx: error: unrecognized arguments: Paenibacillus_sp._1-49

I also have tried with simpler data, but no output has been generated (I also have utilized the -o option).

[Matt Shirley]

For future online interactions, it is usually preferred to tell people what you did (what commands you tried), the error, and show some of the data you're working with. At this point I see no reason why you should not have a solution that works.

[fhsantanna]

Firstly, thanks. I had two solutions, as you can see. I just tried to give you a feedback, since you expend some time helping me.

I have tested your script with the data I provided above (the only difference is the presence of underscores instead of spaces). The command is exactly what you send me (the file names are the same):

xargs faidx rdp_download_2882seqs.fasta --regex  < organisms.txt 2> missing.log > filtered.out

Here is the error message:

faidx: error: unrecognized arguments: Paenibacillus_sp._1-18 Paenibacillus_sp._1-49 Paenibacillus_sp._A9 Paenibacillus_sp._Aloe-11

Certainly your script recognized "Paenibacilus_borealis), but the output was empty.

I tested the same data with fasgrep, and I had a positive result.