Mercurial > repos > mahtabm > ensembl
comparison variant_effect_predictor/Bio/EnsEMBL/DBSQL/CompressedSequenceAdaptor.pm @ 0:1f6dce3d34e0
Uploaded
author | mahtabm |
---|---|
date | Thu, 11 Apr 2013 02:01:53 -0400 |
parents | |
children |
comparison
equal
deleted
inserted
replaced
-1:000000000000 | 0:1f6dce3d34e0 |
---|---|
1 =head1 LICENSE | |
2 | |
3 Copyright (c) 1999-2012 The European Bioinformatics Institute and | |
4 Genome Research Limited. All rights reserved. | |
5 | |
6 This software is distributed under a modified Apache license. | |
7 For license details, please see | |
8 | |
9 http://www.ensembl.org/info/about/code_licence.html | |
10 | |
11 =head1 CONTACT | |
12 | |
13 Please email comments or questions to the public Ensembl | |
14 developers list at <dev@ensembl.org>. | |
15 | |
16 Questions may also be sent to the Ensembl help desk at | |
17 <helpdesk@ensembl.org>. | |
18 | |
19 =cut | |
20 | |
21 =head1 NAME | |
22 | |
23 Bio::EnsEMBL::DBSQL::CompressedSequenceAdaptor - Facilitates DB storage and retrieval of compressed sequence | |
24 | |
25 =head1 SYNOPSIS | |
26 | |
27 $seq_adptr = $database_adaptor->get_SequenceAdaptor(); | |
28 | |
29 $dna = | |
30 ${ $seq_adptr->fetch_by_Slice_start_end_strand( $slice, 1, 1000, | |
31 -1 ) }; | |
32 | |
33 =head1 DESCRIPTION | |
34 | |
35 An adaptor for the retrieval of compressed DNA sequence from the EnsEMBL | |
36 database | |
37 | |
38 =head1 METHODS | |
39 | |
40 =cut | |
41 | |
42 package Bio::EnsEMBL::DBSQL::CompressedSequenceAdaptor; | |
43 | |
44 use vars qw(@ISA); | |
45 use strict; | |
46 | |
47 use Bio::EnsEMBL::DBSQL::SequenceAdaptor; | |
48 | |
49 @ISA = qw(Bio::EnsEMBL::DBSQL::SequenceAdaptor); | |
50 | |
51 | |
52 sub _fetch_seq { | |
53 my $self = shift; | |
54 my $seq_region_id = shift; | |
55 my $start = shift; | |
56 my $len = shift; | |
57 | |
58 #calculate the offset and start in the compressed sequence | |
59 my $comp_start = ($start-1 >> 2) + 1; | |
60 my $comp_len = ($len >> 2) + 2; | |
61 | |
62 my ($bvector, $nline); | |
63 | |
64 my $sth = $self->prepare( | |
65 "SELECT SUBSTRING( d.sequence, ?, ?), n_line | |
66 FROM dnac d | |
67 WHERE d.seq_region_id = ?"); | |
68 $sth->bind_param(1,$comp_start,SQL_INTEGER); | |
69 $sth->bind_param(2,$comp_len ,SQL_INTEGER); | |
70 $sth->bind_param(3,$seq_region_id,SQL_INTEGER); | |
71 $sth->execute(); | |
72 $sth->bind_columns(\$bvector, \$nline); | |
73 $sth->fetch(); | |
74 $sth->finish(); | |
75 | |
76 #convert sequence from binary string to 0123 string | |
77 my $bitlen = length($bvector) << 2; | |
78 my $str = ''; | |
79 for(my $i=0; $i < $bitlen; $i++) { | |
80 $str .= vec($bvector, $i, 2); | |
81 } | |
82 | |
83 #convert from 0123 to ACTG | |
84 $str =~ tr/0123/ACTG/; | |
85 | |
86 $str = substr($str, ($start-1)%4, $len); | |
87 | |
88 #expand the nlines and place them back in the sequence | |
89 my @nlines = split(/:/, $nline); | |
90 foreach my $nl (@nlines) { | |
91 my ($offset,$char,$nlen) = $nl =~ /(\d+)(\D)(\d+)/; | |
92 | |
93 #skip nlines entirely out of range | |
94 next if(($offset+$nlen-1) < $start || $offset > ($start+$len-1)); | |
95 | |
96 #obtain relative offset into requested region | |
97 $offset = $offset - $start + 1; | |
98 | |
99 #nlines that partially overlap requested region have to be shrunk | |
100 if($offset < 1) { | |
101 $nlen = $nlen - (1-$offset); | |
102 $offset = 1; | |
103 } | |
104 if($offset + $nlen > $start+$len) { | |
105 $nlen = $len - $offset + 1; | |
106 } | |
107 | |
108 substr($str,$offset-1,$nlen) = $char x $nlen; | |
109 } | |
110 | |
111 return \$str; | |
112 } | |
113 | |
114 | |
115 =head2 store | |
116 | |
117 Arg [1] : string $seq_region_id the id of the sequence region this dna | |
118 will be associated with. | |
119 Arg [2] : string reference $sequence the dna sequence to be stored in | |
120 the database | |
121 Example : $dbID = $seq_adaptor->store(12,\'ACTGGGTACCAAACAAACACAACA'); | |
122 Description: stores a dna sequence in the databases dna table and returns the | |
123 database identifier for the new record. | |
124 Returntype : int | |
125 Exceptions : none | |
126 Caller : Bio::EnsEMBL::DBSQL::RawContigAdaptor::store | |
127 Status : Stable | |
128 | |
129 =cut | |
130 | |
131 sub store { | |
132 my ($self, $seq_region_id, $sequence) = @_; | |
133 | |
134 if(!$seq_region_id) { | |
135 throw('seq_region_id is required'); | |
136 } | |
137 | |
138 $sequence = uc($sequence); | |
139 | |
140 my $bvector = ''; | |
141 | |
142 #convert sequence to 0s,1s,2s and 3s | |
143 $sequence =~ tr/ACTG/0123/; | |
144 | |
145 #nlines cover sequence which is not ACTG such as N | |
146 #nline format is a set of colon delimited int, char, int triplets: | |
147 #<offset><code><length> | |
148 my($nline_char,$nline_len,$nline_off); | |
149 my @nlines; | |
150 | |
151 my $len = length($sequence); | |
152 for(my $i=0; $i < $len; $i++) { | |
153 my $char = substr($sequence,$i,1); | |
154 | |
155 #quickly check if this character was an A,C,T or G (and was converted to | |
156 # a 0,1,2,3) | |
157 if($char =~ /[0-3]/) { | |
158 vec($bvector, $i,2) = $char; | |
159 if($nline_char) { | |
160 #end of an nline | |
161 push @nlines, "$nline_off$nline_char$nline_len"; | |
162 $nline_char = undef; | |
163 $nline_len = 0; | |
164 $nline_off = 0; | |
165 } | |
166 } else { | |
167 #this was not an ACTG | |
168 if($nline_char) { | |
169 if($nline_char eq $char) { | |
170 #continuation of an nline | |
171 $nline_len++; | |
172 } else { | |
173 #end of a previous nline and start of a new one | |
174 push @nlines, "$nline_off$nline_char$nline_len"; | |
175 $nline_char = $char; | |
176 $nline_len = 1; | |
177 $nline_off = $i+1; | |
178 } | |
179 } else { | |
180 #start of a new nline | |
181 $nline_char = $char; | |
182 $nline_len = 1; | |
183 $nline_off = $i+1; | |
184 } | |
185 $char = 0; #need to put numeric val into bitvector despite nline | |
186 } | |
187 | |
188 vec($bvector, $i,2) = $char; | |
189 } | |
190 | |
191 my $nline = join(':', @nlines); | |
192 my $statement = $self->prepare( | |
193 "INSERT INTO dnac(seq_region_id, sequence, n_line) VALUES(?,?,?)"); | |
194 | |
195 $statement->bind_param(1,$seq_region_id,SQL_INTEGER); | |
196 $statement->bind_param(2,$bvector,SQL_BLOB); | |
197 $statement->bind_param(3,$nline,SQL_LONGVARCHAR); | |
198 $statement->execute(); | |
199 | |
200 $statement->finish(); | |
201 return; | |
202 } | |
203 | |
204 | |
205 1; |