Extract sequence from multi fasta file with PERL

[andreitudor]

Hello,

I want to make a script that will extract a sequence from a multifasta file. I have made a file with 19 genomes from NCBI and I want to be able to extract a genome that I need from it. I have played a little with Bio::SeqIO and learned how to create a object that has the hole file in it and then right the whole information in another file.
What i do not know how to do is how to select only the fasta sequence from the genome that I want and write it onto a file.
I think I have to create a third object that would contain only that genome, but I still don't know how to read that part from my sequence object...

Thanks in advance.

[boetsie]

Hi,

I've attached a script which can do this. If i understand it correctly you have a file like;

>chr1
AGCTGATGATAGT...
>chr2
ACAAAATAGTCGAT....
>chr3
....

And your perl script would be something like;

perl extractSequence.pl genomefile.fa chr1

where 'chr1' corresponds to a sequence named chr1 (indicated by chr1)?

Say you have a more complicated file like;

>chr1_coverage1000_length100
AGATGTATGTTAGA

You can do something like;

perl extractSequence.pl genomefile.fa chr1_.

which will extract all the sequences containing the header chr1_

To store the results, do;

perl extractSequence.pl genomefile.fa chr1 > filename.txt

If this is what you want, you can use my script.

Boetsie

[maasha]

This can be done with Biopieces (www.biopieces.org) like this:

read_fasta -i ncbi_genomes.fna | grab -p my_genome | write_fasta -o my_genome.fna -x

Cheers


Martin

[NM_010117]

This script might be a little easier to wrap your head around (there's less code anyway) and does the same thing as boetsie's:

#!/usr/bin/perl -w
use strict;
my $file = shift;
my $pattern = shift;
my $data = `cat $file`;
my ($fa) = $data =~ /(>$pattern[^>]+)/s;
print $fa;

Mind you boetsie's implementation has the benefit of not having to read in the entire file every time (i.e. her implementation "slurps" in the file to avoid needing to use a large amount of memory). However, this of course isn't an issue with small files (e.g. bacterial genomes).

[ssully]

another way (if you're already set up for local BLAST searches) -- use 'fastacmd' to extract whatever records you need from a fasta db you've created using 'formatdb'

(this refers to the old NCBI BLAST -- the newer NCBI 'blast+' uses different but analogous toolset:

http://www.ncbi.nlm.nih.gov/staff/ta.../pc_setup.html)

[apc2010]

If you need sequences extracted from a multi-FASTA and are open to using a pre-existing tool, I would also suggest either the faSomeRecords or faOneRecord command line utilities from UCSC.

They have versions of this tool for OSX and Linux. Here is a link to the executable downloads:

http://hgdownload.cse.ucsc.edu/admin/exe/

The difference between the two: faOneRecord takes the sequence name to extract from the command line, faSomeRecords reads in a file of 1 or more sequence names to extract from the multi-FASTA.

Usage:

================================================================
========   faOneRecord   ====================================
================================================================
faOneRecord - Extract a single record from a .FA file
usage:
   faOneRecord in.fa recordName

================================================================
========   faSomeRecords   ====================================
================================================================
faSomeRecords - Extract multiple fa records
usage:
   faSomeRecords in.fa listFile out.fa
options:
   -exclude - output sequences not in the list file.

[kmcarr]

There was a previous thread which explained how to do this with the BLAST tools (old and new)

http://seqanswers.com/forums/showthread.php?t=8805&highlight=fastacmd

The problem with using Bio::SeqIO (or a simple perl script in general) is that there is no indexing of the file. This means that to pull out a specific sequence you have to scan through the file from the beginning until you find the sequence you are interested in. If your input file isn't too huge and you only need to do it for a small number of sequences or very infrequently that would work. But if the input file is ginormous (e.g. if your 19 genomes from NCBI are all big, mammalian genomes) and you want to repeatedly extract individual sequences then an indexed method like the BLAST tools is advisable.

[krobison]

There's also Bio:B::Fasta module for Perl. One additional advantage to this tool is that you can easily get a specific range of a sequence as a subsequence.

I'm guilty of frequently doing this as a one-liner -- slow on a large database and a pain for a long list of sequences (or if you get the regex wrong), but sometimes it's easier than remembering a more specific tool.

perl -n -e '$on=(/^>(idA|idB|idC)/) if (/^>/); print $_ if ($on);' sequences.fa

[andreitudor]

Hey guys,

Thanks for your answers. I do not want to use the tools provided from blast+ because I want to do it on my own. I looked into BIO:: DB:Fasta and I decided to use it, since it has a function that indexes your fasta file and then you can search by ID.
Now I have one problem, when I create the indexed file BIO:: DB:Fasta gives me this error :

indexing was interrupted, so unlinking reference.fna.index at /usr/share/perl5/Bio/DB/Fasta.pm line 1053.
Was not able to open files!!
Full error is


------------- EXCEPTION: Bio::Root::Exception -------------
MSG: Each line of the fasta entry must be the same length except the last.
Line above #46357 '
..' is 69 != 71 chars.
STACK: Error::throw
STACK: Bio::Root::Root::throw /usr/share/perl5/Bio/Root/Root.pm:368
STACK: Bio:B::Fasta::calculate_offsets /usr/share/perl5/Bio/DB/Fasta.pm:770
STACK: Bio:B::Fasta::index_file /usr/share/perl5/Bio/DB/Fasta.pm:680
STACK: Bio:B::Fasta::new /usr/share/perl5/Bio/DB/Fasta.pm:491
STACK: extract.pl:22
-----------------------------------------------------------

From this I understand that the lines in my fasta file are not equal. The only problem is that the line that is not the same length is the ID line. What can I do about this??

Andrei

[krobison]

Check the line lengths; the program is complaining that there is a sequence line of different length; it could well be some stray spaces; take a look at line #46357 with a text editor

This little program will find the numbers for all your "illegal" lines
[CODE]
#!/usr/bin/perl
use strict;
## find & print line number for all "illegal" FASTA lines
my %lengths=();
my $lastId=undef;
my $lineCount=0;
my $len=undef;
while ($_ = <>)
{
$lineCount++;
chomp;
if (/^>(\S+)/)
{
$lastId=$1;
if (defined $len) # last line of entry allowed to be ragged
{
delete $lengths{$len}->{$lineCount-1};
}
}
else
{
$len=length($_);
$lengths{$len}={} unless (ref $lengths{$len});
$lengths{$len}->{$lineCount}=$lastId;
}
}

my @lengths=sort {scalar(keys %{$lengths{$b}})<=>scalar(keys %{$lengths{$a}})} keys %lengths;
my $modalLength=shift(@lengths);
print "# modalLength=$modalLength ",scalar(keys %{$lengths{$modalLength}}),"\n";
foreach my $len(@lengths)
{
foreach my $lineCount(sort {$a<=>$b} keys %{$lengths{$len}})
{
print join("\t",$len,$lineCount,$lengths{$len}->{$lineCount}),"\n";
}
}
[/CODE

[mghita]

Hi

Hi,

I might sound dumb to people around here , but how do I install/open/use the command? I tried faSomeRecords in.fa listFile out.fa with my files and it says "-bash: faSomeRecords: command not found", I added it to my path too. I'm new in the area and have no clue

Madalina

If you need sequences extracted from a multi-FASTA and are open to using a pre-existing tool, I would also suggest either the faSomeRecords or faOneRecord command line utilities from UCSC.

They have versions of this tool for OSX and Linux. Here is a link to the executable downloads:

http://hgdownload.cse.ucsc.edu/admin/exe/

The difference between the two: faOneRecord takes the sequence name to extract from the command line, faSomeRecords reads in a file of 1 or more sequence names to extract from the multi-FASTA.

Usage:

================================================================
========   faOneRecord   ====================================
================================================================
faOneRecord - Extract a single record from a .FA file
usage:
   faOneRecord in.fa recordName

================================================================
========   faSomeRecords   ====================================
================================================================
faSomeRecords - Extract multiple fa records
usage:
   faSomeRecords in.fa listFile out.fa
options:
   -exclude - output sequences not in the list file.

[GenoMax]

Hi,

I might sound dumb to people around here , but how do I install/open/use the command? I tried faSomeRecords in.fa listFile out.fa with my files and it says "-bash: faSomeRecords: command not found", I added it to my path too. I'm new in the area and have no clue

Madalina

Madalina,

I will assume that you have downloaded the correct version of the compiled program for your operating system and saved it in a folder (or home dir).

You may need to make this file executable by adding the correct unix permission: chmod u+x /path_to/faSomeRecords (if you are in the same directory where this program is then use chmod u+x ./faSomeRecords).

Now you should be able to execute the file by following the example form earlier post.
e.g. /path_to/faSomeRecords your_data_file

[mghita]

I have added the program to my path and I set the permission right, but now I have another issue:
"You need the Rosetta software to run faSomeRecords. The Rosetta installer is in Optional Installs on your Mac OS X installation disc."

and I don't have Rosetta installed, or the CD for installation, so I don't know how to handle this problem. Any suggestions?


Thanks,
Madalina

Madalina,

I will assume that you have downloaded the correct version of the compiled program for your operating system and saved it in a folder (or home dir).

You may need to make this file executable by adding the correct unix permission: chmod u+x /path_to/faSomeRecords (if you are in the same directory where this program is then use chmod u+x ./faSomeRecords).

Now you should be able to execute the file by following the example form earlier post.
e.g. /path_to/faSomeRecords your_data_file

[GenoMax]

I have added the program to my path and I set the permission right, but now I have another issue:
"You need the Rosetta software to run faSomeRecords. The Rosetta installer is in Optional Installs on your Mac OS X installation disc."

and I don't have Rosetta installed, or the CD for installation, so I don't know how to handle this problem. Any suggestions?


Thanks,
Madalina

Madalina,

If you are connected to the internet you should automatically be offered the option to download rosetta and install it.

Do you have a PowerPC- or an intel-based Mac? What OS are you running?

[mghita]

Madalina,

If you are connected to the internet you should automatically be offered the option to download rosetta and install it.

Do you have a PowerPC- or an intel-based Mac? What OS are you running?

I have Mac OS X 10.6.8, 3.06 GHz. I just get that message in bash, I don't get any install option. I tried to download it, but it doesn't work.