code_perl/Maasha/Seq.pm

   1 package Maasha::Seq;
   2
   3 # Copyright (C) 2007 Martin A. Hansen.
   4
   5 # This program is free software; you can redistribute it and/or
   6 # modify it under the terms of the GNU General Public License
   7 # as published by the Free Software Foundation; either version 2
   8 # of the License, or (at your option) any later version.
   9
  10 # This program is distributed in the hope that it will be useful,
  11 # but WITHOUT ANY WARRANTY; without even the implied warranty of
  12 # MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  13 # GNU General Public License for more details.
  14
  15 # You should have received a copy of the GNU General Public License
  16 # along with this program; if not, write to the Free Software
  17 # Foundation, Inc., 51 Franklin Street, Fifth Floor, Boston, MA 02110-1301, USA.
  18
  19 # http://www.gnu.org/copyleft/gpl.html
  20
  21
  22 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> DESCRIPTION <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  23
  24
  25 # yak yak yak
  26
  27
  28 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  29
  30
  31 use warnings;
  32 use strict;
  33 use POSIX qw( islower );
  34 use Maasha::Common;
  35 use Data::Dumper;
  36 use IPC::Open2;
  37 use List::Util qw( shuffle );
  38 use Time::HiRes qw( gettimeofday );
  39
  40 use vars qw ( @ISA @EXPORT );
  41
  42 @ISA = qw( Exporter );
  43
  44
  45 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
  46
  47
  48 sub seq_guess_type
  49 {
  50     # Martin A. Hansen, May 2007.
  51
  52     # Makes a qualified guess on the type of a given squence.
  53
  54     my ( $seq,   # sequence to check
  55        ) = @_;
  56
  57     # returns string.
  58
  59     my ( $check_seq, $count );
  60
  61     if ( length $seq > 100 ) {
  62         $check_seq = substr $seq, 0, 100;
  63     } else {
  64         $check_seq = $seq;
  65     }
  66
  67     if ( $count = $check_seq =~ tr/FLPQIEflpqie// and $count > 0 ) {
  68         return "protein";
  69     } elsif ( $count = $check_seq =~ tr/Uu// and $count > 0 ) {
  70         return "rna";
  71     } else {
  72         return "dna";
  73     }
  74 }
  75
  76
  77 sub wrap
  78 {
  79     # Martin A. Hansen, July 2007.
  80
  81     # Wraps a given string reference accoring to given width.
  82
  83     my ( $strref,   # ref to sting to wrap
  84          $wrap,     # wrap width
  85        ) = @_;
  86
  87     # Returns nothing.
  88
  89     ${ $strref } =~ s/(.{$wrap})/$1\n/g;
  90
  91     chomp ${ $strref };
  92 }
  93
  94
  95 sub dna_revcomp
  96 {
  97     # Niels Larsen
  98     # modified Martin A. Hansen, March 2005.
  99
 100     # Returns the reverse complement of a dna sequence with preservation of case
 101     # according to this mapping,
 102     #
 103     # AGCUTRYWSMKHDVBNagcutrywsmkhdvbn
 104     # TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn
 105
 106     my ( $seq,   # seq
 107        ) = @_;
 108
 109     # returns string
 110
 111     $seq = reverse $seq;
 112
 113     $seq =~ tr/AGCUTRYWSMKHDVBNagcutrywsmkhdvbn/TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn/;
 114
 115     return $seq;
 116 }
 117
 118
 119 sub rna_revcomp
 120 {
 121     # Niels Larsen
 122     # modified Martin A. Hansen, March 2005.
 123
 124     # Returns the complement of a rna sequence with preservation of case
 125     # according to this mapping,
 126     #
 127     # AGCUTRYWSMKHDVBNagcutrywsmkhdvbn
 128     # UCGAAYRWSKMDHBVNucgaayrwskmdhbvn
 129
 130     my ( $seq,   # seq
 131        ) = @_;
 132
 133     $seq = reverse $seq;
 134
 135     $seq =~ tr/AGCUTRYWSMKHDVBNagcutrywsmkhdvbn/UCGAAYRWSKMDHBVNucgaayrwskmdhbvn/;
 136
 137     return $seq;
 138 }
 139
 140
 141 sub dna_comp
 142 {
 143     # Niels Larsen
 144     # modified Martin A. Hansen, March 2005.
 145
 146     # Returns the reverse complement of a dna sequence with preservation of case
 147     # according to this mapping,
 148     #
 149     # AGCUTRYWSMKHDVBNagcutrywsmkhdvbn
 150     # TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn
 151
 152     my ( $seqref,   # seqref
 153        ) = @_;
 154
 155     # Returns nothing.
 156
 157     ${ $seqref } =~ tr/AGCUTRYWSMKHDVBNagcutrywsmkhdvbn/TCGAAYRWSKMDHBVNtcgaayrwskmdhbvn/;
 158 }
 159
 160
 161 sub rna_comp
 162 {
 163     # Niels Larsen
 164     # modified Martin A. Hansen, March 2005.
 165
 166     # Returns the complement of a rna sequence with preservation of case
 167     # according to this mapping,
 168     #
 169     # AGCUTRYWSMKHDVBNagcutrywsmkhdvbn
 170     # UCGAAYRWSKMDHBVNucgaayrwskmdhbvn
 171
 172     my ( $seqref,   # seqref
 173        ) = @_;
 174
 175     # Returns nothing.
 176
 177     ${ $seqref } =~ tr/AGCUTRYWSMKHDVBNagcutrywsmkhdvbn/UCGAAYRWSKMDHBVNucgaayrwskmdhbvn/;
 178 }
 179
 180
 181 sub dna2rna
 182 {
 183     # Martin A. Hansen, March 2007
 184
 185     # Converts DNA sequence to RNA
 186
 187     my ( $seq,   # nucleotide sequence
 188        ) = @_;
 189
 190     # returns string
 191
 192     $seq =~ tr/Tt/Uu/;
 193
 194     return $seq;
 195 }
 196
 197
 198 sub rna2dna
 199 {
 200     # Martin A. Hansen, March 2007
 201
 202     # Converts RNA sequence to DNA
 203
 204     my ( $seq,   # nucleotide sequence
 205        ) = @_;
 206
 207     # returns string
 208
 209     $seq =~ tr/Uu/Tt/;
 210
 211     return $seq;
 212 }
 213
 214
 215 sub nuc2ambiguity
 216 {
 217     # Martin A. Hansen, March 2005.
 218
 219     # given a string of nucleotides
 220     # returns the corresponding ambiguity code
 221
 222     my ( $str,
 223          $type,    # DNA or RNA - DEFAULT DNA
 224        ) = @_;
 225
 226     my ( %hash, @nts, $key, $code, %nt_hash );
 227
 228     $str = uc $str;
 229
 230     if ( not $type or $type =~ /dna/i )
 231     {
 232         $str =~ s/N/ACGT/g;
 233         $str =~ s/B/CGT/g;
 234         $str =~ s/D/AGT/g;
 235         $str =~ s/H/ACT/g;
 236         $str =~ s/V/ACG/g;
 237         $str =~ s/K/GT/g;
 238         $str =~ s/Y/CT/g;
 239         $str =~ s/S/CG/g;
 240         $str =~ s/W/AT/g;
 241         $str =~ s/R/AG/g;
 242         $str =~ s/M/AC/g;
 243     }
 244     else
 245     {
 246         $str =~ s/N/ACGU/g;
 247         $str =~ s/B/CGU/g;
 248         $str =~ s/D/AGU/g;
 249         $str =~ s/H/ACU/g;
 250         $str =~ s/V/ACG/g;
 251         $str =~ s/K/GU/g;
 252         $str =~ s/Y/CU/g;
 253         $str =~ s/S/CG/g;
 254         $str =~ s/W/AU/g;
 255         $str =~ s/R/AG/g;
 256         $str =~ s/M/AC/g;
 257     }
 258
 259     @nts = split //, $str;
 260
 261     %nt_hash = map { $_ => 1 } @nts;
 262
 263     @nts = sort keys %nt_hash;
 264
 265     $key = join "", @nts;
 266
 267     %hash = (
 268         'A'    => 'A',
 269         'C'    => 'C',
 270         'G'    => 'G',
 271         'T'    => 'T',
 272         'U'    => 'U',
 273         'AC'   => 'M',
 274         'AG'   => 'R',
 275         'AT'   => 'W',
 276         'AU'   => 'W',
 277         'CG'   => 'S',
 278         'CT'   => 'Y',
 279         'CU'   => 'Y',
 280         'GT'   => 'K',
 281         'GU'   => 'K',
 282         'ACG'  => 'V',
 283         'ACT'  => 'H',
 284         'ACU'  => 'H',
 285         'AGT'  => 'D',
 286         'AGU'  => 'D',
 287         'CGT'  => 'B',
 288         'CGU'  => 'B',
 289         'ACGT' => 'N',
 290         'ACGU' => 'N',
 291     );
 292
 293     $code = $hash{ $key };
 294
 295     warn qq(WARNING: No ambiguity code for key->$key\n) if not $code;
 296
 297     return $code;
 298 }
 299
 300
 301 sub aa2codons
 302 {
 303     # Martin A. Hansen, March 2005.
 304
 305     # given an amino acid, returns a list of corresponding codons
 306
 307     my ( $aa,   # amino acid to translate
 308        ) = @_;
 309
 310     # returns list
 311
 312     my ( %hash, $codons );
 313
 314     $aa = uc $aa;
 315
 316     %hash = (
 317         'F' => [ 'TTT', 'TTC' ],                               # Phe
 318         'L' => [ 'TTA', 'TTG', 'CTT', 'CTC', 'CTA', 'CTG' ],   # Leu
 319         'S' => [ 'TCT', 'TCC', 'TCA', 'TCG', 'AGT', 'AGC' ],   # Ser
 320         'Y' => [ 'TAT', 'TAC' ],                               # Tyr
 321         '*' => [ 'TAA', 'TAG', 'TGA' ],                        # Stop
 322         'X' => [ 'TAA', 'TAG', 'TGA' ],                        # Stop
 323         'C' => [ 'TGT', 'TGC' ],                               # Cys
 324         'W' => [ 'TGG' ],                                      # Trp
 325         'P' => [ 'CCT', 'CCC', 'CCA', 'CCG' ],                 # Pro
 326         'H' => [ 'CAT', 'CAC' ],                               # His
 327         'Q' => [ 'CAA', 'CAG' ],                               # Gln
 328         'R' => [ 'CGT', 'CGC', 'CGA', 'CGG', 'AGA', 'AGG' ],   # Arg
 329         'I' => [ 'ATT', 'ATC', 'ATA' ],                        # Ile
 330         'M' => [ 'ATG' ],                                      # Met
 331         'T' => [ 'ACT', 'ACC', 'ACA', 'ACG' ],                 # Thr
 332         'N' => [ 'AAT', 'AAC' ],                               # Asn
 333         'K' => [ 'AAA', 'AAG' ],                               # Lys
 334         'V' => [ 'GTT', 'GTC', 'GTA', 'GTG' ],                 # Val
 335         'A' => [ 'GCT', 'GCC', 'GCA', 'GCG' ],                 # Ala
 336         'D' => [ 'GAT', 'GAC' ],                               # Asp
 337         'E' => [ 'GAA', 'GAG' ],                               # Glu
 338         'G' => [ 'GGT', 'GGC', 'GGA', 'GGG' ],                 # Gly
 339     );
 340
 341     $codons = $hash{ $aa };
 342
 343     return wantarray ? @{ $codons } : $codons;
 344 }
 345
 346
 347 sub codon2aa
 348 {
 349     # Martin A. Hansen, March 2005.
 350
 351     # given a codon, returns the correponding
 352     # vertebrate amino acid.
 353
 354     my ( $codon,   # codon to translate
 355        ) = @_;
 356
 357     # returns string
 358
 359     my ( %hash, $aa );
 360
 361     die qq(ERROR: Bad codon: "$codon"\n) if not $codon =~ /[ATCGatcg]{3}/;
 362
 363     %hash = (
 364         'TTT' => 'F',   # Phe
 365         'TTC' => 'F',   # Phe
 366         'TTA' => 'L',   # Leu
 367         'TTG' => 'L',   # Leu
 368         'TCT' => 'S',   # Ser
 369         'TCC' => 'S',   # Ser
 370         'TCA' => 'S',   # Ser
 371         'TCG' => 'S',   # Ser
 372         'TAT' => 'Y',   # Tyr
 373         'TAC' => 'Y',   # Tyr
 374         'TAA' => '*',   # Stop
 375         'TAG' => '*',   # Stop
 376         'TGT' => 'C',   # Cys
 377         'TGC' => 'C',   # Cys
 378         'TGA' => '*',   # Stop
 379         'TGG' => 'W',   # Trp
 380         'CTT' => 'L',   # Leu
 381         'CTC' => 'L',   # Leu
 382         'CTA' => 'L',   # Leu
 383         'CTG' => 'L',   # Leu
 384         'CCT' => 'P',   # Pro
 385         'CCC' => 'P',   # Pro
 386         'CCA' => 'P',   # Pro
 387         'CCG' => 'P',   # Pro
 388         'CAT' => 'H',   # His
 389         'CAC' => 'H',   # His
 390         'CAA' => 'Q',   # Gln
 391         'CAG' => 'Q',   # Gln
 392         'CGT' => 'R',   # Arg
 393         'CGC' => 'R',   # Arg
 394         'CGA' => 'R',   # Arg
 395         'CGG' => 'R',   # Arg
 396         'ATT' => 'I',   # Ile
 397         'ATC' => 'I',   # Ile
 398         'ATA' => 'I',   # Ile
 399         'ATG' => 'M',   # Met
 400         'ACT' => 'T',   # Thr
 401         'ACC' => 'T',   # Thr
 402         'ACA' => 'T',   # Thr
 403         'ACG' => 'T',   # Thr
 404         'AAT' => 'N',   # Asn
 405         'AAC' => 'N',   # Asn
 406         'AAA' => 'K',   # Lys
 407         'AAG' => 'K',   # Lys
 408         'AGT' => 'S',   # Ser
 409         'AGC' => 'S',   # Ser
 410         'AGA' => 'R',   # Arg
 411         'AGG' => 'R',   # Arg
 412         'GTT' => 'V',   # Val
 413         'GTC' => 'V',   # Val
 414         'GTA' => 'V',   # Val
 415         'GTG' => 'V',   # Val
 416         'GCT' => 'A',   # Ala
 417         'GCC' => 'A',   # Ala
 418         'GCA' => 'A',   # Ala
 419         'GCG' => 'A',   # Ala
 420         'GAT' => 'D',   # Asp
 421         'GAC' => 'D',   # Asp
 422         'GAA' => 'E',   # Glu
 423         'GAG' => 'E',   # Glu
 424         'GGT' => 'G',   # Gly
 425         'GGC' => 'G',   # Gly
 426         'GGA' => 'G',   # Gly
 427         'GGG' => 'G',   # Gly
 428     );
 429
 430     $aa = $hash{ uc $codon };
 431
 432     return $aa;
 433 }
 434
 435
 436 sub translate
 437 {
 438     # Martin A. Hansen, June 2005.
 439
 440     # translates a dna sequence to protein according to a optional given
 441     # frame.
 442
 443     my ( $dna,     # dna sequence
 444          $frame,   # frame of translation - OPTIONAL
 445        ) = @_;
 446
 447     # returns string
 448
 449     my ( $codon, $pos, $pep );
 450
 451     $frame ||= 1;
 452
 453     if ( $frame =~ /-?[1-3]/ )
 454     {
 455         if ( $frame < 0 ) {
 456             $dna = Maasha::Seq::dna_revcomp( $dna );
 457         }
 458
 459         $frame = abs( $frame ) - 1;
 460
 461         $dna =~ s/^.{${frame}}//;
 462     }
 463     else
 464     {
 465         Maasha::Common::error( qq(Badly formated frame "$frame") );
 466     }
 467
 468     $pos = 0;
 469
 470     while ( $codon = substr $dna, $pos, 3 )
 471     {
 472         last if not length $codon == 3;
 473
 474         $pep .= codon2aa( $codon );
 475
 476         $pos += 3;
 477     }
 478
 479     return $pep;
 480 }
 481
 482
 483 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> RNA FOLDING <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
 484
 485
 486 sub fold_struct_rnafold
 487 {
 488     # Martin A. Hansen, February 2008.
 489
 490     # Given a squence fold this using RNAfold.
 491
 492     my ( $seq,   # sequence to fold
 493        ) = @_;
 494
 495     # Returns a tuple of fold string and free energy.
 496
 497     my ( $pid, $fh_out, $fh_in, @lines, $struct, $energy );
 498
 499     $pid = open2( $fh_out, $fh_in, "RNAfold -noPS" );
 500
 501     Maasha::Fasta::put_entry( [ "RNAfold", $seq ], $fh_in );
 502
 503     close $fh_in;
 504
 505     @lines = <$fh_out>;
 506
 507     close $fh_out;
 508
 509     waitpid $pid, 0;
 510
 511     chomp @lines;
 512
 513     if ( $lines[ - 1 ] =~ /^([^ ]+) \((.+)\)$/ )
 514     {
 515         $struct = $1;
 516         $energy = $2;
 517     }
 518
 519     return wantarray ? ( $struct, $energy ) : [ $struct, $energy ];
 520 }
 521
 522
 523 sub fold_struct_contrastruct
 524 {
 525     # Martin A. Hansen, February 2008.
 526
 527     # Given a sequence fold this using Contrafold.
 528
 529     my ( $seq,       # sequence to fold
 530          $tmp_dir,   # temporary directory - OPTIONAL
 531        ) = @_;
 532
 533     # Returns a tuple of fold string and temp index.
 534
 535     my ( $tmp_file, $out_file1, $out_file2, $fh, $line, $struct, @AoA, $i, $temp, $index );
 536
 537     $tmp_dir ||= $ENV{ 'BP_TMP' };
 538
 539     $tmp_file  = "$tmp_dir/fold.fna";
 540     $out_file1 = "$tmp_dir/fold.out1";
 541     $out_file2 = "$tmp_dir/fold.out2";
 542
 543     Maasha::Fasta::put_entries( [ [ "fold", $seq ] ], $tmp_file );
 544
 545     Maasha::Common::run( "contrafold", "predict --parens $out_file1 --bpseq $out_file2 $tmp_file" );
 546
 547     unlink $tmp_file;
 548
 549     $fh = Maasha::Common::read_open( $out_file1 );
 550
 551     while ( $line = <$fh> )
 552     {
 553         chomp $line;
 554
 555         $struct = $line;
 556     }
 557
 558     close $fh;
 559
 560     unlink $out_file1;
 561
 562     $fh = Maasha::Common::read_open( $out_file2 );
 563
 564     while ( $line = <$fh> )
 565     {
 566         chomp $line;
 567
 568         push @AoA, [ split " ", $line ];
 569     }
 570
 571     close $fh;
 572
 573     unlink $out_file2;
 574
 575     for ( $i = 0; $i < @AoA; $i++ )
 576     {
 577         if ( $AoA[ $i ]->[ 2 ] != 0 )
 578         {
 579             last if $AoA[ $i ]->[ 0 ] > $AoA[ $i ]->[ 2 ];
 580
 581             $temp += base_pair_melting_temp( $AoA[ $i ]->[ 1 ] . $AoA[ $AoA[ $i ]->[ 2 ] - 1 ]->[ 1 ] );
 582         }
 583     }
 584
 585     $index = sprintf( "%.2f", $temp / length $seq );
 586
 587     return wantarray ? ( $struct, $index ) : [ $struct, $index ];
 588 }
 589
 590
 591 sub base_pair_melting_temp
 592 {
 593     # Martin A. Hansen, February 2008.
 594
 595     # Given a basepair, returns the melting temperature.
 596
 597     my ( $bp,   # basepair string
 598        ) = @_;
 599
 600     # Returns integer
 601
 602     my ( %melt_hash );
 603
 604     %melt_hash = (
 605         AA => 0,
 606         AT => 2,
 607         AC => 0,
 608         AG => 0,
 609         AU => 2,
 610         TA => 2,
 611         TT => 0,
 612         TC => 0,
 613         TG => 1, #
 614         TU => 0,
 615         CA => 0,
 616         CT => 0,
 617         CC => 0,
 618         CG => 4,
 619         CU => 0,
 620         GA => 0,
 621         GT => 1, #
 622         GC => 4,
 623         GG => 0,
 624         GU => 1, #
 625         UA => 2,
 626         UT => 0,
 627         UC => 0,
 628         UG => 1, #
 629         UU => 0,
 630     );
 631
 632     return $melt_hash{ uc $bp };
 633 }
 634
 635
 636 sub generate_dna_oligos
 637 {
 638     # Martin A. Hansen, April 2007.
 639
 640     # Generates all possible DNA oligos of a given wordsize.
 641
 642     # alternative way: perl -MData::Dumper -e '@CONV = glob( "{A,T,C,G}" x 4 ); print Dumper( \@CONV )'
 643
 644
 645     my ( $wordsize,   # size of DNA oligos
 646        ) = @_;
 647
 648     # Returns list
 649
 650     my ( @alph, @oligos, $oligo, $char, @list );
 651
 652     @alph   = ( qw( A T C G N ) );
 653     @oligos = ( '' );
 654
 655     for ( 1 .. $wordsize )
 656     {
 657         foreach $oligo ( @oligos )
 658         {
 659             foreach $char ( @alph ) {
 660                 push @list, $oligo . $char;
 661             }
 662         }
 663
 664         @oligos = @list;
 665
 666         undef @list;
 667     }
 668
 669     return wantarray ? @oligos : \@oligos;
 670 }
 671
 672
 673 sub seq2oligos
 674 {
 675     # Martin A. Hansen, April 2007
 676
 677     # Given a sequence and a wordsize,
 678     # breaks the sequence into overlapping
 679     # oligos of that wordsize.
 680
 681     my ( $seq,       # sequence reference
 682          $wordsize,  # wordsize
 683        ) = @_;
 684
 685     # returns list
 686
 687     my ( $i, $oligo, @oligos );
 688
 689     for ( $i = 0; $i < length( ${ $seq } ) - $wordsize + 1; $i++ )
 690     {
 691         $oligo = substr ${ $seq }, $i, $wordsize;
 692
 693         push @oligos, $oligo;
 694     }
 695
 696     return wantarray ? @oligos : \@oligos;
 697 }
 698
 699
 700 sub seq2oligos_uniq
 701 {
 702     # Martin A. Hansen, April 2007
 703
 704     # Given a sequence and a wordsize,
 705     # breaks the sequence into overlapping
 706     # oligos of that wordsize and return
 707     # only unique words.
 708
 709     my ( $seq,       # sequence reference
 710          $wordsize,  # wordsize
 711        ) = @_;
 712
 713     # returns list
 714
 715     my ( $i, $oligo, %lookup, @oligos );
 716
 717     for ( $i = 0; $i < length( ${ $seq } ) - $wordsize + 1; $i++ )
 718     {
 719         $oligo = substr ${ $seq }, $i, $wordsize;
 720
 721         if ( not exists $lookup{ $oligo } )
 722         {
 723             push @oligos, $oligo;
 724             $lookup{ $oligo } = 1;
 725         }
 726     }
 727
 728     return wantarray ? @oligos : \@oligos;
 729 }
 730
 731
 732 sub oligo_freq
 733 {
 734     # Martin A. Hansen, August 2007.
 735
 736     # Given a hashref with oligo=>count, calculates
 737     # a frequency table. Returns a list of hashes
 738
 739     my ( $oligo_freq,   # hashref
 740        ) = @_;
 741
 742     # Returns data structure
 743
 744     my ( @freq_table, $total );
 745
 746     $total = 0;
 747
 748     map { push @freq_table, { OLIGO => $_, COUNT => $oligo_freq->{ $_ } }; $total += $oligo_freq->{ $_ } } keys %{ $oligo_freq };
 749
 750     @freq_table = sort { $b->{ "COUNT" } <=> $a->{ "COUNT" } or $a->{ "OLIGO" } cmp $b->{ "OLIGO" } } @freq_table;
 751
 752     map { $_->{ "FREQ" } = sprintf( "%.4f", $_->{ "COUNT" } / $total ) } @freq_table;
 753
 754     return wantarray ? return @freq_table : \@freq_table;
 755 }
 756
 757
 758 sub seq_generate
 759 {
 760     # Martin A. Hansen, May 2007
 761
 762     # Generates a random sequence given a sequence length
 763     # and a alphabet.
 764
 765     my ( $len,    # sequence length
 766          $alph,   # sequence alphabet
 767        ) = @_;
 768
 769     # Returns a string
 770
 771     my ( $alph_len, $i, $seq );
 772
 773     $alph_len = scalar @{ $alph };
 774
 775     for ( $i = 0; $i < $len; $i++ ) {
 776         $seq .= $alph->[ int( rand( $alph_len ) ) ];
 777     }
 778
 779     return $seq;
 780 }
 781
 782
 783 sub seq_insert
 784 {
 785     # Martin A. Hansen, June 2009.
 786
 787     # Randomly duplicates a given number of residues in a given sequence.
 788
 789     my ( $seq,          # sequence to mutate
 790          $insertions,   # number of residues to insert
 791        ) = @_;
 792
 793     # Returns a string.
 794
 795     my ( $i, $pos );
 796
 797     for ( $i = 0; $i < $insertions; $i++ )
 798     {
 799         $pos = int( rand( length $seq ) );
 800
 801         substr $seq, $pos, 0, substr $seq , $pos, 1;
 802     }
 803
 804     return $seq;
 805 }
 806
 807
 808 sub seq_delete
 809 {
 810     # Martin A. Hansen, June 2009.
 811
 812     # Randomly deletes a given number of residues from a given sequence.
 813
 814     my ( $seq,         # sequence to mutate
 815          $deletions,   # number of residues to delete
 816        ) = @_;
 817
 818     # Returns a string.
 819
 820     my ( $i, $pos );
 821
 822     for ( $i = 0; $i < $deletions; $i++ )
 823     {
 824         $pos = int( rand( length $seq ) );
 825
 826         substr $seq, $pos, 1, '';
 827     }
 828
 829     return $seq;
 830 }
 831
 832
 833 sub seq_mutate
 834 {
 835     # Martin A. Hansen, June 2009.
 836
 837     # Introduces a given number of random mutations in a
 838     # given sequence of a specified alphabet.
 839
 840     my ( $seq,         # sequence to mutate
 841          $mutations,   # number of mutations
 842          $alph,        # alphabet of sequence
 843        ) = @_;
 844
 845     # Returns a string.
 846
 847     my ( $i, $pos, %lookup_hash );
 848
 849     $i = 0;
 850
 851     while ( $i < $mutations )
 852     {
 853         $pos = int( rand( length $seq ) );
 854
 855         if ( not exists $lookup_hash{ $pos } )
 856         {
 857             substr $seq, $pos, 1, res_mutate( substr( $seq , $pos, 1 ), $alph );
 858
 859             $lookup_hash{ $pos } = 1;
 860
 861             $i++;
 862         }
 863     }
 864
 865     return $seq;
 866 }
 867
 868
 869 sub res_mutate
 870 {
 871     # Martin A. Hansen, June 2009.
 872
 873     # Mutates a given residue to another from a given alphabet.
 874
 875     my ( $res,    # residue to mutate
 876          $alph,   # alphabet
 877        ) = @_;
 878
 879     # Returns a char.
 880
 881     my ( $alph_len, $new );
 882
 883     $alph_len = scalar @{ $alph };
 884     $new      = $res;
 885
 886     while ( uc $new eq uc $res ) {
 887         $new = $alph->[ int( rand( $alph_len ) ) ];
 888     }
 889
 890     return POSIX::islower( $res ) ? lc $new : uc $new;
 891 }
 892
 893
 894 sub seq_shuffle
 895 {
 896     # Martin A. Hansen, December 2007.
 897
 898     # Shuffles sequence of a given string.
 899
 900     my ( $seq,   # sequence string
 901        ) = @_;
 902
 903     # Returns a string.
 904
 905     my ( @list );
 906
 907     @list = split "", $seq;
 908
 909     return join "", shuffle( @list );
 910 }
 911
 912
 913 sub seq_alph
 914 {
 915     # Martin A. Hansen, May 2007.
 916
 917     # returns a requested alphabet
 918
 919     my ( $type,   # alphabet type
 920        ) = @_;
 921
 922     # returns list
 923
 924     my ( @alph );
 925
 926     if ( $type =~ /^dna$/i ) {
 927         @alph = qw( A T C G );
 928     } elsif ( $type =~ /^rna$/i ) {
 929         @alph = qw( A U C G );
 930     } elsif ( $type =~ /^prot/i ) {
 931         @alph = qw( F L S Y C W P H Q R I M T N K V A D E G );
 932     } else {
 933         die qq(ERROR: Unknown alphabet type: "$type"\n);
 934     }
 935
 936     return wantarray ? @alph : \@alph;
 937 }
 938
 939
 940 sub seq_analyze
 941 {
 942     # Martin A. Hansen, August 2007.
 943
 944     # Analyses the sequence composition of a given sequence.
 945
 946     my ( $seq,   # sequence to analyze
 947        ) = @_;
 948
 949     # Returns hash
 950
 951     my ( %analysis, @chars, @chars_lc, $char, %char_hash, $gc, $at, $lc, $max, $res_sum, @indels, %indel_hash );
 952
 953     $analysis{ "SEQ_TYPE" } = uc Maasha::Seq::seq_guess_type( $seq );
 954     $analysis{ "SEQ_LEN" }  = length $seq;
 955
 956     @indels = qw( - ~ . _ );
 957
 958     if ( $analysis{ "SEQ_TYPE" } eq "DNA" )
 959     {
 960         @chars    = split //, "AGCUTRYWSMKHDVBNagcutrywsmkhdvbn";
 961         @chars_lc = split //, "agcutrywsmkhdvbn";
 962     }
 963     elsif ( $analysis{ "SEQ_TYPE" } eq "RNA" )
 964     {
 965         @chars    = split //, "AGCUTRYWSMKHDVBNagcutrywsmkhdvbn";
 966         @chars_lc = split //, "agcutrywsmkhdvbn";
 967     }
 968     else
 969     {
 970         @chars    = split //, "FLSYCWPHQRIMTNKVADEGflsycwphqrimtnkvadeg";
 971         @chars_lc = split //, "flsycwphqrimtnkvadeg";
 972     }
 973
 974     @char_hash{ @chars }   = map { eval "scalar \$seq =~ tr/$_//" } @chars;
 975     @indel_hash{ @indels } = map { eval "scalar \$seq =~ tr/$_//" } @indels;
 976
 977     if ( $analysis{ "SEQ_TYPE" } =~ /DNA|RNA/ )
 978     {
 979         $gc = $char_hash{ "g" } + $char_hash{ "G" } + $char_hash{ "c" } + $char_hash{ "C" };
 980         $at = $char_hash{ "a" } + $char_hash{ "A" } + $char_hash{ "t" } + $char_hash{ "T" } + $char_hash{ "u" } + $char_hash{ "U" };
 981
 982         $analysis{ "GC%" } = sprintf( "%.2f", 100 * $gc / $analysis{ "SEQ_LEN" } );
 983
 984         map { $lc += $char_hash{ lc $_ } } @chars_lc;
 985
 986         $analysis{ "SOFT_MASK%" } = sprintf( "%.2f", 100 * $lc / $analysis{ "SEQ_LEN" } );
 987         $analysis{ "HARD_MASK%" } = sprintf( "%.2f", 100 * ( $char_hash{ "n" } + $char_hash{ "N" } ) / $analysis{ "SEQ_LEN" } );
 988     }
 989
 990     $max = 0;
 991
 992     foreach $char ( @chars_lc )
 993     {
 994         $char = uc $char;
 995
 996         $char_hash{ $char } += $char_hash{ lc $char };
 997
 998         $analysis{ "RES[$char]" } = $char_hash{ $char };
 999
1000         $max = $char_hash{ $char } if $char_hash{ $char } > $max;
1001
1002         $analysis{ "RES_SUM" } += $char_hash{ $char };
1003     }
1004
1005     map { $analysis{ "RES[$_]" } = $indel_hash{ $_ } } @indels;
1006
1007     $analysis{ "MIX_INDEX" } = sprintf( "%.2f", $max / $analysis{ "SEQ_LEN" } );
1008     $analysis{ "MELT_TEMP" } = sprintf( "%.2f", 4 * $gc + 2 * $at );
1009
1010     return wantarray ? %analysis : \%analysis;
1011 }
1012
1013
1014 sub seq_complexity
1015 {
1016     # Martin A. Hansen, May 2008.
1017
1018     # Given a sequence computes a complexity index
1019     # as the most common di-residue over
1020     # the sequence length. Return ~1 if the entire
1021     # sequence is homopolymeric. Above 0.4 indicates
1022     # low complexity sequence.
1023
1024     my ( $seq,   # sequence
1025        ) = @_;
1026
1027     # Returns float.
1028
1029     my ( $len, $i, $max, $di, %hash );
1030
1031     $seq = uc $seq;
1032     $len = length $seq;
1033     $max = 0;
1034
1035     for ( $i = 0; $i < $len - 1; $i++ ) {
1036         $hash{ substr $seq, $i, 2  }++;
1037     }
1038
1039     foreach $di ( keys %hash ) {
1040         $max = $hash{ $di } if $hash{ $di } > $max;
1041     }
1042
1043     return $max / $len;
1044 }
1045
1046
1047 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> SEQLOGO <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1048
1049
1050 sub seqlogo_calc
1051 {
1052     # Martin A. Hansen, January 2007.
1053
1054     # given max bit size and a list of aligned entries
1055     # in FASTA format, calculates for each sequence position
1056     # the height of the letters in bits.
1057     # returns a data structure with [ letter, height ] tuples
1058     # for all letters at each position.
1059
1060     my ( $bit_max,   # maximum bit height
1061          $entries,   # FASTA entries
1062        ) = @_;
1063
1064     # returns data structure
1065
1066     my ( $logo_len, $char_tot, $i, %char_hash, $bit_height, $bit_diff, $char_heights, @logo );
1067
1068     $logo_len = length $entries->[ 0 ]->[ 1 ];
1069     $char_tot = scalar @{ $entries };
1070
1071     for ( $i = 0; $i < $logo_len; $i++ )
1072     {
1073         undef %char_hash;
1074
1075         map { $char_hash{ uc substr( $_->[ 1 ], $i, 1 ) }++ } @{ $entries };
1076
1077         delete $char_hash{ "-" };
1078         delete $char_hash{ "_" };
1079         delete $char_hash{ "~" };
1080         delete $char_hash{ "." };
1081
1082         $bit_height = seqlogo_calc_bit_height( \%char_hash, $char_tot );
1083
1084         $bit_diff = $bit_max - $bit_height;
1085
1086         $char_heights = seqlogo_calc_char_heights( \%char_hash, $char_tot, $bit_diff );
1087
1088         push @logo, $char_heights;
1089     }
1090
1091     return wantarray ? @logo : \@logo;
1092 }
1093
1094
1095 sub seqlogo_calc_bit_height
1096 {
1097     # Martin A. Hansen, January 2007.
1098
1099     # calculates the bit height using Shannon's famous
1100     # general formula for uncertainty as documentet:
1101     # http://www.ccrnp.ncifcrf.gov/~toms/paper/hawaii/latex/node5.html
1102
1103     my ( $char_hash,   # hashref with chars and frequencies
1104          $tot,         # total number of chars
1105        ) = @_;
1106
1107     # returns float
1108
1109     my ( $char, $freq, $bit_height );
1110
1111     foreach $char ( keys %{ $char_hash } )
1112     {
1113         $freq = $char_hash->{ $char } / $tot;
1114
1115         $bit_height += $freq * ( log( $freq ) / log( 2 ) );
1116     }
1117
1118     $bit_height *= -1;
1119
1120     return $bit_height;
1121 }
1122
1123
1124 sub seqlogo_calc_char_heights
1125 {
1126     # Martin A. Hansen, January 2007.
1127
1128     # calculates the hight of each char in bits, and sorts
1129     # according to height.
1130
1131     my ( $char_hash,   # hashref with chars and frequencies
1132          $tot,         # tot number of chars
1133          $bit_diff,    # information gained from uncertainties
1134        ) = @_;
1135
1136     # returns list of tuples
1137
1138     my ( $char, $freq, $char_height, @char_heights );
1139
1140     foreach $char ( keys %{ $char_hash } )
1141     {
1142         $freq = $char_hash->{ $char } / $tot;
1143
1144         $char_height = $freq * $bit_diff;   # char height in bits
1145
1146         push @char_heights, [ $char, $char_height ];
1147     }
1148
1149     @char_heights = sort { $a->[ 1 ] <=> $b->[ 1 ] } @char_heights;
1150
1151     return wantarray ? @char_heights : \@char_heights;
1152 }
1153
1154
1155 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> RESIDUE COLORS <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1156
1157
1158 sub color_pep
1159 {
1160     # Martin A. Hansen, October 2005.
1161
1162     # color scheme for proteins as defined in Mview.
1163     # given a char returns the appropriate color.
1164     # The amino acids are colored according to physicochemical properties:
1165     # bright green = hydrophobic; dark green = large hydrophobic;
1166     # bright blue = negative charge; red = positive charge;
1167     # dull blue = small alcohol; purple = polar; yellow = cysteine.
1168
1169     my ( $char,   # char to decorate
1170        ) = @_;
1171
1172     # returns string
1173
1174     my ( %hash, $color_set );
1175
1176     %hash = (
1177         K   => "bright-red",
1178         R   => "bright-red",
1179         H   => "dark-green",
1180         D   => "bright-blue",
1181         E   => "bright-blue",
1182         S   => "dull-blue",
1183         T   => "dull-blue",
1184         N   => "purple",
1185         Q   => "purple",
1186         A   => "bright-green",
1187         V   => "bright-green",
1188         I   => "bright-green",
1189         L   => "bright-green",
1190         M   => "bright-green",
1191         F   => "dark-green",
1192         Y   => "dark-green",
1193         W   => "dark-green",
1194         C   => "yellow",
1195         G   => "bright-green",
1196         P   => "bright-green",
1197         Z   => "dark-gray",
1198         B   => "dark-gray",
1199         "?" => "light-gray",
1200         "~" => "light-gray",
1201         "*" => "dark-gray",
1202     );
1203
1204     if ( exists $hash{ uc $char } ) {
1205         $color_set = $hash{ uc $char };
1206     } else {
1207         $color_set = "black";
1208     }
1209
1210     return $color_set;
1211 }
1212
1213
1214 sub color_nuc
1215 {
1216     # Martin A. Hansen, October 2005.
1217
1218     # color scheme for nucleotides as defined in Mview.
1219     # given a char returns the appropriate color
1220     # according to physical/chemical proterties.
1221
1222     my ( $char,   # char to decorate
1223        ) = @_;
1224
1225     # returns string
1226
1227     my ( %hash, $color_set );
1228
1229     %hash = (
1230         A => "bright-red",
1231         G => "yellow",
1232         C => "blue",
1233         T => "green",
1234         U => "green",
1235     );
1236
1237     if ( exists $hash{ uc $char } ) {
1238         $color_set = $hash{ uc $char };
1239     } else {
1240         $color_set = "black";
1241     }
1242
1243     return $color_set;
1244 }
1245
1246
1247 sub color_palette
1248 {
1249     # Martin A. Hansen, October 2005.
1250
1251     # hash table with color-names and color-hex.
1252
1253     my ( $color,   # common color name
1254        ) = @_;
1255
1256     # returns string
1257
1258     my ( %hash );
1259
1260     %hash = (
1261         "black"                => "#000000",
1262         "white"                => "#ffffff",
1263         "red"                  => "#ff0000",
1264         "green"                => "#00ff00",
1265         "blue"                 => "#0000ff",
1266         "cyan"                 => "#00ffff",
1267         "magenta"              => "#ff00ff",
1268 #        "yellow"               => "#ffff00",
1269         "yellow"               => "#ffc800",
1270         "purple"               => "#6600cc",
1271         "dull-blue"            => "#0099ff",
1272         "dark-green-blue"      => "#33cccc",
1273         "medium-green-blue"    => "#00ffcc",
1274         "bright-blue"          => "#0033ff",
1275         "dark-green"           => "#009900",
1276         "bright-green"         => "#33cc00",
1277         "orange"               => "#ff3333",
1278         "orange-brown"         => "#cc6600",
1279         "bright-red"           => "#cc0000",
1280         "light-gray"           => "#999999",
1281         "dark-gray"            => "#666666",
1282         "gray0"                => "#ffffff",
1283         "gray1"                => "#eeeeee",
1284         "gray2"                => "#dddddd",
1285         "gray3"                => "#cccccc",
1286         "gray4"                => "#bbbbbb",
1287         "gray5"                => "#aaaaaa",
1288         "gray6"                => "#999999",
1289         "gray7"                => "#888888",
1290         "gray8"                => "#777777",
1291         "gray9"                => "#666666",
1292         "gray10"               => "#555555",
1293         "gray11"               => "#444444",
1294         "gray12"               => "#333333",
1295         "gray13"               => "#222222",
1296         "gray14"               => "#111111",
1297         "gray15"               => "#000000",
1298         "clustal-red"          => "#ff1111",
1299         "clustal-blue"         => "#1155ff",
1300         "clustal-green"        => "#11dd11",
1301         "clustal-cyan"         => "#11ffff",
1302         "clustal-yellow"       => "#ffff11",
1303         "clustal-orange"       => "#ff7f11",
1304         "clustal-pink"         => "#ff11ff",
1305         "clustal-purple"       => "#6611cc",
1306         "clustal-dull-blue"    => "#197fe5",
1307         "clustal-dark-gray"    => "#666666",
1308         "clustal-light-gray"   => "#999999",
1309         "lin-A"                => "#90fe23",
1310         "lin-R"                => "#fe5e2d",
1311         "lin-N"                => "#2e3d2d",
1312         "lin-D"                => "#00903b",
1313         "lin-C"                => "#004baa",
1314         "lin-Q"                => "#864b00",
1315         "lin-E"                => "#3fa201",
1316         "lin-G"                => "#10fe68",
1317         "lin-H"                => "#b2063b",
1318         "lin-I"                => "#04ced9",
1319         "lin-L"                => "#4972fe",
1320         "lin-K"                => "#c4a100",
1321         "lin-M"                => "#2a84dd",
1322         "lin-F"                => "#a60ade",
1323         "lin-P"                => "#fe61fe",
1324         "lin-S"                => "#f7e847",
1325         "lin-T"                => "#fefeb3",
1326         "lin-W"                => "#4a007f",
1327         "lin-Y"                => "#e903a8",
1328         "lin-V"                => "#5bfdfd",
1329     );
1330
1331     if ( exists $hash{ $color } ) {
1332         return $hash{ $color };
1333     } else {
1334         print STDERR qq(WARNING: color "$color" not found in palette!\n);
1335     }
1336 }
1337
1338
1339 sub color_contrast
1340 {
1341     # Martin A. Hansen, October 2005.
1342
1343     # Hash table with contrast colors to be used for frontground
1344     # text on a given background color.
1345
1346     my ( $color,    # background color
1347        ) = @_;
1348
1349     # returns string
1350
1351     my ( %hash );
1352
1353     %hash = (
1354         "black"                => "white",
1355         "white"                => "black",
1356         "red"                  => "white",
1357         "green"                => "white",
1358         "blue"                 => "white",
1359         "cyan"                 => "white",
1360         "magenta"              => "white",
1361         "yellow"               => "black",
1362         "purple"               => "white",
1363         "dull-blue"            => "white",
1364         "dark-green-blue"      => "white",
1365         "medium-green-blue"    => "white",
1366         "bright-blue"          => "white",
1367         "dark-green"           => "white",
1368         "bright-green"         => "black",
1369         "orange"               => "",
1370         "orange-brown"         => "",
1371         "bright-red"           => "white",
1372         "light-gray"           => "black",
1373         "dark-gray"            => "white",
1374         "gray0"                => "",
1375         "gray1"                => "",
1376         "gray2"                => "",
1377         "gray3"                => "",
1378         "gray4"                => "",
1379         "gray5"                => "",
1380         "gray6"                => "",
1381         "gray7"                => "",
1382         "gray8"                => "",
1383         "gray9"                => "",
1384         "gray10"               => "",
1385         "gray11"               => "",
1386         "gray12"               => "",
1387         "gray13"               => "",
1388         "gray14"               => "",
1389         "gray15"               => "",
1390         "clustal-red"          => "black",
1391         "clustal-blue"         => "black",
1392         "clustal-green"        => "black",
1393         "clustal-cyan"         => "black",
1394         "clustal-yellow"       => "black",
1395         "clustal-orange"       => "black",
1396         "clustal-pink"         => "black",
1397         "clustal-purple"       => "black",
1398         "clustal-dull-blue"    => "black",
1399         "clustal-dark-gray"    => "black",
1400         "clustal-light-gray"   => "black",
1401         "lin-A"                => "",
1402         "lin-R"                => "",
1403         "lin-N"                => "",
1404         "lin-D"                => "",
1405         "lin-C"                => "",
1406         "lin-Q"                => "",
1407         "lin-E"                => "",
1408         "lin-G"                => "",
1409         "lin-H"                => "",
1410         "lin-I"                => "",
1411         "lin-L"                => "",
1412         "lin-K"                => "",
1413         "lin-M"                => "",
1414         "lin-F"                => "",
1415         "lin-P"                => "",
1416         "lin-S"                => "",
1417         "lin-T"                => "",
1418         "lin-W"                => "",
1419         "lin-Y"                => "",
1420         "lin-V"                => "",
1421     );
1422
1423     if ( exists $hash{ $color } ) {
1424         return $hash{ $color };
1425     } else {
1426         print STDERR qq(WARNING: color "$color" not found in palette!\n);
1427     }
1428 }
1429
1430
1431 sub seq_word_pack
1432 {
1433     # Martin A. Hansen, April 2008.
1434
1435     # Packs a sequence word into a binary number.
1436
1437     my ( $word,        # Word to be packed
1438        ) = @_;
1439
1440     # Returns integer.
1441
1442     my ( %hash, $bin, $word_size, $pad );
1443
1444     %hash = (
1445         'A' => '000',
1446         'T' => '001',
1447         'C' => '010',
1448         'G' => '100',
1449         'N' => '011',
1450         '-' => '101',
1451         '.' => '110',
1452         '~' => '111',
1453     );
1454
1455     map { $bin .= pack "B3", $hash{ $_ } } split //, $word;
1456
1457     $word_size = length $word;
1458
1459     $pad = ( 3 * $word_size ) / 8;
1460
1461     if ( $pad =~ /\./ )
1462     {
1463         $pad = ( ( int $pad + 1 ) * 8 ) - 3 * $word_size;
1464
1465         $bin .= pack "B$pad", 0 x $pad;
1466     }
1467
1468     return $bin;
1469 }
1470
1471
1472 sub seq_word_unpack
1473 {
1474     # Martin A. Hansen, April 2008.
1475
1476     # Unpacks a binary sequence word to ASCII.
1477
1478     my ( $bin,         # Binary sequence word
1479          $word_size,   # Size of word
1480        ) = @_;
1481
1482     # Returns string.
1483
1484     my ( %hash, $word );
1485
1486     %hash = (
1487         '000' => 'A',
1488         '001' => 'T',
1489         '010' => 'C',
1490         '100' => 'G',
1491         '011' => 'N',
1492         '101' => '-',
1493         '110' => '.',
1494         '111' => '~',
1495     );
1496
1497     map { $word .= $hash{ $_ } } unpack "(B3)$word_size", $bin;
1498
1499     return $word;
1500 }
1501
1502
1503 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>> ADAPTOR LOCATING <<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<
1504
1505
1506 ###############################    REDUNDANT   ##############################
1507
1508 # these functions have been replaced by index_m and match_m in Common.pm
1509
1510 sub find_adaptor
1511 {
1512     # Martin A. Hansen & Selene Fernandez, August 2008
1513
1514     # Locates an adaptor sequence in a given sequence
1515     # starting from a given offset position allowing a given
1516     # number of mismatches (no indels).
1517
1518     my ( $adaptor,       # list of chars
1519          $tag,           # list of chars
1520          $adaptor_len,   # length of adaptor sequence
1521          $tag_len,       # length of tag sequence
1522          $tag_offset,    # offset position
1523          $max_match,     # number of matches indicating success
1524          $max_mismatch,  # number of mismatches
1525        ) = @_;
1526
1527     # Returns an integer.
1528
1529     my ( $i, $j, $match, $mismatch );
1530
1531     $i = $tag_offset;
1532
1533     while ( $i < $tag_len - ( $max_match + $max_mismatch ) + 1 )
1534     {
1535         $j = 0;
1536
1537         while ( $j < $adaptor_len - ( $max_match + $max_mismatch ) + 1 )
1538         {
1539             return $i if check_diag( $adaptor, $tag, $adaptor_len, $tag_len, $j, $i, $max_match, $max_mismatch );
1540
1541             $j++
1542         }
1543
1544         $i++;
1545     }
1546
1547     return -1;
1548 }
1549
1550
1551 sub check_diag
1552 {
1553     # Martin A. Hansen & Selene Fernandez, August 2008
1554
1555     # Checks the diagonal starting at a given coordinate
1556     # of a search space constituted by an adaptor and a tag sequence.
1557     # Residues in the diagonal are compared between the sequences allowing
1558     # for a given number of mismatches. We terminate when search space is
1559     # exhausted or if max matches or mismatches is reached.
1560
1561     my ( $adaptor,       # list of chars
1562          $tag,           # list of chars
1563          $adaptor_len,   # length of adaptor sequence
1564          $tag_len,       # length of tag sequence
1565          $adaptor_beg,   # adaptor begin coordinate
1566          $tag_beg,       # tag begin coordinate
1567          $max_match,     # number of matches indicating success
1568          $max_mismatch,  # number of mismatches
1569        ) = @_;
1570
1571     # Returns boolean.
1572
1573     my ( $match, $mismatch );
1574
1575     $match    = 0;
1576     $mismatch = 0;
1577
1578     while ( $adaptor_beg <= $adaptor_len and $tag_beg <= $tag_len )
1579     {
1580         if ( $adaptor->[ $adaptor_beg ] eq $tag->[ $tag_beg ] )
1581         {
1582             $match++;
1583
1584             return 1 if $match >= $max_match;
1585         }
1586         else
1587         {
1588             $mismatch++;
1589
1590             return 0 if $mismatch > $max_mismatch;
1591         }
1592
1593         $adaptor_beg++;
1594         $tag_beg++;
1595     }
1596
1597     return 0;
1598 }
1599
1600
1601 # >>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>>><<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<<