view release on metacpan or  search on metacpan

lib/Bio/BPWrapper/SeqManipulations.pm  view on Meta::CPAN

                my @prots = Bio::SeqUtils->translate_3frames($seq);
		foreach (@prots) {
		    my $id = $_->id();
		    $id =~ /^(\S+)-(\d)F$/;
		    my ($oriId, $fm) = ($1, $2);
		    $_->id($oriId . "|+" . ($fm+1));
		    $out->write_seq($_);
		}
        } elsif ($frame == 6) {
                my @prots = Bio::SeqUtils->translate_6frames($seq);
		foreach (@prots) {
		    my $id = $_->id();
		    $id =~ /^(\S+)-(\d)([RF])$/;
		    my ($oriId, $fm, $dir) = ($1, $2, $3);
		    if ($dir eq 'F') {
			$_->id($oriId . "|+" . ($fm+1));
		    } else {
			$_->id($oriId . "|-" . ($fm+1));
		    }
		    $out->write_seq($_);
		}
        } else { warn "Accepted frame arguments: 1, 3, and 6\n"}
    }
}

=head2 restrict_coord()

Note: This function is currently DEPRECATED.

Finds digestion coordinates by a specified restriction enzyme
specified in C<$opts{restrinct}> set via L<C<#initilize(\%opts)>|/initialize>.

An input file with sequences is expected. Wraps
L<Bio::Restriction::Analysis-E<gt>cut()|https://github.com/bioperl/Bio-Restriction>.

Outputs coordinates of overhangs in BED format.


sub restrict_coord {
    use Bio::Restriction::Analysis;
    use Bio::Restriction::EnzymeCollection;

    my $enz = $opts{"restrict-coord"};
    my $re = Bio::Restriction::EnzymeCollection->new()->get_enzyme($enz);
    my $len = length($re->overhang_seq());

    while ( $seq = $in->next_seq() ) {
        my $seq_str = $seq->seq();
        die "Not a DNA sequence\n" unless $seq_str =~ /^[ATCGRYSWKMBDHVN]+$/i;
        my $ra = Bio::Restriction::Analysis->new(-seq=>$seq);
        foreach my $pos ($ra->positions($enz)) {
	    print $seq->id()."\t".($pos-$len)."\t".$pos."\n";
        }
    }
}

=head2 restrict_digest()

Note: This function is currently DEPRECATED.

Predicted fragments from digestion by a specified restriction enzyme
specified in C<$opts{restrinct}> set via L<C<#initilize(\%opts)>|/initialize>.

An input file with sequences is expected. Wraps
L<Bio::Restriction::Analysis-E<gt>cut()|https://metacpan.org/pod/Bio::Restriction::Analysis#cut>.


sub restrict_digest {
    my $enz = $opts{"restrict"};
    use Bio::Restriction::Analysis;
    while ( $seq = $in->next_seq() ) {
	my $seq_str = $seq->seq();
	die "Not a DNA sequence\n" unless $seq_str =~ /^[ATCGRYSWKMBDHVN]+$/i;
	my $ra = Bio::Restriction::Analysis->new(-seq=>$seq);
	foreach my $frag ($ra->fragment_maps($enz)) {
	    my $seq_obj = Bio::Seq->new(
		-id=>$seq->id().'|'.$frag->{start}.'-'.$frag->{end}.'|'.($frag->{end}-$frag->{start}+1),
		-seq=>$frag->{seq});
	    $out->write_seq($seq_obj)
	}
    }
}

=cut

=head2 anonymize()

Replace sequence IDs with serial IDs I<n> characters long, as specified in
C<$opts{'anonymize'}> set via L<C<#initilize(\%opts)>|/initialize>.
For example if C<$opts{'anonymize'}>, the first ID will be C<S0001>.
leading 'S' The length of the serial idea

A sed script file is produced with a F<.sed> suffix that may be used
with sed's C<'-f'> argument. If the filename is F<'-'>, the sed file
is named C<STDOUT.sed> instead. A message containing the sed filename is
written to C<STDERR>.

=cut

sub anonymize {
    my $char_len = $opts{"anonymize"} // die "Tried to use option 'prefix' without using option 'anonymize'. Exiting...\n";
    my $prefix = (defined($opts{"prefix"})) ? $opts{"prefix"} : "S";

    pod2usage(1) if $char_len < 1;

    my $ct = 1;
    my %serial_name;
    my $length_warn = 0;
    while ($seq = $in->next_seq()) {
        my $serial = $prefix . sprintf "%0" . ($char_len - length($prefix)) . "s", $ct;
        $length_warn = 1 if length($serial) > $char_len;
        $serial_name{$serial} = $seq->id();
        $seq->id($serial);
        $out->write_seq($seq);
        $ct++
    }

    _make_sed_file($filename, %serial_name);
    warn "Anonymization map:\n";
    while (my ($k, $v) = each %serial_name) { warn "$k => $v\n" }

    warn "WARNING: Anonymized ID length exceeded requested length: try a different length or prefix.\n" if $length_warn
}

=head2 shred_seq()

Break into individual sequences writing a FASTA file for each sequence.

=cut


sub shred_seq {
    while ($seq = $in->next_seq()) {
        my $newid = $seq->id();
        $newid =~ s/[\s\|]/_/g;
        print $newid, "\n";
	my $suffix = $out_format || "fas";  
        my $newout = Bio::SeqIO->new(-format => $out_format, -file => ">" . $newid . ".$suffix");