0
|
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;
|