[Bio] / FigKernelPackages / gjoseqlib.pm Repository:
ViewVC logotype

Diff of /FigKernelPackages/gjoseqlib.pm

Parent Directory Parent Directory | Revision Log Revision Log | View Patch Patch

revision 1.6, Sun Jun 10 17:24:38 2007 UTC revision 1.9, Thu Nov 1 20:07:42 2007 UTC
# Line 51  Line 51 
51  #  #
52  #  @entry  = subseq_DNA_entry( @seq_entry, $from, $to [, $fix_id] );  #  @entry  = subseq_DNA_entry( @seq_entry, $from, $to [, $fix_id] );
53  #  @entry  = subseq_RNA_entry( @seq_entry, $from, $to [, $fix_id] );  #  @entry  = subseq_RNA_entry( @seq_entry, $from, $to [, $fix_id] );
54    #  $DNAseq = DNA_subseq(  $seq, $from, $to );
55    #  $DNAseq = DNA_subseq( \$seq, $from, $to );
56    #  $RNAseq = RNA_subseq(  $seq, $from, $to );
57    #  $RNAseq = RNA_subseq( \$seq, $from, $to );
58  #  @entry  = complement_DNA_entry( @seq_entry [, $fix_id] );  #  @entry  = complement_DNA_entry( @seq_entry [, $fix_id] );
59  #  @entry  = complement_RNA_entry( @seq_entry [, $fix_id] );  #  @entry  = complement_RNA_entry( @seq_entry [, $fix_id] );
60  #  $DNAseq = complement_DNA_seq( $NA_seq );  #  $DNAseq = complement_DNA_seq( $NA_seq );
# Line 67  Line 71 
71  #  User-supplied genetic code must be upper case index and match the  #  User-supplied genetic code must be upper case index and match the
72  #  DNA versus RNA type of sequence  #  DNA versus RNA type of sequence
73  #  #
74    #  $seq = translate_seq_with_user_code( $seq, $gen_code_hash [, $met_start] )
75    #
76  #  Locations (= oriented intervals) are ( id, start, end )  #  Locations (= oriented intervals) are ( id, start, end )
77  #  Intervals are ( id, left, right )  #  Intervals are ( id, left, right )
78  #  #
# Line 81  Line 87 
87  #  Convert GenBank locations to SEED locations  #  Convert GenBank locations to SEED locations
88  #  #
89  #  @seed_locs = gb_location_2_seed( $contig, @gb_locs )  #  @seed_locs = gb_location_2_seed( $contig, @gb_locs )
90    #
91    #  Read quality scores from a fasta-like file:
92    #
93    #  @seq_entries = read_qual( )               #  STDIN
94    # \@seq_entries = read_qual( )               #  STDIN
95    #  @seq_entries = read_qual( \*FILEHANDLE )
96    # \@seq_entries = read_qual( \*FILEHANDLE )
97    #  @seq_entries = read_qual(  $filename )
98    # \@seq_entries = read_qual(  $filename )
99    #
100    
101  use strict;  use strict;
102    use Carp;
 use gjolib qw( wrap_text );  
103    
104  #  Exported global variables:  #  Exported global variables:
105    
# Line 131  Line 145 
145    
146          subseq_DNA_entry          subseq_DNA_entry
147          subseq_RNA_entry          subseq_RNA_entry
148            DNA_subseq
149            RNA_subseq
150          complement_DNA_entry          complement_DNA_entry
151          complement_RNA_entry          complement_RNA_entry
152          complement_DNA_seq          complement_DNA_seq
# Line 152  Line 168 
168          reverse_intervals          reverse_intervals
169    
170          gb_location_2_seed          gb_location_2_seed
171    
172            read_qual
173          );          );
174    
175  our @EXPORT_OK = qw(  our @EXPORT_OK = qw(
# Line 430  Line 448 
448      my ( $fh, undef, $close, $unused ) = output_filehandle( shift );      my ( $fh, undef, $close, $unused ) = output_filehandle( shift );
449      ( unshift @_, $unused ) if $unused;      ( unshift @_, $unused ) if $unused;
450    
451      ( ref( $_[0] ) eq "ARRAY" ) or die "Bad sequence entry passed to print_alignment_as_fasta\n";      ( ref( $_[0] ) eq "ARRAY" ) or confess "Bad sequence entry passed to print_alignment_as_fasta\n";
452    
453      #  Expand the sequence entry list if necessary:      #  Expand the sequence entry list if necessary:
454    
# Line 633  Line 651 
651    
652    
653  #-----------------------------------------------------------------------------  #-----------------------------------------------------------------------------
654    #  Return a string with text wrapped to defined line lengths:
655    #
656    #     $wrapped_text = wrap_text( $str )                  # default len   =  80
657    #     $wrapped_text = wrap_text( $str, $len )            # default ind   =   0
658    #     $wrapped_text = wrap_text( $str, $len, $indent )   # default ind_n = ind
659    #     $wrapped_text = wrap_text( $str, $len, $indent_1, $indent_n )
660    #-----------------------------------------------------------------------------
661    sub wrap_text {
662        my ($str, $len, $ind, $indn) = @_;
663    
664        defined($str)  || die "wrap_text called without a string\n";
665        defined($len)  || ($len  =   80);
666        defined($ind)  || ($ind  =    0);
667        ($ind  < $len) || die "wrap error: indent greater than line length\n";
668        defined($indn) || ($indn = $ind);
669        ($indn < $len) || die "wrap error: indent_n greater than line length\n";
670    
671        $str =~ s/\s+$//;
672        $str =~ s/^\s+//;
673        my ($maxchr, $maxchr1);
674        my (@lines) = ();
675    
676        while ($str) {
677            $maxchr1 = ($maxchr = $len - $ind) - 1;
678            if ($maxchr >= length($str)) {
679                push @lines, (" " x $ind) . $str;
680                last;
681            }
682            elsif ($str =~ /^(.{0,$maxchr1}\S)\s+(\S.*)$/) { # no expr in {}
683                push @lines, (" " x $ind) . $1;
684                $str = $2;
685            }
686            elsif ($str =~ /^(.{0,$maxchr1}-)(.*)$/) {
687                push @lines, (" " x $ind) . $1;
688                $str = $2;
689            }
690            else {
691                push @lines, (" " x $ind) . substr($str, 0, $maxchr);
692                $str = substr($str, $maxchr);
693            }
694            $ind = $indn;
695        }
696    
697        return join("\n", @lines);
698    }
699    
700    
701    #-----------------------------------------------------------------------------
702  #  Build an index from seq_id to pointer to sequence entry: (id, desc, seq)  #  Build an index from seq_id to pointer to sequence entry: (id, desc, seq)
703  #  #
704  #     my \%seq_ind  = index_seq_list(  @seq_list );  #     my \%seq_ind  = index_seq_list(  @seq_list );
# Line 794  Line 860 
860  }  }
861    
862    
863    sub DNA_subseq
864    {
865        my ( $seq, $from, $to ) = @_;
866    
867        my $len = ref( $seq ) eq 'SCALAR' ? length( $$seq )
868                                          : length(  $seq );
869        if ( ( $from eq '$' ) || ( $from eq "" ) ) { $from = $len }
870        if ( ( $to   eq '$' ) || ( ! $to       ) ) { $to   = $len }
871    
872        my $left  = ( $from < $to ) ? $from : $to;
873        my $right = ( $from < $to ) ? $to   : $from;
874        if ( ( $right < 1 ) || ( $left > $len ) ) { return "" }
875        if ( $right > $len ) { $right = $len }
876        if ( $left  < 1    ) { $left  =    1 }
877    
878        my $subseq = ref( $seq ) eq 'SCALAR' ? substr( $$seq, $left-1, $right-$left+1 )
879                                             : substr(  $seq, $left-1, $right-$left+1 );
880    
881        if ( $from > $to )
882        {
883            $subseq = reverse $subseq;
884            $subseq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
885                         [TGCAAMKYSWRVHDBNtgcaamkyswrvhdbn];
886        }
887    
888        $subseq
889    }
890    
891    
892    sub RNA_subseq
893    {
894        my ( $seq, $from, $to ) = @_;
895    
896        my $len = ref( $seq ) eq 'SCALAR' ? length( $$seq )
897                                          : length(  $seq );
898        if ( ( $from eq '$' ) || ( $from eq "" ) ) { $from = $len }
899        if ( ( $to   eq '$' ) || ( ! $to       ) ) { $to   = $len }
900    
901        my $left  = ( $from < $to ) ? $from : $to;
902        my $right = ( $from < $to ) ? $to   : $from;
903        if ( ( $right < 1 ) || ( $left > $len ) ) { return "" }
904        if ( $right > $len ) { $right = $len }
905        if ( $left  < 1    ) { $left  =    1 }
906    
907        my $subseq = ref( $seq ) eq 'SCALAR' ? substr( $$seq, $left-1, $right-$left+1 )
908                                             : substr(  $seq, $left-1, $right-$left+1 );
909    
910        if ( $from > $to )
911        {
912            $subseq = reverse $subseq;
913            $subseq =~ tr[ACGTUKMRSWYBDHVNacgtukmrswybdhvn]
914                         [UGCAAMKYSWRVHDBNugcaamkyswrvhdbn];
915        }
916    
917        $subseq
918    }
919    
920    
921  sub complement_DNA_entry {  sub complement_DNA_entry {
922      my ($id, $desc, $seq, $fix_id) = @_;      my ($id, $desc, $seq, $fix_id) = @_;
923      $fix_id ||= 0;     #  fix undef values      $fix_id ||= 0;     #  fix undef values
# Line 1001  Line 1125 
1125  }  }
1126    
1127    
1128  #  Construct the genetic code with selanocysteine by difference:  #  Construct the genetic code with selenocysteine by difference:
1129    
1130  %genetic_code_with_U = map { $_ => $genetic_code{ $_ } } keys %genetic_code;  %genetic_code_with_U = map { $_ => $genetic_code{ $_ } } keys %genetic_code;
1131  $genetic_code_with_U{ TGA } = "U";  $genetic_code_with_U{ TGA } = "U";
# Line 1676  Line 1800 
1800  }  }
1801    
1802    
1803    #-----------------------------------------------------------------------------
1804    #  Read qual.
1805    #
1806    #  Save the contents in a list of refs to arrays: [ $id, $descript, \@qual ]
1807    #
1808    #     @seq_entries = read_qual( )               #  STDIN
1809    #    \@seq_entries = read_qual( )               #  STDIN
1810    #     @seq_entries = read_qual( \*FILEHANDLE )
1811    #    \@seq_entries = read_qual( \*FILEHANDLE )
1812    #     @seq_entries = read_qual(  $filename )
1813    #    \@seq_entries = read_qual(  $filename )
1814    #-----------------------------------------------------------------------------
1815    sub read_qual {
1816        my ( $fh, $name, $close, $unused ) = input_filehandle( $_[0] );
1817        $unused && die "Bad reference type (" . ref( $unused ) . ") passed to read_qual\n";
1818    
1819        my @quals = ();
1820        my ($id, $desc, $qual) = ("", "", []);
1821    
1822        while ( <$fh> ) {
1823            chomp;
1824            if (/^>\s*(\S+)(\s+(.*))?$/) {        #  new id
1825                if ($id && @$qual) { push @quals, [ $id, $desc, $qual ] }
1826                ($id, $desc, $qual) = ($1, $3 ? $3 : "", []);
1827            }
1828            else {
1829                push @$qual, split;
1830            }
1831        }
1832        close( $fh ) if $close;
1833    
1834        if ($id && @$qual) { push @quals, [ $id, $desc, $qual ] }
1835        return wantarray ? @quals : \@quals;
1836    }
1837    
1838    
1839  1;  1;

Legend:
Removed from v.1.6  
changed lines
  Added in v.1.9

MCS Webmaster
ViewVC Help
Powered by ViewVC 1.0.3