-
Notifications
You must be signed in to change notification settings - Fork 1
Expand file tree
/
Copy pathMapRoutines.pm
More file actions
65 lines (56 loc) · 2.58 KB
/
Copy pathMapRoutines.pm
File metadata and controls
65 lines (56 loc) · 2.58 KB
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
package MapRoutines;
use strict;
use warnings;
use Bio::SeqIO;
use IO::String;
use File::Slurp;
use Exporter;
use base qw( Exporter );
our @EXPORT_OK = qw( get_gene translateToProt );
sub get_gene{
my ($genomeHash,$gbkHash,$ID,$ext_start,$ext_stop) = @_;
#The soubroutine will retrieve a subseq (e.g. a gene) from a genome
#Arguments: 1 = path to genome, 2 = %gbkhash, 3 = protein_id, 4 = extension of startposition backwards (optional), 5 = extension of stop position (optional)
my $specieID = $gbkHash->{$ID}[5];
if (not defined $specieID) {
die "Error: The protein id $ID was not found in any annotation file";
}
my $genomeHash_value = $genomeHash->{$specieID};
if (not defined $genomeHash_value) {
die "Error: The species id $specieID was not found in any genome file";
}
my $genome_string = ">$specieID\n" . $genomeHash_value;
my $io = IO::String->new($genome_string);
my $genome = Bio::SeqIO->new(-fh => $io, -format => 'fasta');
my $gene;
while (my $seq = $genome->next_seq){
#Start and stop positions from %gbkhash
my ($start, $stop);
#If gene is on complementary strand, take the reverse complement of the seq
if($gbkHash->{$ID}[2] == 0){
$start = $gbkHash->{$ID}[0] - $ext_start;
$stop = $gbkHash->{$ID}[1] + $ext_stop;
$gene = $seq->subseq($start, $stop);
}
else{
$start = $gbkHash->{$ID}[0] - $ext_stop;
$stop = $gbkHash->{$ID}[1] + $ext_start;
my $gene_obj = $seq->trunc($start, $stop);
my $reversed = $gene_obj->revcom;
$gene = $reversed->seq;
}
}
return $gene;
}
sub translateToProt {
my $sequence = $_[0];
my $DNA = Bio::PrimarySeq->new ( -seq => $sequence ,
#-id => 'YP_989117.1',
#-accession_number => 'X78121',
-alphabet => 'dna',
-is_circular => 0 );
my $protein = $DNA->translate;
#print substr $protein->seq,0,-1 . "\n";
return substr $protein->seq,0,-1;
}
1;