<commit_before>use strict;
use warnings;

use Getopt::Long;
use Fcntl ':flock';
use DBI;
use DBH;
use FindBin qw( $Bin );

use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Utils::Exception qw(warning throw verbose);
use Bio::EnsEMBL::Utils::Sequence qw(reverse_comp expand);
use Bio::EnsEMBL::Variation::ConsequenceType;
use Bio::EnsEMBL::Utils::TranscriptAlleles qw(type_variation);  #function to calculate the consequence type of a Variation in a transcript
use ImportUtils qw(dumpSQL debug load create_and_load);

my ($TMP_DIR, $TMP_FILE, $LIMIT,$status_file);

{

  my ($species,$limit,$num_processes);

  GetOptions('species=s'   => \$species,
	     'tmpdir=s'  => \$ImportUtils::TMP_DIR,
	     'tmpfile=s' => \$ImportUtils::TMP_FILE,
	     'limit=s'   => \$limit,
 	     'num_processes=i' => \$num_processes,
 	     'status_file=s' => \$status_file);

  $num_processes ||= 1;

  $LIMIT = ($limit) ? " $limit " : ''; #the limit will refer to the slices the process must used (1,4-5,7,.....)

  usage('-num_processes must at least be 1') if ($num_processes == 0);
  usage('-status_file argument is required') if (!$status_file);

  my $registry_file ||= $Bin . "/ensembl.registry";

  Bio::EnsEMBL::Registry->load_all( $registry_file );

  my $cdba = Bio::EnsEMBL::Registry->get_DBAdaptor($species,'core');
  my $vdba = Bio::EnsEMBL::Registry->get_DBAdaptor($species,'variation');
  my $fdba = Bio::EnsEMBL::Registry->get_DBAdaptor($species,'funcgen');

  my $dbVar = $vdba->dbc;
  my $dbCore = $cdba;
  my $dbFunc = $fdba;

  $TMP_DIR  = $ImportUtils::TMP_DIR;
  $TMP_FILE = $ImportUtils::TMP_FILE;


  transcript_variation($dbCore, $dbVar, $dbFunc);

  open STATUS, ">>$TMP_DIR/$status_file"
    or throw("Could not open tmp file: $TMP_DIR/$status_file\n"); 
  flock(STATUS,LOCK_EX);
  seek(STATUS, 0, 2); #move to the end of the file
  print STATUS "process finished\n";
  flock(STATUS,LOCK_UN);
  close STATUS;
  #check if it is the last process
  my $processes = `cat $TMP_DIR/$status_file | wc -l`;
  if ($processes == $num_processes){
      #if is the last process, finish it
    debug("Last process: ready to import data");
    last_process($dbCore,$dbVar);
  }
}

#
# Loads the transcript variation table.  Retrieves every transcript in the
# the database and types all of the snps in the vicinity of the transcript.
# The amino acid changes for coding snps is also set.
#
#

sub transcript_variation {
  my $dbCore = shift;
  my $dbVar  = shift;
  my $dbFunc = shift;

  my $UPSTREAM = 5000;
  my $DNSTREAM = 5000;
  my $MAX_FEATURE_LENGTH  = 1000;
 
  my $sth = $dbVar->prepare
    (qq{SELECT vf.variation_feature_id, vf.seq_region_start, vf.seq_region_end,
               vf.seq_region_strand, vf.allele_string
        FROM   variation_feature vf
        WHERE  vf.seq_region_id = ?
        AND    vf.seq_region_end >= ?
        AND    vf.seq_region_start <= ?
        AND    vf.seq_region_start >= ?
	});

  my $sa = $dbCore->get_SliceAdaptor();

  my $dbname = $dbVar->dbname(); #get the name of the database to create the file
  my $host = `hostname`;
  chop $host;
  open FH, ">$TMP_DIR/$dbname.transcript_variation_$host\:$$\.txt";
  #open FH, ">/tmp/$dbname.transcript_variation_$host\_$$\.txt" or die "Could not open file: $!\n";
  my $inc_non_ref = 1;

  my $slices = $sa->fetch_all('toplevel', undef, $inc_non_ref);
  #order the slices
  my @slices_ordered = sort {$a->seq_region_name cmp $b->seq_region_name }  @{$slices};

  my ($offset,$length) = split /,/,$LIMIT; #get the offset and length of the elements we have to get from the slice
  # assumes that variation features have already been pushed to toplevel

  my $rfa = $dbFunc->get_ExternalFeatureAdaptor();
  my $gene_adaptor = $dbCore->get_GeneAdaptor();
  my $transcript_adaptor = $dbCore->get_TranscriptAdaptor();
					      
#  my $rfa = $dbCore->get_RegulatoryFeatureAdaptor(); this information comes from the FuncGen database
  my %done;
  foreach my $slice (splice (@slices_ordered,$offset,$length)) {
    debug("Processing transcript variations for ",
          $slice->seq_region_name(), "\n");
    #first consider the regulatory region that overlap with SNPs
    foreach my $rf (@{$rfa->fetch_all_by_Slice($slice)}) {
	next if (($rf->feature_set->name !~ /miRNA/) && ($rf->feature_set->name !~ /cisRED/));
#    foreach my $rf (@{fa->fetch_all_by_Slice($slice)}) {
      # request all variations which lie in the region of a regulate feature
      $sth->execute($slice->get_seq_region_id(),
                    $rf->seq_region_start(),
                    $rf->seq_region_end(),
                    $rf->seq_region_start());
      my $rows = $sth->fetchall_arrayref();
      my ($start,$end, $rf_start, $rf_end); #start, end of the variation feature in the slice

      foreach my $row (@$rows) {
	$start = $row->[1];
	$end = $row->[2];
	$rf_start = $rf->seq_region_start();
	$rf_end = $rf->seq_region_end();
	#print "start is $start,end is $end,rf_start is $rf_start and rf_end is $rf_end\n";
	if ($end >= $rf_start and $start <= $rf_end) {
#	  foreach my $g (@{$rf->regulated_genes()}) {
	  foreach my $dbEntry (@{$rf->get_all_DBEntries('core_gene','MISC')}) {
	      my $g = $gene_adaptor->fetch_by_stable_id($dbEntry->primary_id); #get the gene for the stable_id
	      next if(!defined $g); #some of the genes do not seem to be in the core database
	      my $g_start = $g->seq_region_start;
	      my $g_end = $g->seq_region_end;
	    #SNPs needs to be out side gene+flanking region,otherwise will be considered later
	    if ($end < $g_start - $UPSTREAM or $start > $g_end + $DNSTREAM ) {
	      foreach my $tr (@{$g->get_all_Transcripts()}) {
		my @arr = ($tr->dbID,$row->[0],'\N','\N','\N','\N','\N','REGULATORY_REGION');
		if (! $done{$row->[0]}{$tr->dbID}) {#get rid of duplicated transcript entries
		  print FH join("\t", @arr), "\n";
		  $done{$row->[0]}{$tr->dbID}=1;
		}
	      }
	    }
	  }
#	  foreach my $tr (@{$rf->regulated_transcripts()}) {
	  foreach my $dbEntry (@{$rf->get_all_DBEntries('core_transcript')}) {	     
	      my $tr = $transcript_adaptor->fetch_by_stable_id($dbEntry->primary_id); #get transcript for stable_id
	      next if(!defined $tr); #some of the transcripts do not seem to be in the core database
	    my $tr_start = $tr->seq_region_start;
	    my $tr_end = $tr->seq_region_end;
	    #SNPs needs to be out side gene+flanking region,otherwise will be considered later
	    if ($end < $tr_start - $UPSTREAM or $start > $tr_end + $DNSTREAM ) {
	      my @arr = ($tr->dbID,$row->[0],'\N','\N','\N','\N','\N','REGULATORY_REGION');
	      if (! $done{$row->[0]}{$tr->dbID}) {#get rid of duplicated transcript entries
		print FH join("\t", @arr), "\n";
		$done{$row->[0]}{$tr->dbID}=1;
	      }
	    }
	  }
	}
      }<commit_msg>only for human to calculate regulatory feature consequence<commit_after>use strict;
use warnings;

use Getopt::Long;
use Fcntl ':flock';
use DBI;
use DBH;
use FindBin qw( $Bin );

use Bio::EnsEMBL::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Utils::Exception qw(warning throw verbose);
use Bio::EnsEMBL::Utils::Sequence qw(reverse_comp expand);
use Bio::EnsEMBL::Variation::ConsequenceType;
use Bio::EnsEMBL::Utils::TranscriptAlleles qw(type_variation);  #function to calculate the consequence type of a Variation in a transcript
use ImportUtils qw(dumpSQL debug load create_and_load);

my ($TMP_DIR, $TMP_FILE, $LIMIT,$status_file);

{

  my ($species,$limit,$num_processes);

  GetOptions('species=s'   => \$species,
	     'tmpdir=s'  => \$ImportUtils::TMP_DIR,
	     'tmpfile=s' => \$ImportUtils::TMP_FILE,
	     'limit=s'   => \$limit,
 	     'num_processes=i' => \$num_processes,
 	     'status_file=s' => \$status_file);

  $num_processes ||= 1;

  $LIMIT = ($limit) ? " $limit " : ''; #the limit will refer to the slices the process must used (1,4-5,7,.....)

  usage('-num_processes must at least be 1') if ($num_processes == 0);
  usage('-status_file argument is required') if (!$status_file);

  my $registry_file ||= $Bin . "/ensembl.registry";

  Bio::EnsEMBL::Registry->load_all( $registry_file );

  my $cdba = Bio::EnsEMBL::Registry->get_DBAdaptor($species,'core');
  my $vdba = Bio::EnsEMBL::Registry->get_DBAdaptor($species,'variation');
  my $fdba = Bio::EnsEMBL::Registry->get_DBAdaptor($species,'funcgen') if $species =~ /hum|hom/i;

  my $dbVar = $vdba->dbc;
  my $dbCore = $cdba;
  my $dbFunc = $fdba;

  $TMP_DIR  = $ImportUtils::TMP_DIR;
  $TMP_FILE = $ImportUtils::TMP_FILE;


  transcript_variation($dbCore, $dbVar, $dbFunc);

  open STATUS, ">>$TMP_DIR/$status_file"
    or throw("Could not open tmp file: $TMP_DIR/$status_file\n"); 
  flock(STATUS,LOCK_EX);
  seek(STATUS, 0, 2); #move to the end of the file
  print STATUS "process finished\n";
  flock(STATUS,LOCK_UN);
  close STATUS;
  #check if it is the last process
  my $processes = `cat $TMP_DIR/$status_file | wc -l`;
  if ($processes == $num_processes){
      #if is the last process, finish it
    debug("Last process: ready to import data");
    last_process($dbCore,$dbVar);
  }
}

#
# Loads the transcript variation table.  Retrieves every transcript in the
# the database and types all of the snps in the vicinity of the transcript.
# The amino acid changes for coding snps is also set.
#
#

sub transcript_variation {
  my $dbCore = shift;
  my $dbVar  = shift;
  my $dbFunc = shift;

  my $UPSTREAM = 5000;
  my $DNSTREAM = 5000;
  my $MAX_FEATURE_LENGTH  = 1000;
 
  my $sth = $dbVar->prepare
    (qq{SELECT vf.variation_feature_id, vf.seq_region_start, vf.seq_region_end,
               vf.seq_region_strand, vf.allele_string
        FROM   variation_feature vf
        WHERE  vf.seq_region_id = ?
        AND    vf.seq_region_end >= ?
        AND    vf.seq_region_start <= ?
        AND    vf.seq_region_start >= ?
	});

  my $sa = $dbCore->get_SliceAdaptor();

  my $dbname = $dbVar->dbname(); #get the name of the database to create the file
  my $host = `hostname`;
  chop $host;
  open FH, ">$TMP_DIR/$dbname.transcript_variation_$host\:$$\.txt";
  #open FH, ">/tmp/$dbname.transcript_variation_$host\_$$\.txt" or die "Could not open file: $!\n";
  my $inc_non_ref = 1;

  my $slices = $sa->fetch_all('toplevel', undef, $inc_non_ref);
  #order the slices
  my @slices_ordered = sort {$a->seq_region_name cmp $b->seq_region_name }  @{$slices};

  my ($offset,$length) = split /,/,$LIMIT; #get the offset and length of the elements we have to get from the slice
  # assumes that variation features have already been pushed to toplevel


  my %done;
  foreach my $slice (splice (@slices_ordered,$offset,$length)) {
    debug("Processing transcript variations for ",
	  $slice->seq_region_name(), "\n");
    #if it is human, $dbFunc is defined, therefore go through below to calculate regulatory feature
    if ($dbFunc) {
      my $rfa = $dbFunc->get_ExternalFeatureAdaptor();
      my $gene_adaptor = $dbCore->get_GeneAdaptor();
      my $transcript_adaptor = $dbCore->get_TranscriptAdaptor();
      #first consider the regulatory region that overlap with SNPs
      foreach my $rf (@{$rfa->fetch_all_by_Slice($slice)}) {
	next if (($rf->feature_set->name !~ /miRNA/) && ($rf->feature_set->name !~ /cisRED/));

	# request all variations which lie in the region of a regulate feature
	$sth->execute($slice->get_seq_region_id(),
		      $rf->seq_region_start(),
		      $rf->seq_region_end(),
		      $rf->seq_region_start());
	my $rows = $sth->fetchall_arrayref();
	my ($start,$end, $rf_start, $rf_end); #start, end of the variation feature in the slice

	foreach my $row (@$rows) {
	  $start = $row->[1];
	  $end = $row->[2];
	  $rf_start = $rf->seq_region_start();
	  $rf_end = $rf->seq_region_end();
	  #print "start is $start,end is $end,rf_start is $rf_start and rf_end is $rf_end\n";
	  if ($end >= $rf_start and $start <= $rf_end) {
	    foreach my $dbEntry (@{$rf->get_all_DBEntries('core_gene','MISC')}) {
	      my $g = $gene_adaptor->fetch_by_stable_id($dbEntry->primary_id); #get the gene for the stable_id
	      next if(!defined $g); #some of the genes do not seem to be in the core database
	      my $g_start = $g->seq_region_start;
	      my $g_end = $g->seq_region_end;
	      #SNPs needs to be out side gene+flanking region,otherwise will be considered later
	      if ($end < $g_start - $UPSTREAM or $start > $g_end + $DNSTREAM ) {
		foreach my $tr (@{$g->get_all_Transcripts()}) {
		  my @arr = ($tr->dbID,$row->[0],'\N','\N','\N','\N','\N','REGULATORY_REGION');
		  if (! $done{$row->[0]}{$tr->dbID}) {#get rid of duplicated transcript entries
		    print FH join("\t", @arr), "\n";
		    $done{$row->[0]}{$tr->dbID}=1;
		  }
		}
	      }
	    }
	    foreach my $dbEntry (@{$rf->get_all_DBEntries('core_transcript')}) {	     
	      my $tr = $transcript_adaptor->fetch_by_stable_id($dbEntry->primary_id); #get transcript for stable_id
	      next if(!defined $tr); #some of the transcripts do not seem to be in the core database
	      my $tr_start = $tr->seq_region_start;
	      my $tr_end = $tr->seq_region_end;
	      #SNPs needs to be out side gene+flanking region,otherwise will be considered later
	      if ($end < $tr_start - $UPSTREAM or $start > $tr_end + $DNSTREAM ) {
		my @arr = ($tr->dbID,$row->[0],'\N','\N','\N','\N','\N','REGULATORY_REGION');
		if (! $done{$row->[0]}{$tr->dbID}) {#get rid of duplicated transcript entries
		  print FH join("\t", @arr), "\n";
		  $done{$row->[0]}{$tr->dbID}=1;
		}
	      }
	    }
	  }
	}
      }
{'fraction_non_alphanumeric': 0.15944295501936254, 'fraction_numerical': 0.004766160262138814, 'mean_word_length': 3.750265298903431, 'pattern_counts': {'":': 0, '<': 11, '<?xml version=': 0, '>': 143, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 0, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '3306215', 'n_tokens_mistral': 5074, 'n_tokens_neox': 4770, 'n_words': 1237}