Bio::Tools::SeqStats to count bases

classic Classic list List threaded Threaded
8 messages Options
Reply | Threaded
Open this post in threaded view
|

Bio::Tools::SeqStats to count bases

Slym-2
The result I get is:

Number of bases of type A =
Number of bases of type C =
Number of bases of type G =
Number of bases of type T =

i.e. There's no expected values.
Please help!

#! /usr/bin/perl

use Bio::Tools::SeqStats;
use Bio::Seq;

open (FILE, "seq.fasta");
@array = <FILE>;

# Removing first line of fasta

shift (@array);
$array = join('',@array);
open (FILE2, ">>seq2.fasta");
print FILE2 "$array";

$seqobj = Bio::PrimarySeq->new( -file => "sekw2.fasta",
- alphabet => 'dna',);


my $seq_stats = Bio::Tools::SeqStats->new(-seq=>$seqobj);

my $monomer_ref = $seq_stats->count_monomers();

foreach $base (sort keys %$monomer_ref) {
print "Liczba zasad typu ", $base," = ", $monomer_ref{$base},"\n";
}


_______________________________________________
Bioperl-l mailing list
[hidden email]
http://lists.open-bio.org/mailman/listinfo/bioperl-l
Reply | Threaded
Open this post in threaded view
|

Re: Bio::Tools::SeqStats to count bases

Slym-2
Thanks Roy,

It still doesn't seem to produce anything. :/
_______________________________________________
Bioperl-l mailing list
[hidden email]
http://lists.open-bio.org/mailman/listinfo/bioperl-l
Reply | Threaded
Open this post in threaded view
|

Re: Bio::Tools::SeqStats to count bases

Roy Chaudhuri-3
Sorry, I'd missed another problem in your code - you are trying to
load a fasta file using Bio::PrimarySeq. To read sequence data from a
file you should use Bio::SeqIO, see:

http://www.bioperl.org/wiki/HOWTO:Beginners#Retrieving_a_sequence_from_a_file
http://www.bioperl.org/wiki/HOWTO:SeqIO

Cheers,
Roy.
_______________________________________________
Bioperl-l mailing list
[hidden email]
http://lists.open-bio.org/mailman/listinfo/bioperl-l
Reply | Threaded
Open this post in threaded view
|

Re: Bio::Tools::SeqStats to count bases

Adam Sjøgren
In reply to this post by Slym-2
On Mon, 4 Feb 2013 07:39:19 -0800 (PST), Slym wrote:

> #! /usr/bin/perl

> use Bio::Tools::SeqStats;
> use Bio::Seq;

It can be a good idea to add "use strict; use warnings;" to the top of
your script. At least two problems in your program would have been
caught by perl if you had.

> open (FILE, "seq.fasta");

Using (global) literal filehandles and the two parameter open() is
somewhat outdated, a more current way to do it could be:

  open my $fh, '<', 'seq.fasta';

> @array = <FILE>;

> # Removing first line of fasta

> shift (@array);
> $array = join('',@array);
> open (FILE2, ">>seq2.fasta");
> print FILE2 "$array";

Note that you are writing just the sequence to your seq2.fasta file
here, so the new file isn't really a fasta file.

> $seqobj = Bio::PrimarySeq->new( -file => "sekw2.fasta",
> - alphabet => 'dna',);

Bio::PrimarySeq doesn't take a '-file' parameter. Also, note that the
filename is different than before "sekw2" vs. "seq2"!

Either you should use Bio::SeqIO with a '-file' parameter, or you can
use Bio::PrimarySeq with a '-seq' parameter.

> my $seq_stats = Bio::Tools::SeqStats->new(-seq=>$seqobj);

> my $monomer_ref = $seq_stats->count_monomers();

> foreach $base (sort keys %$monomer_ref) {
> print "Liczba zasad typu ", $base," = ", $monomer_ref{$base},"\n";

Here you wanted $monomer_ref->{$base}, as %monomer_ref isn't mentioned
anywhere else.

> }

Here is a complete version of your script - I chose to use Bio::SeqIO -
that works:

  #!/usr/bin/perl

  use strict;
  use warnings;

  use Bio::SeqIO;
  use Bio::Tools::SeqStats;

  my $io=Bio::SeqIO->new(-file=>'seq.fasta', -alphabet=>'dna');
  my $seqobj=$io->next_seq; # Get the first sequence from the file

  my $seq_stats = Bio::Tools::SeqStats->new(-seq=>$seqobj);
  my $monomer_ref = $seq_stats->count_monomers();
  foreach my $base (sort keys %$monomer_ref) {
      print "Liczba zasad typu ", $base," = ", $monomer_ref->{$base},"\n";
  }

E.g.:

  $ cat seq.fasta
  >test
  aaaacccggt
  $ ./slym.pl
  Liczba zasad typu A = 4
  Liczba zasad typu C = 3
  Liczba zasad typu G = 2
  Liczba zasad typu T = 1
  $


  Best regards,

    Adam

--
 "Grittings. Ma nam is Kahlfin."                              Adam Sjøgren
                                                         [hidden email]

_______________________________________________
Bioperl-l mailing list
[hidden email]
http://lists.open-bio.org/mailman/listinfo/bioperl-l
Reply | Threaded
Open this post in threaded view
|

Re: Bio::Tools::SeqStats to count bases

Slym-2
In reply to this post by Roy Chaudhuri-3
The thing is, if I use Bio::SeqIO then  Bio::Tools::SeqStats produces an
error (saying that it wants input provided by Bio::PrimarySeq).
(btw in this line
 $seqobj = Bio::PrimarySeq->new( -file => "sekw2.fasta", - alphabet =>
'dna',);
there's a typo "sekw2" instead of "seq2" but this is correct in my original
code).

_______________________________________________
Bioperl-l mailing list
[hidden email]
http://lists.open-bio.org/mailman/listinfo/bioperl-l
Reply | Threaded
Open this post in threaded view
|

Re: Bio::Tools::SeqStats to count bases

Fields, Christopher J
Please make sure and read both Roy's and Adam's responses all the way through; Bio::SeqIO is not a sequence object but the front-end for format parsing (e.g. FASTA, etc).  Bio::PrimarySeq does not have a '-file' parameter, Bio::SeqIO does.  

If SeqStats truly doesn't work with Bio::Seq we can fix that, but according to Adam he has tested using Bio::SeqIO out and it seems to work.

chris

On Feb 4, 2013, at 12:02 PM, Slym <[hidden email]>
 wrote:

> The thing is, if I use Bio::SeqIO then  Bio::Tools::SeqStats produces an
> error (saying that it wants input provided by Bio::PrimarySeq).
> (btw in this line
> $seqobj = Bio::PrimarySeq->new( -file => "sekw2.fasta", - alphabet =>
> 'dna',);
> there's a typo "sekw2" instead of "seq2" but this is correct in my original
> code).
>
> _______________________________________________
> Bioperl-l mailing list
> [hidden email]
> http://lists.open-bio.org/mailman/listinfo/bioperl-l


_______________________________________________
Bioperl-l mailing list
[hidden email]
http://lists.open-bio.org/mailman/listinfo/bioperl-l
Reply | Threaded
Open this post in threaded view
|

Re: Bio::Tools::SeqStats to count bases

Adam Sjøgren
In reply to this post by Slym-2
On Mon, 4 Feb 2013 10:02:29 -0800 (PST), Slym wrote:

> The thing is, if I use Bio::SeqIO then  Bio::Tools::SeqStats produces an
> error (saying that it wants input provided by Bio::PrimarySeq).

That sounds like you forgot to call ->next_seq() on the Bio::SeqIO
object - to get a sequence object - please see the complete, working
example I sent earlier.


  Best regards,

    Adam

--
 "Denial springs eternal."                                    Adam Sjøgren
                                                         [hidden email]

_______________________________________________
Bioperl-l mailing list
[hidden email]
http://lists.open-bio.org/mailman/listinfo/bioperl-l
Reply | Threaded
Open this post in threaded view
|

Re: Bio::Tools::SeqStats to count bases

Slym-2
Everything's working now! Thank you very much, especially to you Adam!


>
_______________________________________________
Bioperl-l mailing list
[hidden email]
http://lists.open-bio.org/mailman/listinfo/bioperl-l