Bio-AutomatedAnnotation
view release on metacpan or search on metacpan
lib/Bio/AutomatedAnnotation/Prokka.pm view on Meta::CPAN
# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
# read in sequences; remove small contigs; replace ambig with N
my $in = $self->assembly_file;
$self->msg("Loading and checking input file: $in");
my $fin = Bio::SeqIO->new( -file => $in, -format => 'fasta' );
my $fout = Bio::SeqIO->new( -file => ">$outdir/$prefix.fna", -format => 'fasta' );
my $ncontig = 0;
while ( my $seq = $fin->next_seq ) {
if ( $seq->length < $mincontig ) {
$self->msg( "Skipping short (<$mincontig bp) contig:", $seq->display_id );
next;
}
$ncontig++;
my $id;
# http://www.ncbi.nlm.nih.gov/genomes/static/Annotation_pipeline_README.txt
if($self->keep_original_order_and_names)
{
$id = $seq->display_id;
}
else
{
$id = sprintf "contig%06d", $ncontig;
$id = "$contig_uniq_id|$centre|$id" if $centre;
$seq->display_id($id);
}
my $s = $seq->seq;
$s = uc($s);
$s =~ s/[^ACTG]/N/g;
$seq->seq($s);
$seq->desc(undef);
$fout->write_seq($seq);
$seq{$id}{DNA} = $seq;
$seq{$id}{order} = $ncontig;
}
$self->msg("Wrote $ncontig contigs");
#$self->msg(sort keys %seq); exit;
# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
# tRNA + tmRNA
$self->msg("Predicting tRNAs and tmRNAs");
my $cmd = "aragorn -gc$gcode -w $outdir/$prefix.fna"; # -t/-m
$self->msg("Running: $cmd");
my $num_trna = 0;
open TRNA, "$cmd |";
my $sid;
while (<TRNA>) {
chomp;
if (m/^>end/) {
last;
}
if (m/^>(\S+)/) {
$sid = $1;
next;
}
my @x = split m/\s+/;
next unless @x == 5 and $x[0] =~ m/^\d+$/;
# and $x[4] =~ m/^\([ATCG]{3}\)$/i;
#$self->msg($_);
$self->msg("@x");
$x[2] =~ m/(c)?\[(\d+),(\d+)\]/;
my ( $revcom, $start, $end ) = ( $1, $2, $3 );
# bug fix for aragorn when revcom trna ends at start of contig!
# if (defined $revcom and $start > $end) {
# $self->msg("Activating kludge for Aragorn bug for tRNA end at contig start");
# $start = 1;
# }
if ( $start > $end ) {
$self->msg("tRNA $x[2] has start($start) > end ($end) - skipping.");
next;
}
if ( abs( $end - $start ) > 500 ) {
$self->msg("tRNA/tmRNA $x[2] is too big (>500bp) - skipping.");
next;
}
# end kludge
$num_trna++;
my $ftype = 'tRNA';
my $product = $x[1] . $x[4];
my @gene = ();
if ( $x[1] eq 'tmRNA' ) {
$ftype = $x[1];
$product = "transfer-messenger RNA, SsrA";
@gene = ( 'gene' => 'ssrA' );
}
my $tool = "Aragorn:" . $tools{aragorn}->{VERSION};
push @{ $seq{$sid}{FEATURE} }, Bio::SeqFeature::Generic->new(
-primary => $ftype, # tRNA or tmRNA
-seq_id => $sid,
-source => $tool,
-start => $start,
-end => $end,
-strand => ( defined $revcom ? -1 : +1 ),
-score => undef,
-frame => 0,
-tag => {
'product' => $product,
'inference' => "COORDINATES:profile:$tool",
@gene,
}
);
}
$self->msg("Found $num_trna tRNAs");
# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
# rRNA
if ( $kingdom ne 'Viruses' ) {
$self->msg("Predicting Ribosomal RNAs");
my $rnammerfn = "$outdir/rnammer.xml";
my $num_rrna = 0;
lib/Bio/AutomatedAnnotation/Prokka.pm view on Meta::CPAN
# same contig, overlapping (could check same strand too? not sure)
if ( $rna->seq_id eq $sid and $cds->overlaps($rna) ) {
$overlap = $rna;
last;
}
}
if ($overlap) {
$self->msg("Not including CDS which overlaps existing RNA at $sid:$1..$2 on $3 strand");
}
else {
$num_cds++;
push @{ $seq{$sid}{FEATURE} }, $cds;
}
}
}
$self->msg("Found $num_cds CDS");
# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
# Connect features to their parent sequences
$self->msg("Connecting features back to sequences");
for my $sid ( sort keys %seq ) {
for my $f ( @{ $seq{$sid}{FEATURE} } ) {
$f->attach_seq( $seq{$sid}{DNA} );
}
}
# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
# Find signal peptide leader sequences
my $sigpver = substr $tools{signalp}{VERSION}, 0, 1; # first char, expect 3 or 4
if ( $kingdom eq 'Bacteria' and $sigpver == 3 || $sigpver == 4 ) {
if ($gram) {
$gram = $gram =~ m/\+|[posl]/i ? 'gram+' : 'gram-';
$self->msg("Looking for signal peptides at start of predicted proteins");
$self->msg("Treating $kingdom as $gram");
my $spoutfn = "$outdir/signalp.faa";
my $spout = Bio::SeqIO->new( -file => ">$spoutfn", -format => 'fasta' );
my %cds;
my $count = 0;
for my $sid ( sort keys %seq ) {
for my $f ( @{ $seq{$sid}{FEATURE} } ) {
next unless $f->primary_tag eq 'CDS';
$cds{ ++$count } = $f;
my $seq = $f->seq->translate;
$seq->display_id($count);
$spout->write_seq($seq);
}
}
my $opts = $sigpver == 3 ? '-m hmm' : '';
my $cmd = "signalp -t $gram -f short $opts $spoutfn 2> /dev/null";
$self->msg("Running: $cmd");
my $tool = "SignalP:" . $tools{signalp}->{VERSION};
my $num_sigpep = 0;
open SIGNALP, "$cmd |";
while (<SIGNALP>) {
my @x = split m/\s+/;
if ( $sigpver == 3 ) {
next unless @x == 7 and $x[6] eq 'Y'; # has sig_pep
my $parent = $cds{ $x[0] };
my $prob = $x[5];
my $cleave = $x[3];
my $start = $parent->strand > 0 ? $parent->start : $parent->end;
my $end = $start + $parent->strand * ( $cleave - 1 );
my $sigpep = Bio::SeqFeature::Generic->new(
-seq_id => $parent->seq_id,
-source_tag => $tool,
-primary => 'sig_peptide',
-start => min( $start, $end ),
-end => max( $start, $end ),
-strand => $parent->strand,
-frame => 0, # PHASE: compulsory for peptides, can't be '.'
-tag => {
# 'ID' => $ID,
# 'Parent' => $x[0], # don't have proper IDs yet....
'product' => "putative signal peptide",
'inference' => "ab initio prediction:$tool",
'note' => "predicted cleavage at residue $x[3] with probability $prob",
}
);
push @{ $seq{ $parent->seq_id }{FEATURE} }, $sigpep;
$num_sigpep++;
}
else {
# $self->msg("sigp$sigpver: @x");
next unless @x == 12 and $x[9] eq 'Y'; # has sig_pep
my $parent = $cds{ $x[0] };
my $cleave = $x[2];
my $start = $parent->strand > 0 ? $parent->start : $parent->end;
my $end = $start + $parent->strand * ( $cleave - 1 );
my $sigpep = Bio::SeqFeature::Generic->new(
-seq_id => $parent->seq_id,
-source_tag => $tool,
-primary => 'sig_peptide',
-start => min( $start, $end ),
-end => max( $start, $end ),
-strand => $parent->strand,
-frame => 0, # PHASE: compulsory for peptides, can't be '.'
-tag => {
# 'ID' => $ID,
# 'Parent' => $x[0], # don't have proper IDs yet....
'product' => "putative signal peptide",
'inference' => "ab initio prediction:$tool",
'note' => "predicted cleavage at residue $x[2]",
}
);
push @{ $seq{ $parent->seq_id }{FEATURE} }, $sigpep;
$num_sigpep++;
}
}
$self->msg("Found $num_sigpep signal peptides");
$self->delfile($spoutfn);
}
else {
$self->msg("Option --gram not specified, will NOT check for signal peptides.");
lib/Bio/AutomatedAnnotation/Prokka.pm view on Meta::CPAN
$self->msg("Preparing user-supplied primary annotation source: $proteins");
$self->runcmd("makeblastdb -dbtype prot -in '$proteins' -out $outdir/proteins -logfile /dev/null");
unshift @database,
{
DB => "$outdir/proteins",
SRC => 'similar to AA sequence:RefSeq:',
FMT => 'blast',
CMD => $BLASTPCMD,
};
}
if ($fast) {
$self->msg("Option --fast enabled, so skipping CDS similarity searches");
}
else {
$self->msg("Annotating CDS, please be patient.");
my $paropts = $cpus > 0 ? " -j $cpus" : "";
$self->msg( "Will use", ( $cpus > 0 ? $cpus : 'all available' ), "CPUs for similarity searching." );
my $num_cleaned = 0;
my %cds;
my $count = 0;
for my $sid ( sort keys %seq ) {
for my $f ( @{ $seq{$sid}{FEATURE} } ) {
next unless $f->primary_tag eq 'CDS';
next if $f->has_tag('product');
$cds{ ++$count } = $f;
}
}
if ( $count > 0 ) {
#Â Minimise the number of files created at a time. Tradeoff with efficiency of parallelisation.
# This creates X files per CPU.
my $slice_size = ( ( $cpus > 0 ) ? $cpus : 1 ) * ( ( $files_per_chunk <= 0 ) ? 10 : $files_per_chunk );
my @cds_counter = sort( keys %cds );
for ( my $i = 0 ; $i < ceil( (@cds_counter) / $slice_size ) ; $i++ ) {
for ( my $j = $slice_size * $i ; $j < @cds_counter && $j < $slice_size * ( $i + 1 ) ; $j++ ) {
$self->create_cds_sequences_in_file( $tempdir, $cds_counter[$j], $cds{ $cds_counter[$j] } );
}
for my $db (@database) {
my $cmd = $db->{CMD};
$cmd =~ s/%i/{}/g;
$cmd =~ s/%o/{}.out/g;
$cmd =~ s/%e/$evalue/g;
$cmd =~ s,%d,$db->{DB},g;
$self->msg( $db->{FMT}, "$count (of $num_cds) proteins against", $db->{DB} );
$self->runcmd("nice parallel$paropts $cmd ::: $tempdir/*.seq");
for ( my $j = $slice_size * $i ; $j < @cds_counter && $j < $slice_size * ( $i + 1 ) ; $j++ ) {
my $pid = $cds_counter[$j];
my $bls = Bio::SearchIO->new( -file => "$tempdir/$pid.seq.out", -format => $db->{FMT}, -version => $db->{VERSION} );
my $res = $bls->next_result or next;
my $hit = $res->next_hit or next;
my ( $prod, $gene, $EC ) = ( $hit->description, '', '' );
if ( $prod =~ m/~~~/ ) {
( $EC, $gene, $prod ) = split m/~~~/, $prod;
$EC =~ s/n\d+/-/g; # collapse transitionary EC numbers
}
my $cleanprod = $prod;
if ( $self->cleanup_prod ) {
$cleanprod = $self->cleanup_product($prod);
if ( $cleanprod ne $prod ) {
$self->msg("Modify product: $prod => $cleanprod");
if ( $cleanprod eq $HYPO ) {
$cds{$pid}->add_tag_value( 'note', $prod );
$cds{$pid}->remove_tag('gene') if $cds{$pid}->has_tag('gene');
$cds{$pid}->remove_tag('EC_number') if $cds{$pid}->has_tag('EC_number');
}
$num_cleaned++;
}
}
$cds{$pid}->add_tag_value( 'product', $cleanprod );
$cds{$pid}->add_tag_value( 'EC_number', $EC ) if $EC;
if ( defined($gene) && $gene ne "" && !$cds{$pid}->has_tag('gene') ) {
$cds{$pid}->add_tag_value( 'gene', $gene );
}
$cds{$pid}->add_tag_value( 'inference', $db->{SRC} . $hit->name );
unlink "$tempdir/$pid.seq.out";
}
}
unlink map { "$tempdir/$_.seq" } keys %cds;
unlink map { "$tempdir/$_.seq.out" } keys %cds;
}
unlink map { "$tempdir/$_.seq" } keys %cds;
unlink map { "$tempdir/$_.seq.out" } keys %cds;
$self->msg("Cleaned $num_cleaned /product names") if $num_cleaned > 0;
}
}
if ($proteins) {
$self->delfile( map { "$outdir/proteins.$_" } qw(psq phr pin) );
}
# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
# Label unannotated proteins as 'hypothetical protein'
my $empty_label = $fast ? 'unannotated protein' : $HYPO;
my $num_hypo = 0;
for my $sid ( sort keys %seq ) {
for my $f ( @{ $seq{$sid}{FEATURE} } ) {
if ( $f->primary_tag eq 'CDS' and not $f->has_tag('product') ) {
$f->add_tag_value( 'product', $empty_label );
$num_hypo++;
}
}
}
$self->msg("Labelling remaining $num_hypo proteins as '$empty_label'") if $num_hypo > 0;
# . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
# Look for possible /pseudo genes - adjacent with same annotation
( run in 1.302 second using v1.01-cache-2.11-cpan-af0e5977854 )