view release on metacpan or  search on metacpan

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

		foreach (@{$ref->{$mono}}) { push @pos, $_} 
	    }
#	    print STDERR join "\t", sort {$a <=> $b} @pos;
#	    print STDERR "\n"; 
	}
    }
    exit;
}

sub _count_gap {
    my ($str, $type) = @_;
    my @mono = split('', $str);
    my %seen_gaps;
    my $num_gaps = 0;
    my $ct = 0;
    foreach my $mon (@mono) {
	$ct++;
	next if ($type eq 'dna' && $mon =~ /[atcg]/i) or ($type eq 'protein' && $mon =~ /[ACDEFGHIKLMNPQRSTVWY]/i) or ($type eq 'protein' && $mon =~ /\*\s*$/);
	$num_gaps ++;
	if ($seen_gaps{$mon}) {
	    push @{$seen_gaps{$mon}}, $ct;
	} else {
	    $seen_gaps{$mon} = [$ct]
	}
    }
#    print Dumper(\%seen_gaps) if $num_gaps;
    return ($num_gaps, \%seen_gaps);
}

sub rename_id {
    my $optStr = $opts{rename};
    my %names;

    if ($optStr =~  /^id:(\S+);(\S+)$/) {
	$names{$1} = $2;
    } else {
	open NAME, "<", $optStr || die "a file with old-tab-new needed\n";
	while(<NAME>) {
	    chomp;
	    my ($oldN, $newN) = split;
	    $names{$oldN} = $newN;
	}
	close NAME;
    }

    while( my $seqobj  = $in->next_seq() ) {
	my $id = $seqobj->display_id();
	if ($names{$id}) {
#	    $seqobj->id($id . "|" . $names{$id});
#	    warn "$id appended by $names{$id}\n";
	    $seqobj->id($names{$id});
	    warn "$id replaced by $names{$id}\n";
	} else {
	    warn "$id not changed\n";
	}
	$out->write_seq($seqobj);
    }
}


sub update_longest_orf {
    while( my $seqobj  = $in->next_seq() ) {
	my $pep_string = $seqobj->translate( undef, undef, 0 )->seq();
	unless ($pep_string =~ /\*[A-Z]/) { # no internal stop; don't proceed
	    my $id = $seqobj->id();
	    $seqobj->id($id . "|+1");
	    $seqobj = &_trim_end_to_frame($seqobj);
	    $out->write_seq($seqobj);
#	    warn $seqobj->id, ": +1 ok\n";
	    next;
	}

	unless ($opts{"no-revcom"}) { # do not search in revcom
	    my $pep_rev = $seqobj->revcom()->translate( undef, undef, 0 )->seq();
	    unless ($pep_rev =~ /\*[A-Z]/) { # no internal stop for revcom
		my $id = $seqobj->id();
		$seqobj->id($id . "|-1");
		my $rev = $seqobj->revcom();
		$rev = &_trim_end_to_frame($rev);
		$out->write_seq($rev);
#	    warn $seqobj->id(), ": -1 ok\n";
		next;
	    }
	}

	my $longest = {
	    'aa_start'  => 1,
	    'aa_end'    => 1,
	    'aa_length' => 1,
	    'nt_start'  => 1,
	    'nt_end'    => 1,
	    'nt_seq'    => $seqobj->seq(),
	    'frame'     => 1,
	};

	foreach my $fm (1, 2, 3, -1, -2, -3 ) {
#	    warn "checking frame $fm ...\n";
	    next if $opts{"no-revcom"} && $fm < 0; # do not search in revcom
	    my $new_seqobj = Bio::Seq->new(
		-id  => $seqobj->id() . "|$fm",
		-seq => $fm > 0 ? $seqobj->subseq( $fm, $seqobj->length() ) : $seqobj->revcom()->subseq( abs($fm), $seqobj->length() )
		);    # chop seq to frame first
	    
	    &_get_longest($new_seqobj, $longest, $fm);
#	    warn "longest ORF:", $longest->{aa_length}, "\n";

	}
#	warn "start codon not M/V/L:", $seqobj->id() unless substr( $longest->{nt_seq}, 0, 3 ) =~ /[atg|gt[atcg]|ct[atcg]|tt[ag]/i;
#	print ">", $seqobj->id, "|f", $longest->{frame}, "|longest-orf\n", $longest->{nt_seq}, "\n";
	my $fid = $longest->{frame} > 0 ? "+" . $longest->{frame} : $longest->{frame};
	my $longest_seq = Bio::Seq->new(-id => $seqobj->id . "|" . $fid, -seq => $longest->{nt_seq} );
	$longest_seq = &_trim_end_to_frame($longest_seq);
	$out->write_seq($longest_seq);
    }
}

sub _trim_end_to_frame {
    my $seqobj = shift;
    my $seq_len = $seqobj->length();
    my $remainder = $seq_len % 3;
    my $seqstr = $seqobj->subseq(1, $seq_len - $remainder);