[Bioperl-l] counting gaps in a column of alignment

subha kalyanamoorthy sksweety24 at gmail.com
Thu Jul 25 00:37:02 EDT 2013


Hi there,

I am a new bioperl user.  I made a script to count the nucleotide
composition in each column of the alignment. I am able to get the count for
the nucleotides, but not for the gap characters from the following program.
I would greatly appreciate any suggestions to correct this program.

Thanks.

#***************************My Program**********************************

#!/bin/perl -w
use strict;
use warnings;

use List::Util 'max';
use Bio::SimpleAlign;
use Bio::Align::AlignI;
use Bio::AlignIO;
use Bio::SeqIO;

my $in= Bio::AlignIO->new( -file => "seq.fst", -format => "fasta");

my $align = $in->next_aln();

print "column\tA's\tT's\C's\G's\n";

for (my $i = 1; $i <= $align->length; $i++) {

 my %count;

 my $seqs = $align->slice($i,$i);

my $gap_char = $seqs->gap_char();

    my $count_A=0;
    my $count_C=0;
    my $count_T=0;
    my $count_G=0;
    my $count_N=0;
    my $count_gap=0;

    foreach my $seq ($seqs->each_seq) {

    my $col=$seq->seq;

        if ($col eq 'A'){
        $count_A++;
        }elsif ($col eq 'C'){
        $count_C++;
        }elsif ($col eq 'T'){
        $count_T++;
        }elsif ($col eq 'G'){
        $count_G++;
        }elsif ( $col eq 'N'){
        $count_N++;
        }elsif ($col =~ m/^\Q$gap_char\E$/){
        $count_gap++;
}
      $count{$seq->seq} += 1;
 }

print"$i\t$count_A\t$count_T\t$count_C\t$count_G\n";

}

#***********************************************************************`


More information about the Bioperl-l mailing list