lib/CXGN/Tools/GetGaps.pm

   1
   2 =head1 NAME
   3
   4 CXGN::Tools::GetGaps - a script to obtain information about gaps (indicated as N's) in a fasta file (such as the tomato genome)
   5
   6 =head1 DESCRIPTION
   7
   8 Implemented as a MooseX::Runnable script to be run as follows:
   9
  10 mx-run CXGN::Tools::GetGaps --fasta_file sequences.fasta --min_gap_size 20 > output.txt
  11
  12 =head1 AUTHOR
  13
  14 Lukas Mueller
  15
  16
  17 =cut
  18
  19
  20
  21 package CXGN::Tools::GetGaps;
  22 use Moose;
  23
  24 with 'MooseX::Runnable';
  25 with 'MooseX::Getopt';
  26
  27 use Bio::SeqIO;
  28
  29 has 'min_gap_size' => (is => 'rw',
  30                        isa => 'Int',
  31                        traits=> ['Getopt'],
  32                        default=>10,
  33     );
  34
  35
  36 has 'fasta_file' => (is => 'rw',
  37                      isa => 'Str',
  38                      traits=>['Getopt'],
  39     );
  40
  41
  42
  43 sub run {
  44     my $self = shift;
  45
  46     my $io = Bio::SeqIO->new(-format=>'largefasta', -file=>$self->fasta_file());
  47
  48     my $gap_no = 1;
  49     my $old_id = "";
  50
  51     while (my $s = $io->next_seq()) {
  52         my $seq = $s->seq();
  53         my $id = $s->id();
  54
  55         if ($id ne $old_id) {
  56             $gap_no =1;
  57         }
  58
  59         warn "Processing sequence $id (".$s->length()." nucleotides)...\n";
  60
  61         my $n_region_start = 0;
  62         my $n_region_end = 0;
  63         foreach my $i (1..$s->length()) {
  64             my $nuc = $s->subseq($i, $i);
  65
  66             if ($nuc=~/n/i) {
  67                 if (!$n_region_start) { $n_region_start=$i; }
  68             }
  69
  70             else {
  71                 if ($n_region_start) { $n_region_end = $i; }
  72
  73                 my $gap_size = $n_region_end - $n_region_start + 1;
  74                 if ($gap_size >= $self->min_gap_size()) {
  75
  76                     print "$id\_"; printf "%06d", "$gap_no"; print "\t$id\t$n_region_start\t$n_region_end\t$gap_size\n";
  77                     $gap_no++;
  78
  79                 }
  80
  81                 $n_region_start = 0;
  82                 $n_region_end = 0;
  83
  84             }
  85
  86
  87         }
  88     }
  89 }
  90
  91
  92
  93 1;