$Bio::EnsEMBL::DBSQL::DensityFeatureAdaptor::VERSION
=
'112.0_56'
;
$Bio::EnsEMBL::DBSQL::DensityFeatureAdaptor::VERSION
=
'112.056'
;
@ISA
=
qw(Bio::EnsEMBL::DBSQL::BaseFeatureAdaptor)
;
our
$DENSITY_FEATURE_CACHE_SIZE
= 20;
sub
new {
my
$caller
=
shift
;
my
$class
=
ref
(
$caller
) ||
$caller
;
my
$self
=
$class
->SUPER::new(
@_
);
my
%cache
;
tie
(
%cache
,
'Bio::EnsEMBL::Utils::Cache'
,
$DENSITY_FEATURE_CACHE_SIZE
);
$self
->{
'_density_feature_cache'
} = \
%cache
;
return
$self
;
}
sub
fetch_all_by_Slice {
my
(
$self
,
$slice
,
$logic_name
,
$num_blocks
,
$interpolate
,
$max_ratio
) =
@_
;
if
(
defined
(
$num_blocks
) &&
$num_blocks
< 1) {
warning(
"Invalid number of density blocks [$num_blocks] requested.\n"
.
"Returning empty list."
);
return
[];
}
$num_blocks
||= 50;
my
$length
=
$slice
->
length
();
my
$wanted_block_size
= POSIX::ceil(
$length
/
$num_blocks
);
my
$dta
=
$self
->db()->get_DensityTypeAdaptor();
my
@dtypes
= @{
$dta
->fetch_all_by_logic_name(
$logic_name
)};
if
( !
@dtypes
){
my
@all_dtypes
= @{
$dta
->fetch_all() };
@all_dtypes
or warning(
"No DensityTypes in $dta"
) &&
return
[];
my
$valid_list
=
join
(
", "
,
map
{
$_
->analysis->logic_name}
@all_dtypes
);
warning(
"No DensityTypes for logic name $logic_name. "
.
"Select from $valid_list"
);
return
[];
}
my
$best_ratio
=
undef
;
my
$density_type
=
undef
;
my
$best_ratio_large
=
undef
;
my
$density_type_large
=
undef
;
my
%dt_ratio_hash
;
foreach
my
$dt
(
@dtypes
) {
my
$ratio
;
if
(
$dt
->block_size() > 0 ) {
$ratio
=
$wanted_block_size
/
$dt
->block_size();
}
else
{
my
$block_size
=
$slice
->seq_region_length() /
$dt
->region_features();
$ratio
=
$wanted_block_size
/
$block_size
;
}
if
(
$ratio
< 1 ) {
$ratio
= 5/
$ratio
;
}
$dt_ratio_hash
{
$ratio
} =
$dt
;
}
$best_ratio
= (
sort
{
$a
<=>
$b
} (
keys
%dt_ratio_hash
))[0];
if
(!
defined
(
$best_ratio
) ||
(
defined
(
$max_ratio
) &&
$best_ratio
>
$max_ratio
)) {
return
[];
}
$density_type
=
$dt_ratio_hash
{
$best_ratio
};
my
$constraint
=
"df.density_type_id = "
.
$density_type
->dbID();
my
@features
=
@{
$self
->fetch_all_by_Slice_constraint(
$slice
,
$constraint
)};
return
\
@features
if
(!
$interpolate
);
my
@out
;
@features
=
sort
( {
$a
->start() <=>
$b
->start() }
@features
);
my
$start
= 1;
my
$end
=
$start
+
$wanted_block_size
-1;
my
$value_type
=
$density_type
->value_type();
$density_type
= Bio::EnsEMBL::DensityType->new
(
-VALUE_TYPE
=>
$value_type
,
-ANALYSIS
=>
$density_type
->analysis(),
-BLOCK_SIZE
=>
$wanted_block_size
);
while
(
$start
<
$length
) {
my
$density_value
= 0.0;
my
(
$f
,
$fstart
,
$fend
,
$portion
);
my
@dvalues
;
while
((
$f
=
shift
(
@features
)) &&
$end
>
$f
->{
'start'
}) {
$fstart
= (
$f
->{
'start'
} <
$start
) ?
$start
:
$f
->{
'start'
};
$fend
= (
$f
->{
'end'
} >
$end
) ?
$end
:
$f
->{
'end'
};
$fend
=
$length
if
(
$fend
>
$length
);
if
(
$value_type
eq
'sum'
) {
$portion
= (
$fend
-
$fstart
+1)/
$f
->
length
();
$density_value
+=
$portion
*
$f
->{
'density_value'
};
}
elsif
(
$value_type
eq
'ratio'
) {
push
(
@dvalues
,[
$fend
-
$fstart
+1,
$f
->{
'density_value'
}]);
}
else
{
throw(
"Unknown density value type [$value_type]."
);
}
last
if
(
$fend
<
$f
->{
'end'
});
}
if
(
defined
(
$f
) && (!
defined
(
$fend
) ||
$fend
<
$f
->{
'end'
})) {
unshift
(
@features
,
$f
);
}
if
(
$value_type
eq
'ratio'
) {
my
$total_len
=
$end
-
$start
+ 1;
if
(
$total_len
> 0) {
foreach
my
$pair
(
@dvalues
) {
my
(
$dlen
,
$dval
) =
@$pair
;
$density_value
+=
$dval
* (
$dlen
/
$total_len
);
}
}
}
push
@out
, Bio::EnsEMBL::DensityFeature->new
(
-seq_region
=>
$slice
,
-start
=>
$start
,
-end
=>
$end
,
-density_type
=>
$density_type
,
-density_value
=>
$density_value
);
$start
=
$end
+ 1;
$end
+=
$wanted_block_size
;
}
return
\
@out
;
}
sub
fetch_all {
my
$self
=
shift
;
my
$logic_name
=
shift
;
my
(
$sth
,
$seq_region_id
,
$start
,
$end
,
$density_value
,
$density_type_id
);
my
(
$slice
,
$density_type
,
@out
);
my
$sa
=
$self
->db()->get_SliceAdaptor();
my
$dta
=
$self
->db()->get_DensityTypeAdaptor();
if
(
$logic_name
) {
$sth
=
$self
->prepare(
"SELECT df.seq_region_id, df.seq_region_start, df.seq_region_end, df.density_value, df.density_type_id "
.
"FROM density_feature df, density_type dt, analysis a "
.
"WHERE df.density_type_id = dt.density_type_id "
.
"AND dt.analysis_id = a.analysis_id "
.
"AND logic_name = ?"
);
$sth
->bind_param(1,
$logic_name
, SQL_VARCHAR);
$sth
->execute();
}
else
{
$sth
=
$self
->prepare(
"SELECT df.seq_region_id, df.seq_region_start, df.seq_region_end, df.density_value, df.density_type_id "
.
"FROM density_feature df, density_type dt "
.
"WHERE df.density_type_id = dt.density_type_id"
);
$sth
->execute();
}
$sth
->bind_columns(\(
$seq_region_id
,
$start
,
$end
,
$density_value
,
$density_type_id
));
while
(
$sth
->fetch()) {
$slice
=
$sa
->fetch_by_seq_region_id(
$seq_region_id
);
$density_type
=
$dta
->fetch_by_dbID(
$density_type_id
);
push
(
@out
, Bio::EnsEMBL::DensityFeature->new
(
-seq_region
=>
$slice
,
-start
=>
$start
,
-end
=>
$end
,
-density_type
=>
$density_type
,
-density_value
=>
$density_value
));
}
return
\
@out
;
}
sub
fetch_all_by_Slice_constraint {
my
$self
=
shift
;
my
$slice
=
shift
;
my
$constraint
=
shift
;
my
@result
;
my
$bind_params
=
$self
->bind_param_generic_fetch();
if
( !
ref
(
$slice
)
|| !(
$slice
->isa(
'Bio::EnsEMBL::Slice'
)
or
$slice
->isa(
'Bio::EnsEMBL::LRGSlice'
) ) )
{
throw(
"Bio::EnsEMBL::Slice argument expected."
);
}
$constraint
||=
''
;
if
( !
defined
(
$constraint
) ) {
return
[] }
my
$key
;
my
$cache
;
if
(
!(
defined
(
$self
->db()->no_cache() ) &&
$self
->db()->no_cache() ) )
{
if
(
$slice
->isa(
'Bio::EnsEMBL::Variation::StrainSlice'
) ) {
my
$string
=
$self
->dbc()->db_handle()->quote(
$slice
->strain_name() );
if
(
$constraint
ne
""
) {
$constraint
.=
" AND $string = $string "
;
}
else
{
$constraint
.=
" $string = $string "
;
}
}
$key
=
uc
(
join
(
':'
,
$slice
->name(),
$constraint
) );
if
(
defined
(
$bind_params
) ) {
$key
.=
':'
.
join
(
':'
,
map
{
$_
->[0] .
'/'
.
$_
->[1] } @{
$bind_params
} );
}
$cache
=
$self
->_slice_feature_cache();
if
(
exists
(
$cache
->{
$key
} ) ) {
$self
->{
'_bind_param_generic_fetch'
} = ();
return
$cache
->{
$key
};
}
}
$self
->_bind_param_generic_fetch(
$bind_params
);
my
$features
=
$self
->_slice_fetch(
$slice
,
$constraint
);
if
(
defined
(
$key
) ) {
$cache
->{
$key
} =
$features
;
}
return
$features
;
}
sub
_tables {
my
$self
=
shift
;
return
([
'density_feature'
,
'df'
]);
}
sub
_columns {
my
$self
=
shift
;
return
qw( df.density_feature_id
df.seq_region_id df.seq_region_start df.seq_region_end
df.density_value df.density_type_id )
;
}
sub
_objs_from_sth {
my
(
$self
,
$sth
,
$mapper
,
$dest_slice
) =
@_
;
my
$sa
=
$self
->db()->get_SliceAdaptor();
my
$dta
=
$self
->db()->get_DensityTypeAdaptor();
my
@features
;
my
%dtype_hash
;
my
%slice_hash
;
my
%sr_name_hash
;
my
%sr_cs_hash
;
my
(
$density_feature_id
,
$seq_region_id
,
$seq_region_start
,
$seq_region_end
,
$density_value
,
$density_type_id
);
$sth
->bind_columns(\(
$density_feature_id
,
$seq_region_id
,
$seq_region_start
,
$seq_region_end
,
$density_value
,
$density_type_id
));
my
$dest_slice_start
;
my
$dest_slice_end
;
my
$dest_slice_strand
;
my
$dest_slice_length
;
my
$dest_slice_cs
;
my
$dest_slice_sr_name
;
my
$dest_slice_sr_id
;
my
$asma
;
if
(
$dest_slice
) {
$dest_slice_start
=
$dest_slice
->start();
$dest_slice_end
=
$dest_slice
->end();
$dest_slice_strand
=
$dest_slice
->strand();
$dest_slice_length
=
$dest_slice
->
length
();
$dest_slice_cs
=
$dest_slice
->coord_system();
$dest_slice_sr_name
=
$dest_slice
->seq_region_name();
$dest_slice_sr_id
=
$dest_slice
->get_seq_region_id();
$asma
=
$self
->db->get_AssemblyMapperAdaptor();
}
FEATURE:
while
(
$sth
->fetch()) {
my
$density_type
=
$dtype_hash
{
$density_type_id
} ||=
$dta
->fetch_by_dbID(
$density_type_id
);
$seq_region_id
=
$self
->get_seq_region_id_internal(
$seq_region_id
);
my
$slice
=
$slice_hash
{
"ID:"
.
$seq_region_id
};
if
(!
$slice
) {
$slice
=
$sa
->fetch_by_seq_region_id(
$seq_region_id
);
$slice_hash
{
"ID:"
.
$seq_region_id
} =
$slice
;
$sr_name_hash
{
$seq_region_id
} =
$slice
->seq_region_name();
$sr_cs_hash
{
$seq_region_id
} =
$slice
->coord_system();
}
if
(!
$mapper
&&
$dest_slice
&& !
$dest_slice_cs
->equals(
$slice
->coord_system)) {
$mapper
=
$asma
->fetch_by_CoordSystems(
$dest_slice_cs
,
$slice
->coord_system);
}
my
$sr_name
=
$sr_name_hash
{
$seq_region_id
};
my
$sr_cs
=
$sr_cs_hash
{
$seq_region_id
};
if
(
$mapper
) {
my
$len
=
$seq_region_end
-
$seq_region_start
+ 1;
if
(
defined
$dest_slice
&&
$mapper
->isa(
'Bio::EnsEMBL::ChainedAssemblyMapper'
) ) {
(
$seq_region_id
,
$seq_region_start
,
$seq_region_end
) =
$mapper
->
map
(
$sr_name
,
$seq_region_start
,
$seq_region_end
, 1,
$sr_cs
, 0,
$dest_slice
);
}
else
{
(
$seq_region_id
,
$seq_region_start
,
$seq_region_end
) =
$mapper
->
map
(
$sr_name
,
$seq_region_start
,
$seq_region_end
, 1,
$sr_cs
);
}
next
FEATURE
if
(!
defined
(
$seq_region_id
));
if
(
$density_type
->value_type() eq
'sum'
) {
my
$newlen
=
$seq_region_end
-
$seq_region_start
+1;
$density_value
*=
$newlen
/
$len
if
(
$newlen
!=
$len
);
}
$slice
=
$slice_hash
{
"ID:"
.
$seq_region_id
} ||=
$sa
->fetch_by_seq_region_id(
$seq_region_id
);
}
if
(
defined
(
$dest_slice
)) {
my
$seq_region_len
=
$dest_slice
->seq_region_length();
if
(
$dest_slice_strand
== 1 ) {
$seq_region_start
=
$seq_region_start
-
$dest_slice_start
+ 1;
$seq_region_end
=
$seq_region_end
-
$dest_slice_start
+ 1;
if
(
$dest_slice
->is_circular ) {
if
(
$seq_region_start
>
$seq_region_end
) {
if
(
$seq_region_end
>
$dest_slice_start
) {
$seq_region_start
-=
$seq_region_len
;
}
if
(
$seq_region_end
< 0 ) {
$seq_region_end
+=
$seq_region_len
;
}
}
else
{
if
(
$dest_slice_start
>
$dest_slice_end
&&
$seq_region_end
< 0) {
$seq_region_start
+=
$seq_region_len
;
$seq_region_end
+=
$seq_region_len
;
}
}
}
}
else
{
my
$start
=
$dest_slice_end
-
$seq_region_end
+ 1;
my
$end
=
$dest_slice_end
-
$seq_region_start
+ 1;
if
(
$dest_slice
->is_circular()) {
if
(
$dest_slice_start
>
$dest_slice_end
) {
if
(
$seq_region_start
>=
$dest_slice_start
) {
$end
+=
$seq_region_len
;
$start
+=
$seq_region_len
if
$seq_region_end
>
$dest_slice_start
;
}
elsif
(
$seq_region_start
<=
$dest_slice_end
) {
}
elsif
(
$seq_region_end
>=
$dest_slice_start
) {
$start
+=
$seq_region_len
;
$end
+=
$seq_region_len
;
}
elsif
(
$seq_region_end
<=
$dest_slice_end
) {
$end
+=
$seq_region_len
if
$end
< 0;
}
elsif
(
$seq_region_start
>
$seq_region_end
) {
$end
+=
$seq_region_len
;
}
}
else
{
if
(
$seq_region_start
<=
$dest_slice_end
and
$seq_region_end
>=
$dest_slice_start
) {
}
elsif
(
$seq_region_start
>
$seq_region_end
) {
if
(
$seq_region_start
<=
$dest_slice_end
) {
$start
-=
$seq_region_len
;
}
elsif
(
$seq_region_end
>=
$dest_slice_start
) {
$end
+=
$seq_region_len
;
}
}
}
}
$seq_region_start
=
$start
;
$seq_region_end
=
$end
;
}
if
(
$seq_region_end
< 1
||
$seq_region_start
>
$dest_slice_length
|| (
$dest_slice_sr_id
!=
$seq_region_id
)) {
next
FEATURE;
}
$slice
=
$dest_slice
;
}
push
(
@features
,
$self
->_create_feature(
'Bio::EnsEMBL::DensityFeature'
, {
-dbID
=>
$density_feature_id
,
-adaptor
=>
$self
,
-start
=>
$seq_region_start
,
-end
=>
$seq_region_end
,
-seq_region
=>
$slice
,
-density_value
=>
$density_value
,
-density_type
=>
$density_type
} ) );
}
return
\
@features
;
}
sub
list_dbIDs {
my
(
$self
,
$ordered
) =
@_
;
return
$self
->_list_dbIDs(
"density_feature"
,
undef
,
$ordered
);
}
sub
store{
my
(
$self
,
@df
) =
@_
;
if
(
scalar
(
@df
) == 0 ) {
throw(
"Must call store with list of DensityFeatures"
);
}
my
$sth
=
$self
->prepare
(
"INSERT INTO density_feature (seq_region_id, seq_region_start, "
.
"seq_region_end, density_type_id, "
.
"density_value) "
.
"VALUES (?,?,?,?,?)"
);
my
$db
=
$self
->db();
my
$analysis_adaptor
=
$db
->get_AnalysisAdaptor();
FEATURE:
foreach
my
$df
(
@df
) {
if
( !
ref
$df
|| !
$df
->isa(
"Bio::EnsEMBL::DensityFeature"
) ) {
throw(
"DensityFeature must be an Ensembl DensityFeature, "
.
"not a ["
.
ref
(
$df
).
"]"
);
}
next
if
(
$df
->density_value == 0 );
if
(
$df
->is_stored(
$db
)) {
warning(
"DensityFeature ["
.
$df
->dbID.
"] is already stored"
.
" in this database."
);
next
FEATURE;
}
if
(!
defined
(
$df
->density_type)) {
throw(
"A density type must be attached to the features to be stored."
);
}
if
(!
$df
->density_type->is_stored(
$db
)) {
my
$dta
=
$db
->get_DensityTypeAdaptor();
$dta
->store(
$df
->density_type());
}
my
$original
=
$df
;
my
$seq_region_id
;
(
$df
,
$seq_region_id
) =
$self
->_pre_store(
$df
);
$sth
->bind_param(1,
$seq_region_id
,SQL_INTEGER);
$sth
->bind_param(2,
$df
->start,SQL_INTEGER);
$sth
->bind_param(3,
$df
->end,SQL_INTEGER);
$sth
->bind_param(4,
$df
->density_type->dbID,SQL_INTEGER);
$sth
->bind_param(5,
$df
->density_value,SQL_FLOAT);
$sth
->execute();
$original
->dbID(
$self
->last_insert_id(
'density_feature_id'
,
undef
,
'density_feature'
));
$original
->adaptor(
$self
);
}
}
sub
fetch_Featureset_by_Slice {
my
(
$self
,
$slice
,
$logic_name
,
$num_blocks
,
$interpolate
,
$max_ratio
) =
@_
;
my
$key
=
join
(
":"
,
$slice
->name,
$logic_name
,
$num_blocks
|| 50,
$interpolate
|| 0,
$max_ratio
);
unless
(
$self
->{
'_density_feature_cache'
}->{
$key
}) {
my
$dfeats
=
$self
->fetch_all_by_Slice(
$slice
,
$logic_name
,
$num_blocks
,
$interpolate
,
$max_ratio
);
$self
->{
'_density_feature_cache'
}->{
$key
} =
new Bio::EnsEMBL::DensityFeatureSet(
$dfeats
);
}
return
$self
->{
'_density_feature_cache'
}->{
$key
};
}
1;