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 has 'format'     => (is => 'rw',
  42                      isa => 'Str',
  43                      traits => [ 'Getopt' ],
  44                      default => 'txt',
  45     );
  46
  47
  48 sub run {
  49     my $self = shift;
  50
  51     my $io = Bio::SeqIO->new(-format=>'largefasta', -file=>$self->fasta_file());
  52
  53     my $gap_no = 1;
  54     my $old_id = "";
  55
  56     while (my $s = $io->next_seq()) {
  57         my $seq = $s->seq();
  58         my $id = $s->id();
  59
  60         if ($id ne $old_id) {
  61             $gap_no =1;
  62         }
  63
  64         warn "Processing sequence $id (".$s->length()." nucleotides)...\n";
  65
  66         my $n_region_start = 0;
  67         my $n_region_end = 0;
  68         foreach my $i (1..$s->length()) {
  69             my $nuc = $s->subseq($i, $i);
  70
  71             if ($nuc=~/n/i) {
  72                 if (!$n_region_start) { $n_region_start=$i; }
  73             }
  74
  75             else {
  76                 if ($n_region_start) { $n_region_end = $i; }
  77
  78                 my $gap_size = $n_region_end - $n_region_start;
  79                 if ($gap_size >= $self->min_gap_size()) {
  80
  81                     if ($self->format=~ /gff/i) {
  82                         print join "\t", ($id, "getgaps", "gap", $n_region_start, $n_region_end, 1, ".", ".", "ID=$id\_$gap_no;Parent=$id"); print "\n";
  83                     }
  84                     else {
  85
  86                         print "$id\_"; printf "%06d", "$gap_no"; print "\t$id\t$n_region_start\t$n_region_end\t$gap_size\n";
  87                     }
  88                         $gap_no++;
  89
  90                 }
  91
  92                 $n_region_start = 0;
  93                 $n_region_end = 0;
  94
  95             }
  96
  97
  98         }
  99             $old_id = $id;
 100     }
 101
 102 }
 103
 104
 105
 106 1;