SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Extract sequence from multi fasta file with PERL andreitudor Bioinformatics 27 07-07-2019 07:45 AM
FASTA sequence From large BAM file mez Bioinformatics 9 01-13-2013 05:42 AM
Perl: get specific base from FASTA file. njh_TO Bioinformatics 6 02-02-2012 05:34 AM
subsetting a bam file with specific alingment length joseph Bioinformatics 4 12-28-2011 07:39 PM
Find all occurrences of a sequence in a fasta file dphansti Bioinformatics 3 12-06-2011 06:11 AM

Reply
 
Thread Tools
Old 07-11-2012, 08:21 AM   #1
entomology
Member
 
Location: zhejiang, china

Join Date: May 2011
Posts: 34
Default how to get specific length sequence from a fasta file

I'm not so familiar with perl, is there script can fetch specific length of sequence from a fasta file? many thanks!
entomology is offline   Reply With Quote
Old 07-11-2012, 09:37 AM   #2
maubp
Peter (Biopython etc)
 
Location: Dundee, Scotland, UK

Join Date: Jul 2009
Posts: 1,543
Default

You mean loop over the records, outputting just those of a certain length? That should just be a couple of lines using BioPerl, Biopython, etc.

If you're not familiar with Perl, what can you program in?
maubp is offline   Reply With Quote
Old 07-11-2012, 04:12 PM   #3
entomology
Member
 
Location: zhejiang, china

Join Date: May 2011
Posts: 34
Default

thanks, I'll try to learn basic perl first, parse with fasta is first step, maybe bioperl is more easy to handle.
entomology is offline   Reply With Quote
Old 07-12-2012, 12:37 PM   #4
rexxi
Member
 
Location: Santiago, Chile

Join Date: Jun 2012
Posts: 20
Default

its name is bp_aacomp.pl, this script can count the aminoacids of a FASTA format file. This doesn't work good with multifasta file, because it takes count of all the aminoacids of the file. If you want to count the aminoacids of 1 sequence I would recomend you. the usage is really simple

perl bp_aacomp.pl yourfile.fasta

if you need to bp, you can try this LINK where the is another script.

Quote:
#!perl
use strict;
use warnings;
use Carp;

use Bio::SeqIO;
use Getopt::Long;
use Bio::SeqUtils;
use Bio::Tools::IUPAC;

my $table = new Bio::SeqUtils;
my @BASES = $table->valid_aa(0);
my %all = $table->valid_aa(2);
my ($file,$format,$help) = ( undef, 'fasta');
GetOptions(
'i|in:s' => \$file,
'f|format:s' => \$format,
'h|help|?' => \$help,
);

my $USAGE = "usage: aacomp [-f format] filename\n\tdefault format is fasta\n";
$file = shift unless $file;

die $USAGE if $help;

my $seqin;
if( defined $file ) {
print "Could not open file [$file]\n$USAGE" and exit unless -e $file;
$seqin = new Bio::SeqIO(-format => $format,
-file => $file);
} else {
$seqin = new Bio::SeqIO(-format => $format,
-fh => \*STDIN);
}

my %composition;
my $total;
foreach my $base ( @BASES ) {
$composition{$base} = 0;
}
while ( my $seq = $seqin->next_seq ) {
if( $seq->alphabet ne 'protein' ) {
confess("Must only provide amino acid sequences to aacomp...skipping this seq");
next;
}
foreach my $base ( split(//,$seq->seq()) ) {
$composition{uc $base}++;
$total++;
}
}

printf("%d aa\n",$total);
printf("%5s %4s\n", 'aa', '#' );
my $ct = 0;
foreach my $base ( @BASES ) {
printf(" %s %s %3d\n", $base, $all{$base}, $composition{$base} );
$ct += $composition{$base};
}
printf( "%6s %s\n", '','-'x5);
printf( "%6s %3d\n", '',$ct);


__END__


=head1 NAME

bp_aacomp - amino acid composition of protein sequences

=head1 SYNOPSIS

bp_aacomp [-f/--format FORMAT] [-h/--help] filename
or
bp_aacomp [-f/--format FORMAT] < filename
or
bp_aacomp [-f/--format FORMAT] -i filename

=head1 DESCRIPTION

This scripts prints out the count of amino acids over all protein
sequences from the input file.

=head1 OPTIONS

The default sequence format is fasta.

The sequence input can be provided using any of the three methods:

=over 3

=item unnamed argument

bp_aacomp filename

=item named argument

bp_aacomp -i filename

=item standard input

bp_aacomp < filename

=back

=head1 FEEDBACK

=head2 Mailing Lists

User feedback is an integral part of the evolution of this and other
Bioperl modules. Send your comments and suggestions preferably to
the Bioperl mailing list. Your participation is much appreciated.

bioperl-l@bioperl.org - General discussion
http://bioperl.org/wiki/Mailing_lists - About the mailing lists

=head2 Reporting Bugs

Report bugs to the Bioperl bug tracking system to help us keep track
of the bugs and their resolution. Bug reports can be submitted via the
web:

https://redmine.open-bio.org/projects/bioperl/

=head1 AUTHOR - Jason Stajich

Email jason@bioperl.org

=head1 HISTORY

Based on aacomp.c from an old version of EMBOSS

=cut
rexxi is offline   Reply With Quote
Old 07-12-2012, 03:37 PM   #5
adaptivegenome
Super Moderator
 
Location: US

Join Date: Nov 2009
Posts: 437
Default

Samtools will do this if you just want a specific sequence from a Fasta file.
adaptivegenome is offline   Reply With Quote
Old 07-12-2012, 03:59 PM   #6
entomology
Member
 
Location: zhejiang, china

Join Date: May 2011
Posts: 34
Default

Thank you all the same!
entomology is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 04:17 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2020, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO