load_feature_mappings.pl File Reference
Go to the source code of this file.
Functions | |
| public | map_slice () |
| public | get_old_new_feature () |
Function Documentation
◆ get_old_new_feature()
| public get_old_new_feature | ( | ) |
Undocumented method
Code:
click to view
sub get_old_new_feature {
my $old_adaptor = shift;
my $new_adaptor = shift;
my $feature_type = shift;
my $chr = shift;
my $new_va = shift;
$support->log_stamped("fetching $feature_type mappings for chromosome $chr\n",1);
#get all genes on the chromosome from the old database
my $old_slice = $old_sa->fetch_by_region('chromosome', $chr, undef, undef,
undef, $support->param('altassembly'));
# get a slice on the old assembly from the new database
my $newdb_oldasm_chr = $new_sa->fetch_by_region('chromosome', $chr, undef, undef,
undef, $support->param('altassembly'));
my @old_features = @{$old_adaptor->fetch_all_by_Slice($old_slice)};
foreach my $old_feature (@old_features) {
my $feature_name;
if ($feature_type eq 'gene') {
$feature_name = $old_feature->stable_id;
} elsif ($feature_type eq 'SNP') {
$feature_name = $old_feature->variation_name;
}
$support->log_verbose("$feature_type $feature_name location in old db: " . $old_feature->start . ":" . $old_feature->end . ":" . $old_feature->strand . " chromosome $chr\n",1);
if ($old_feature->length == 0) {
next;
}
my $newdb_old_asm_slice = Bio::EnsEMBL::Slice->new(-coord_system => $newdb_oldasm_chr->coord_system,
-start => $old_feature->start,
-end => $old_feature->end,
-strand => $old_feature->strand,
-seq_region_name => $chr,
-seq_region_length => $old_feature->length,
-adaptor => $new_sa);
my $mapping_start = $newdb_oldasm_chr->end;
my ($mapping_end, $mapping_length, $mapping_strands, $mapping_chrs, $mapping);
($mapping_start, $mapping_end, $mapping_length, $mapping_strands,$mapping_chrs, $mapping) = map_slice($newdb_old_asm_slice, $mapping_start);
#get the feature from the new db
my $new_feature;
if ($feature_type eq 'gene') {
$new_feature = $new_adaptor->fetch_by_stable_id($feature_name);
} elsif ($feature_type eq 'SNP') {
my $variation = $new_va->fetch_by_name($feature_name);
if ($variation) {
my @new_features = @{$new_adaptor->fetch_all_by_Variation($variation)};
foreach my $vf (@new_features) {
if ($vf->allele_string eq $old_feature->allele_string) {
$new_feature = $vf;
}
}
}
}
my $new_feature_found = 0;
my ($new_chr,$new_start,$new_end,$new_length,$new_strand);
if (defined($new_feature) ) {
$new_chr = $new_feature->slice->seq_region_name;
$new_start = $new_feature->start;
$new_end = $new_feature->end;
$new_strand = $new_feature->strand;
$new_length = $new_feature->end - $new_feature->start + 1;
$new_feature_found = 1;
$support->log_verbose("$feature_type $feature_name found in the new db; location:$new_start:$new_end:$new_strand chromosome $new_chr\n",1);
} else {
$support->log_verbose("$feature_type $feature_name not found in the new db; can't compare $feature_type location to mapped location\n",1);
}
#store mapping in the db
if ( !$support->param('dry_run') ) {
$sth->execute($feature_type,$feature_name,$support->param('altassembly'),$chr,$old_feature->start, $old_feature->end, $old_feature->end - $old_feature->start + 1, $old_feature->strand, $new_feature_found, $new_chr, $new_start, $new_end,$new_length, $new_strand, $support->param('assembly'), $mapping, $mapping_start, $mapping_end, $mapping_length, $mapping_strands, $mapping_chrs );
}
}
}
◆ map_slice()
| public map_slice | ( | ) |
Undocumented method
Code:
click to view
sub map_slice;
sub get_old_new_feature;
Undocumented method
Code:
click to view
sub map_slice {
my $newdb_old_asm_slice = shift;
my $new_asm_start = shift;
$support->log_verbose("projection to new assembly: \n",1);
#project to new assembly
my @segments = @{ $newdb_old_asm_slice->project('chromosome', $support->param('assembly')) };
my $new_asm_end = 0;
my %mapping_strands;
my %mapping_chrs;
my $first_strand;
my $first_chr;
my $mapping_length = 0;
foreach my $segment (@segments) {
my $p_slice = $segment->to_Slice;
if (! defined($first_strand)) {
$first_strand = $p_slice->strand;
}
if (! defined($first_chr) ){
$first_chr = $p_slice->seq_region_name;
}
if ($p_slice->start < $new_asm_start && $p_slice->strand == $first_strand && $p_slice->seq_region_name eq $first_chr) {
$new_asm_start = $p_slice->start;
}
if ($p_slice->end > $new_asm_end && $p_slice->strand == $first_strand && $p_slice->seq_region_name eq $first_chr ) {
$new_asm_end = $p_slice->end;
}
if ($new_asm_end > 0 && $p_slice->strand == $first_strand && $p_slice->seq_region_name eq $first_chr ) {
$mapping_length += $p_slice->length;
}
$mapping_strands{$p_slice->strand} = 1;
$mapping_chrs{$p_slice->seq_region_name} = 1;
$support->log_verbose($p_slice->start() . ":" . $p_slice->end() . ":" . $p_slice->strand . " chromosome ".$p_slice->seq_region_name ."\n",1);
}
my $mapping; #0-no mapping, 1-gaps exist, 2-no gaps
if (scalar(@segments) == 0) {
$support->log_verbose("no mapping available\n",1);
$mapping = 0;
undef $new_asm_start;
undef $new_asm_end;
undef $mapping_length;
} else {
$support->log_verbose("maping start $new_asm_start mapping end $new_asm_end\n",1);
if (scalar(@segments) > 1) {
$support->log_verbose("gaps in new assembly sequence mapping\n",1);
$mapping = 1;
} else {
$support->log_verbose("no gaps in new assembly sequence mapping\n",1);
$mapping = 2;
}
}
my $mapping_strands;
if (scalar(keys %mapping_strands) > 1) {
$mapping_strands = 'both';
}
if (scalar(keys %mapping_strands) == 1) {
my @strands = keys %mapping_strands;
$mapping_strands = $strands[0];
}
my $mapping_chrs;
if (scalar(keys %mapping_chrs) >= 1) {
my @chrs = keys %mapping_chrs;
$mapping_chrs = join(',',@chrs);
}
return($new_asm_start, $new_asm_end, $mapping_length, $mapping_strands,$mapping_chrs, $mapping);
}
