XrefMapper::ProcessMappings Class Reference
Public Member Functions | |
| public | new () |
| public | process_mappings () |
| public | process_map_file () |
Detailed Description
Definition at line 7 of file ProcessMappings.pm.
Member Function Documentation
◆ new()
| public XrefMapper::ProcessMappings::new | ( | ) |
Undocumented method
Code:
click to view
sub new {
my($class, $mapper) = @_;
my $self ={};
bless $self,$class;
$self->core($mapper->core);
$self->xref($mapper->xref);
$self->verbose($mapper->verbose);
return $self;
}
◆ process_map_file()
| public XrefMapper::ProcessMappings::process_map_file | ( | ) |
Undocumented method
Code:
click to view
sub process_map_file{
my ($self, $map_file, $query_cutoff, $target_cutoff, $job_id, $array_number, $dbi) = @_;
my $ret = 1;
my $ensembl_type = "Translation";
if($map_file =~ /dna_/){
$ensembl_type = "Transcript";
}
my $mh;
if(!(open $mh ,"<",$map_file) ){
print STDERR "Could not open file $map_file\n Resubmit this job??\n";
return -1;
}
my $total_lines = 0;
my $root_dir = $self->core->dir;
my $dep_sth = $dbi->prepare("select dependent_xref_id, linkage_annotation from dependent_xref where master_xref_id = ?");
my $start_sth = $dbi->prepare("update mapping_jobs set object_xref_start = ? where job_id = ? and array_number = ?");
my $end_sth = $dbi->prepare("update mapping_jobs set object_xref_end = ? where job_id = ? and array_number = ?");
# my $update_dependent_xref_sth = $dbi->prepare("update dependent_xref set object_xref_id = ? where master_xref_id = ? and dependent_xref_id =?");
my $object_xref_id;
my $sth = $dbi->prepare("select max(object_xref_id) from object_xref");
$sth->execute();
$sth->bind_columns(\$object_xref_id);
$sth->fetch();
$sth->finish;
if(!defined($object_xref_id)){
$object_xref_id = 0;
}
my $object_xref_sth = $dbi->prepare("insert into object_xref (ensembl_id,ensembl_object_type, xref_id, linkage_type, ox_status ) values (?, ?, ?, ?, ?)");
my $object_xref_sth2 = $dbi->prepare("insert into object_xref (ensembl_id,ensembl_object_type, xref_id, linkage_type, ox_status, master_xref_id ) values (?, ?, ?, ?, ?, ?)");
my $get_object_xref_id_sth = $dbi->prepare("select object_xref_id from object_xref where ensembl_id = ? and ensembl_object_type = ? and xref_id = ? and linkage_type = ? and ox_status = ?");
my $get_object_xref_id_master_sth = $dbi->prepare("select object_xref_id from object_xref where ensembl_id = ? and ensembl_object_type = ? and xref_id = ? and linkage_type = ? and ox_status = ? and master_xref_id = ?");
local $object_xref_sth->{RaiseError}; #catch duplicates
local $object_xref_sth->{PrintError}; # cut down on error messages
local $object_xref_sth2->{RaiseError}; #catch duplicates
local $object_xref_sth2->{PrintError}; # cut down on error messages
my $identity_xref_sth = $dbi->prepare("insert ignore into identity_xref (object_xref_id, query_identity, target_identity, hit_start, hit_end, translation_start, translation_end, cigar_line, score ) values (?, ?, ?, ?, ?, ?, ?, ?, ?)");
my $ins_dep_ix_sth = $dbi->prepare("insert ignore into identity_xref (object_xref_id, query_identity, target_identity) values(?, ?, ?)");
my $source_name_sth = $dbi->prepare("select s.name from xref x join source s using(source_id) where x.xref_id = ?");
my $biotype_sth = $dbi->prepare("select biotype from transcript_stable_id where internal_id = ?");
my $last_query_id = 0;
my $best_match_found = 0;
my $best_identity = 0;
my $best_score = 0;
my %mRNA_biotypes = ( 'protein_coding' => 1,
'TR_C_gene' => 1,
'IG_V_gene' => 1,
'nonsense_mediated_decay' => 1,
'polymorphic_pseudogene' => 1);
my $first = 1;
while(<$mh>){
my $load_object_xref = 0;
$total_lines++;
chomp();
my ($label, $query_id, $target_id, $identity, $query_length, $target_length, $query_start, $query_end, $target_start, $target_end, $cigar_line, $score) = split(/:/, $_);
if ($last_query_id != $query_id) {
$best_match_found = 0;
$best_score = 0;
$best_identity = 0;
} else {
#ignore mappings with worse identity or score if we already found a good mapping
if ( ($identity < $best_identity || $score < $best_score) && $best_match_found) {
next;
}
}
if ($ensembl_type eq "Translation"){
$load_object_xref = 1;
} else {
#check if source name is RefSeq_ncRNA or RefSeq_mRNA
#if yes check biotype, if ok store object xref
$source_name_sth->execute($query_id);
my ($source_name) = $source_name_sth->fetchrow_array;
if ($source_name && ($source_name =~ /^RefSeq_(m|nc)RNA/ || $source_name =~ /^miRBase/ || $source_name =~ /^RFAM/)) {
#make sure mRNA xrefs are matched to protein_coding biotype only
$biotype_sth->execute($target_id);
my ($biotype) = $biotype_sth->fetchrow_array;
if ($source_name =~ /^RefSeq_mRNA/ && exists($mRNA_biotypes{$biotype})) {
$load_object_xref = 1;
}
if ($source_name =~ /^RefSeq_ncRNA/ && !exists($mRNA_biotypes{$biotype}) ) {
$load_object_xref = 1;
}
if (($source_name =~ /miRBase/ || $source_name =~ /^RFAM/) && $biotype =~ /RNA/) {
$load_object_xref = 1;
}
} else {
$load_object_xref = 1;
}
}
$last_query_id = $query_id;
if ($score > $best_score || $identity > $best_identity) {
$best_score = $score;
$best_identity = $identity;
}
if (!$load_object_xref) {
next;
} else {
$best_match_found = 1;
}
if(!defined($score)){
$end_sth->execute(($object_xref_id),$job_id, $array_number);
die "No score on line. Possible file corruption\n$_\n";
}
# calculate percentage identities
my $query_identity = int (100 * $identity / $query_length);
my $target_identity = int (100 * $identity / $target_length);
my $status = "DUMP_OUT";
# Only keep alignments where both sequences match cutoff
if($query_identity < $query_cutoff or $target_identity < $target_cutoff){
$status = "FAILED_CUTOFF";
}
$object_xref_sth->execute($target_id, $ensembl_type, $query_id, 'SEQUENCE_MATCH', $status) ;
$get_object_xref_id_sth->execute($target_id, $ensembl_type, $query_id, 'SEQUENCE_MATCH', $status);
$object_xref_id = ($get_object_xref_id_sth->fetchrow_array())[0];
if($object_xref_sth->err){
my $err = $object_xref_sth->errstr;
if($err =~ /Duplicate/){
# can get match from same xref and ensembl entity e.g.
# ensembl/ExonerateGappedBest1_dna_569.map:xref:934818:155760:54:1617:9648:73:12:3456:3517: M 61:242
# ensembl/ExonerateGappedBest1_dna_569.map:xref:934818:151735:58:1617:10624:73:6:5329:5397: M 48 D 1 M 19:242
next;
}
else{
$end_sth->execute(($object_xref_id),$job_id, $array_number);
die "Problem loading error is $err\n";
}
}
if($first){
$start_sth->execute($object_xref_id,$job_id, $array_number);
$first = 0;
}
$cigar_line =~ s/ //g;
$cigar_line =~ s/([MDI])(\d+)/$2$1/ig;
if(!$identity_xref_sth->execute($object_xref_id, $query_identity, $target_identity, $query_start+1, $query_end, $target_start+1, $target_end, $cigar_line, $score)){
$end_sth->execute(($object_xref_id),$job_id, $array_number);
die "Problem loading identity_xref";
}
my @master_xref_ids;
push @master_xref_ids, $query_id;
while(my $master_xref_id = pop(@master_xref_ids)){
my ($dep_xref_id, $link);
$dep_sth->execute($master_xref_id);
$dep_sth->bind_columns(\$dep_xref_id, \$link);
while($dep_sth->fetch){
$object_xref_sth2->execute($target_id, $ensembl_type, $dep_xref_id, 'DEPENDENT', $status, $master_xref_id);
$get_object_xref_id_master_sth->execute($target_id, $ensembl_type, $dep_xref_id, 'DEPENDENT', $status, $master_xref_id);
$object_xref_id = ($get_object_xref_id_master_sth->fetchrow_array())[0];
if($object_xref_sth2->err){
my $err = $object_xref_sth->errstr;
if($err =~ /Duplicate/){
# $duplicate_dependent_count++;
next;
}
else{
$end_sth->execute($object_xref_id,$job_id, $array_number);
die "Problem loading error is $err\n";
}
}
if($object_xref_sth2->err){
print STDERR "WARNING: Should not reach here??? object_xref_id = $object_xref_id\n";
}
$ins_dep_ix_sth->execute($object_xref_id, $query_identity, $target_identity);
# now store in object_xref $update_dependent_xref_sth->execute($object_xref_id, $master_xref_id, $dep_xref_id);
push @master_xref_ids, $dep_xref_id; # get the dependent, dependents just in case
}
}
}
close $mh;
$end_sth->execute($object_xref_id,$job_id, $array_number);
$start_sth->finish;
$end_sth->finish;
$dep_sth->finish;
$object_xref_sth->finish;
$identity_xref_sth->finish;
$ins_dep_ix_sth->finish;
$source_name_sth->finish;
$biotype_sth->finish;
return $total_lines;
}
◆ process_mappings()
| public XrefMapper::ProcessMappings::process_mappings | ( | ) |
Undocumented method
Code:
click to view
sub process_mappings {
my ($self) = @_;
# get the jobs from the mapping table
# for i =1 i < jobnum{
# if( Not parsed already see mapping_jobs){
# check the .err, .out and .map files in that order.
# put data into object_xref, identity_xref etc...
# add data to mapping_job table
# }
# }
my %query_cutoff;
my %target_cutoff;
my ($job_id, $percent_query_cutoff, $percent_target_cutoff);
my $dbi = $self->xref->dbc();
my $sth = $dbi->prepare("select job_id, percent_query_cutoff, percent_target_cutoff from mapping");
$sth->execute();
$sth->bind_columns(\$job_id, \$percent_query_cutoff, \$percent_target_cutoff);
my $object_xref_id;
while($sth->fetch){
$query_cutoff{$job_id} = $percent_query_cutoff;
$target_cutoff{$job_id} = $percent_target_cutoff;
}
$sth->finish;
my ($root_dir, $map, $status, $out, $err, $array_number);
my ($map_file, $out_file, $err_file);
my $map_sth = $dbi->prepare("select root_dir, map_file, status, out_file, err_file, array_number, job_id from mapping_jobs");
$map_sth->execute();
$map_sth->bind_columns(\$root_dir, \$map, \$status, \$out, \$err, \$array_number, \$job_id);
my $already_processed_count = 0;
my $processed_count = 0;
my $error_count = 0;
my $empty_count = 0;
my $stat_sth = $dbi->prepare("update mapping_jobs set status = ? where job_id = ? and array_number = ?");
while($map_sth->fetch()){
my $err_file = $root_dir."/".$err;
my $out_file = $root_dir."/".$out;
my $map_file = $root_dir."/".$map;
if($status eq "SUCCESS"){
$already_processed_count++;
}
else{
if(-s $err_file){
$error_count++;
print STDERR "Problem $err_file is non zero\n";
if(open my $eh ,"<", $err_file){
while(<$eh>){
print STDERR "#".$_;
}
close $eh;
}
else{
print STDERR "No file exists $err_file???\n Resubmit this job\n";
}
if($status eq "SUBMITTED"){
$stat_sth->execute('FAILED',$job_id, $array_number);
}
}
else{ #err file checks out so process the mapping file.
if(-e $map_file){
my $count = $self->process_map_file($map_file, $query_cutoff{$job_id}, $target_cutoff{$job_id}, $job_id, $array_number, $dbi);
if( $count > 0){
$processed_count++;
$stat_sth->execute('SUCCESS',$job_id, $array_number);
}
elsif($count ==0){
print STDERR "WARNING $map_file was empty could be okay but if there are alot of these it is probably a problem\n";
$processed_count++;
$empty_count++;
$stat_sth->execute('SUCCESS',$job_id, $array_number);
}
else{
$error_count++;
$stat_sth->execute('FAILED',$job_id, $array_number);
}
}
else{
$error_count++;
print STDERR "Could not open file $map_file???\n Resubmit this job\n";
$stat_sth->execute('FAILED',$job_id, $array_number);
}
}
}
}
$map_sth->finish;
$stat_sth->finish;
print "already processed = $already_processed_count, processed = $processed_count, errors = $error_count, empty = $empty_count\n" if($self->verbose);
if(!$error_count){
my $sth = $dbi->prepare("insert into process_status (status, date) values('mapping_processed',now())");
$sth->execute();
$sth->finish;
}
}
The documentation for this class was generated from the following file:
- ensembl/misc-scripts/xref_mapping/XrefMapper/ProcessMappings.pm
