Functions
public	get_options ()

public	run ()

public	run_on_dba ()

public	get_liftover_mappings ()

public	write_mappings ()

public	build_chain_mappings ()

public	get_assembly_mappings ()

public	get_reverse_assembly_mappings ()

public	get_sql ()

public	print_chains ()

public	create_ucsc_chains ()

public	ensembl_to_ucsc_name ()

Function Documentation

◆ build_chain_mappings()

public build_chain_mappings ( )

Undocumented method

Code:

click to view

sub build_chain_mappings {
  my ($assembly_mappings) = @_;
  my @chain_mappings;
  
  my ($t_name, $t_size, $t_strand, $t_start, $t_end);
  my ($q_name, $q_size, $q_strand, $q_start, $q_end);
  my $chain_id = 1;
  my @chain_gaps;
  
  my $length = scalar(@{$assembly_mappings});
  for (my $i = 0; $i < $length; $i++) {
 
    my $current = $assembly_mappings->[$i];
    my $next = ($i+1 != $length) ? $assembly_mappings->[$i+1] : undef;
 
    my $ori = $current->{ori};
    my ($asm_diff, $cmp_diff);
    if($next) {
      $asm_diff = ($next->{asm_start} - $current->{asm_end})-1;
      # Rev strands means the next cmp region has a lower start than the 
      # current end (because it's running in reverse). Rember length in 1-based
      # coords always is (high-low)-1
      $cmp_diff = ($ori == 1) ? ($next->{cmp_start} - $current->{cmp_end})-1 : ($current->{cmp_start} - $next->{cmp_end})-1;
    }
 
    if(! $t_name) {
      # Reset variables to current
      @chain_gaps = ();
      $chain_id++;
      ($t_name, $t_size, $t_strand, $t_start) = ($current->{asm_name}, $current->{asm_length}, 1, $current->{asm_start});
      ($q_name, $q_size, $q_strand) = ($current->{cmp_name}, $current->{cmp_length}, $current->{ori});
      $q_start = ($ori == 1) ? $current->{cmp_start} : $current->{cmp_end};
    }
 
    # Block that decides we need to start a new chain definition
    #
    # Can mean we have run out of mappings, we are into a new chromsome (both source and target) 
    # or strand has swapped.
    # Final reason is we've had an out-of-order meaning a negative gap was produced 
    # (we're going backwards). In any situation this means the chain is finished
    if( ! defined $next || 
        $t_name ne $next->{asm_name} ||
        $q_name ne $next->{cmp_name} || # we can switch target chromosomes. e.g. cross chromsome mappings in mouse NCBI37->GRCm38
        $ori != $next->{ori} ||
        $asm_diff < 0 ||
        $cmp_diff < 0) {
      # Add the last gap on which is just the length of this alignment
      push(@chain_gaps, [$current->{length}]);
      # Set the ends of the chain since this is the last block
      $t_end = $current->{asm_end};
      $q_end = ($ori == 1) ? $current->{cmp_end} : $current->{cmp_start};
 
      #If strand was negative we need to represent all data as reverse complemented regions
      if($q_strand == -1) {
        # $t_start = ($t_size - $t_start)+1;
        # $t_end = ($t_size - $t_end)+1;
        $q_start = ($q_size - $q_start)+1;
        $q_end = ($q_size - $q_end)+1;
      }
      # Convert to UCSC formats (0-based half-open intervals and +/- strands)
      $t_start--;
      $q_start--;
      $t_strand = ($t_strand == 1) ? '+' : '-';
      $q_strand = ($q_strand == 1) ? '+' : '-';
 
      #Store the chain
      my $chain_score = 1;
      push(@chain_mappings, {
        header => ['chain', $chain_score, $t_name, $t_size, $t_strand, $t_start, $t_end, $q_name, $q_size, $q_strand, $q_start, $q_end, $chain_id],
        gaps => [@chain_gaps]
      });
 
      if(! defined $next) {
        last;
      }
 
      # Clear variables
      ($t_name, $t_size, $t_strand, $t_start) = ();
      ($q_name, $q_size, $q_strand, $q_start) = ();
    }
 
    push(@chain_gaps, [$current->{length}, $asm_diff, $cmp_diff]);
  }
 
  return \@chain_mappings;
}

◆ create_ucsc_chains()

public create_ucsc_chains ( )

Undocumented method

Code:

click to view

◆ ensembl_to_ucsc_name()

public ensembl_to_ucsc_name ( )

Undocumented method

Code:

click to view

◆ get_assembly_mappings()

public get_assembly_mappings ( )

Undocumented method

Code:

click to view

◆ get_liftover_mappings()

public get_liftover_mappings ( )

Undocumented method

Code:

click to view

◆ get_options()

public get_options ( )

Undocumented method

Code:

click to view

sub get_options {
  my ($db_name, $db_host, $db_user, $db_pass, $db_port, $help, @species, $group, $release, $dir, $compress, $ucsc);
  $db_port = 3306;
  $group = 'core';
  $compress=0;
  $ucsc=0;
  $dir = File::Spec->curdir();
 
  GetOptions(
    "db_name|dbname|database=s"       => \$db_name,
    "db_host|dbhost|host=s"           => \$db_host,
    "db_user|dbuser|user|username=s"  => \$db_user,
    "db_pass|dbpass|pass|password=s"  => \$db_pass,
    "db_port|dbport|port=s"           => \$db_port,
    "species=s@"                      => \@species,
    "version|release=i"               => \$release,
    "directory|dir=s"                 => \$dir,
    "compress|gzip|gz!"               => \$compress,
    "ucsc!"                           => \$ucsc,
    "help|h!"                         => \$help,
  );
 
  my %args = (
    -HOST => $db_host, -PORT => $db_port, 
    -USER => $db_user
  );
  $args{-PASS} = $db_pass if $db_pass;
 
  # must be a DB we need
  if($db_name) {
    Bio::EnsEMBL::DBSQL::DBAdaptor->new(
      %args,
      -DBNAME => $db_name,
      -SPECIES => $db_name,
      -GROUP => 'core'
    );
  }
  else {
    Bio::EnsEMBL::Registry->load_registry_from_db(
      %args,
      -DB_VERSION => $release,
    );
  }
  
  $COMPRESS = $compress;
  $UCSC = $ucsc;
  my @dbas;
  my %final_dbas;
    
  if(@species) {
    say "Working against a restricted species list";
    foreach my $s (@species) {
      my $dba = Bio::EnsEMBL::Registry->get_DBAdaptor($s, $group);
      die "Cannot find a DBAdaptor for the species ${s}" unless $dba;
      push(@dbas, $dba);
    }
  }
  else {
    say "Dumping chain file for all available species";
    @dbas = @{Bio::EnsEMBL::Registry->get_all_DBAdaptors(-GROUP => 'core')};
  }
  
  foreach my $dba (@dbas) {
    my $entries = get_liftover_mappings($dba);
    my $prod_name = $dba->get_MetaContainer()->get_production_name();
    $prod_name //= $dba->species();
    $final_dbas{$prod_name} = $dba if @{$entries};
    $dba->dbc->disconnect_if_idle();
  }
  
  # Return including sorting the species by name as it's just nicer that way
  return ($dir, map { $final_dbas{$_} } sort keys %final_dbas);
}

◆ get_reverse_assembly_mappings()

public get_reverse_assembly_mappings ( )

Undocumented method

Code:

click to view

◆ get_sql()

public get_sql ( )

Undocumented method

Code:

click to view

◆ print_chains()

public print_chains ( )

Undocumented method

Code:

click to view

◆ run()

public run ( )

Undocumented method

Code:

click to view

Undocumented method

Code:

click to view

Undocumented method

Code:

click to view

Undocumented method

Code:

click to view

◆ run_on_dba()

public run_on_dba ( )

Undocumented method

Code:

click to view

◆ write_mappings()

public write_mappings ( )

Undocumented method

Code:

click to view

Functions

Function Documentation

◆ build_chain_mappings()

◆ create_ucsc_chains()

◆ ensembl_to_ucsc_name()

◆ get_assembly_mappings()

◆ get_liftover_mappings()

◆ get_options()

◆ get_reverse_assembly_mappings()

◆ get_sql()

◆ print_chains()

◆ run()

◆ run_on_dba()

◆ write_mappings()