Skip to content

Commit

Permalink
In progress
Browse files Browse the repository at this point in the history
  • Loading branch information
bosborne committed Jan 29, 2013
1 parent 44627ba commit b28cef3
Show file tree
Hide file tree
Showing 2 changed files with 37 additions and 25 deletions.
14 changes: 7 additions & 7 deletions lib/diya/BDRD/rpsblastCDS.pm
Expand Up @@ -109,15 +109,15 @@ sub parse {

if ( $feat->primary_tag eq 'CDS' ) {

my @locus = $feat->get_tag_values('locus_tag');
my ($locus) = $feat->get_tag_values('locus_tag');

my $fh;

eval { $fh = $index->get_stream( $locus[0] ); };
eval { $fh = $index->get_stream( $locus ); };

if ($@) {
print "Problem retrieving "
. $locus[0]
. $locus
. " from $blastout.index\n";
next;
}
Expand All @@ -141,7 +141,7 @@ sub parse {

print "Found rpsblast hit "
. $hit->description
. " for $locus[0]\n"
. " for $locus\n"
if $diya->verbose;

# get domain id from hit
Expand All @@ -166,11 +166,11 @@ sub parse {

}

my @products = $feat->get_tag_values('product')
my ($product) = $feat->get_tag_values('product')
if ( $feat->has_tag('product') );

# If this feature is already annotated with a UniRef match
if ( $products[0] =~ /RepID=/ && defined $tags->{'product'} ) {
if ( $product =~ /RepID=/ && defined $tags->{'product'} ) {

my $note = "similar to " . $tags->{'product'};
$note .= " " . $tags->{'cluster'};
Expand Down Expand Up @@ -332,7 +332,7 @@ sub load_cdd_table {
$cddmap{$id} = \%cluster;
}

# 103623 PRK08948 PRK08948 2-octaprenyl-6-methoxyphenyl hydroxylase; Validated 392
# 103623 PRK08948 PRK08948 2-octaprenyl-6-methoxyphenyl hydroxylase; Validated 392
if ( $line =~ /^\d+\s+(PRK\d+)\s+\w+\s+([^;]+)/ ) {
my ( $id, $def ) = ( $1, $2 );
$def =~ s/\s+$//;
Expand Down
48 changes: 30 additions & 18 deletions scripts/extractCDS.pl
Expand Up @@ -7,7 +7,7 @@ =head1 NAME
=head1 DESCRIPTION
Extracts the CDS features from a GenBank file and creates nucleotide and
peptide fasta files containing all genes and translations.
peptide fasta files containing them.
gene complement(32855..33520)
/locus_tag="test0001_380"
Expand All @@ -23,32 +23,44 @@ =head1 DESCRIPTION
use Bio::Seq;
use Data::Dumper;

my ($inputfile,$outputfile) = @ARGV;
my ( $inputfile, $outputfile ) = @ARGV;

die "Need both input file name and output file name"
die "Need both input file name and output file name"
unless ( $inputfile && $outputfile );

my $in = Bio::SeqIO->new( -file => "$inputfile.gbk",
-format => 'genbank' );
my $in = Bio::SeqIO->new(
-file => "$inputfile.gbk",
-format => 'genbank'
);

my $ntout = Bio::SeqIO->new( -file => ">>$outputfile",
-format => 'fasta' );
my $out = Bio::SeqIO->new(
-file => ">>$outputfile",
-format => 'fasta'
);

my $pepout = Bio::SeqIO->new( -file => ">>$outputfile.pep",
-format => 'fasta' );
my $pepout = Bio::SeqIO->new(
-file => ">>$outputfile.pep",
-format => 'fasta'
);

my $gbk = $in->next_seq;

for my $feature ( $gbk->get_SeqFeatures ) {
if ( $feature->primary_tag eq "gene" ) {
my @ids = $feature->get_tag_values("locus_tag");
my $str = $feature->seq->seq;
my $seq = Bio::Seq->new( -seq => $str,
-display_id => "@ids" );
$ntout->write_seq($seq);

$seq = $seq->revcom if ( $feature->location->start > $feature->location->end );
my $pep = $seq->translate(-complete => 1);
$pepout->write_seq($pep);
my @ids = $feature->get_tag_values("locus_tag");
my $str = $feature->seq->seq;
my $seq = Bio::Seq->new(
-seq => $str,
-display_id => "@ids"
);
$out->write_seq($seq);

$seq = $seq->revcom
if ( $feature->location->start > $feature->location->end );
my $pep = $seq->translate( -complete => 1 );
$pepout->write_seq($pep);
}
}

__END__

0 comments on commit b28cef3

Please sign in to comment.