<commit_before># Ensembl module for Bio::EnsEMBL::Analysis::Runnable::Blast
#
# Copyright (c) 2004 Ensembl
#

=head1 NAME

  Bio::EnsEMBL::Analysis::Runnable::BlastRfam

=head1 SYNOPSIS

  my $blast = Bio::EnsEMBL::Analysis::Runnable::BlastRfam->
  new(
      -query => $slice,
      -program => 'wublastn',
      -database => 'embl_vertrna',
      -options => 'hitdist=40 -cpus=1',
      -parser => $bplitewrapper,
      -filter => $featurefilter,
     );
  $blast->run;
  my @output =@{$blast->output};

=head1 DESCRIPTION

Modified blast runnable for specific use with RFAMSEQ.
Use for running BLASTN of genomic vs RFAMSEQ prior to 
ncRNA analysis using Infernal.
Keeps the coverage in the dna_align_feature score field.
Also clusters overlapping hits and picks the one with the lowest 
evalue to represent that cluster.

=head1 CONTACT

Post questions to the Ensembl development list: <EMAIL>


=cut


package Bio::EnsEMBL::Analysis::Runnable::BlastRfam;

use strict;
use warnings;

use Bio::EnsEMBL::Analysis::Runnable;
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
use Bio::EnsEMBL::Utils::Argument qw( rearrange );
use Bio::EnsEMBL::Analysis::Runnable::Blast;
use vars qw(@ISA);

@ISA = qw(Bio::EnsEMBL::Analysis::Runnable::Blast);

=head2 run

  Arg [1]   : Bio::EnsEMBL::Analysis::Runnable
  Arg [2]   : string, directory
  Function  : a generic run method. This checks the directory specifed
  to run it, write the query sequence to file, marks the query sequence
  file and results file for deletion, runs the analysis parses the 
  results and deletes any files
  Returntype: 1
  Exceptions: throws if no query sequence is specified
  Example   : 

=cut


sub run{
  my ($self, $dir) = @_;
  my $coverage = 70; # Coverage filter for very common blast hits
  $self->workdir($dir) if($dir);
  throw("Can't run ".$self." without a query sequence") 
    unless($self->query);
  $self->checkdir();
  my $filename = $self->write_seq_file();
  $self->files_to_delete($filename);
  $self->files_to_delete($self->resultsfile);
  $self->run_ncbi_analysis("-W7 -F F -b 1000000 -v 1000000","/ecs2/work2/sw4/RFAM/Partition/low_copy.fasta");  
  $self->parse_results(20);
  # evil way of emptying the result array!
  $self->{'results_files'} = ();
  $self->database("/ecs2/work2/sw4/RFAM/Partition/high_copy.fasta");
  # NOTE TO SELF need to set the analysis to use wublastn or it will have a benny unless you want to hard code it?
  $self->run_analysis;
  $self->parse_results($coverage);
  $self->delete_files;
  return 1;
}


=head2 parse_results

  Arg [1]   : Bio::EnsEMBL::Analysis::Runnable::Blast
  Function  : override results parsing to allow for 
coverage calculations
  Returntype: none
  Exceptions: none
  Example   : 

=cut


sub parse_results{
  my ($self,$coverage_cutoff) = @_;
  my $results = $self->results_files;
  my @daf_coverage_results;
  my $filtered_output;
  my $bplite = $self->parser->get_parsers($results);
  foreach my $blast (@{$bplite}){
      while( my $subject = $blast->nextSbjct){
	 while (my $hsp = $subject->nextHSP) {
	   my @daf_results;
	   my $hsp_length = $hsp->length."\t";
	   my $subject_length = $subject->{'LENGTH'};
	   $subject_length = 1 unless ($subject_length);
	   my $coverage = $hsp_length/$subject_length*100;
	   $coverage =~ s/\.\d+//;
	   if ($coverage_cutoff){
	     next unless($coverage > $coverage_cutoff);
	   }
	   $subject->name =~ /^(\S+)\/\S+\s+(\w+);\w+/;
	   my $name = $2."-".$1;
	   push @daf_results, $self->parser->split_hsp($hsp,$name);
	   # add coverage into daf score?
	   foreach my $daf(@daf_results){
	     # swaps strands over if blast aligns +strand genomic to -ve strand in 
	     # RFAM file, Thing is RFAM file sequences are the correct orientation
	     # Whereas the genomic can be either, with unspliced DNA it becomes
	     # impossibe to tell automatically
	     if ($daf->hstrand == -1 && $daf->strand == 1){
	       $daf->strand(-1);
	       $daf->hstrand(1);
	     }
	     $daf->score($coverage);
	     push  @daf_coverage_results, $daf;
	   }
	 }
      }
    }
#  print "Before clustering ".scalar( @daf_coverage_results)."\n";
  return undef unless  ( @daf_coverage_results);
  my $output = $self->cluster(\@daf_coverage_results);
  $self->output($output);
}

# foreach hit, look at all the other hits and push into an array anything that overlaps with it
# then sort the array and take the top scoring hits for heach familly that lives in it.
# you have 3 instances of a variable called RFAM in the same method all of which hold different values

sub cluster{
  my ($self,$dafs_ref)=@_;
  my @dafs = @$dafs_ref;
  my $start =0;
  my @representative_sequences;
 DAFS: foreach my $daf (@dafs){
    my %family_cluster;
    $start ++;
    next DAFS unless($daf);	
    my $RFAM = substr($daf->hseqname,0,7);
    push @{$family_cluster{$RFAM}},$daf;
  MATCHES:  for (my $index = $start; $index <= $#dafs ; $index ++){
      next MATCHES unless ($dafs[$index]);
      if ($daf->overlaps($dafs[$index])){
	my $familly = substr($dafs[$index]->hseqname,0,7);
	push @{$family_cluster{$familly}},$dafs[$index];
	$dafs[$index] = undef;
      }
    }
    foreach my $domain (keys %family_cluster){
      my @temp_array = @{$family_cluster{$domain}};
      # get highest scoring alignment for each family to be representative
      @temp_array = sort{$a->p_value <=> $b->p_value} @temp_array;
      push @representative_sequences,shift @temp_array;
    }
  }
  return \@representative_sequences;
}

#####################################################
# USE THIS TO RUN NCBI BLAST IF YOU ARE SO INCLINED #
#####################################################
<commit_msg>Removed debugging statements, Removed storing coverage in dna_align_feature score field<commit_after># Ensembl module for Bio::EnsEMBL::Analysis::Runnable::Blast
#
# Copyright (c) 2004 Ensembl
#

=head1 NAME

  Bio::EnsEMBL::Analysis::Runnable::BlastRfam

=head1 SYNOPSIS

  my $blast = Bio::EnsEMBL::Analysis::Runnable::BlastRfam->
  new(
      -query => $slice,
      -program => 'wublastn',
      -database => 'embl_vertrna',
      -options => 'hitdist=40 -cpus=1',
      -parser => $bplitewrapper,
      -filter => $featurefilter,
     );
  $blast->run;
  my @output =@{$blast->output};

=head1 DESCRIPTION

Modified blast runnable for specific use with RFAMSEQ.
Use for running BLASTN of genomic vs RFAMSEQ prior to 
ncRNA analysis using Infernal.
Keeps the coverage in the dna_align_feature score field.
Also clusters overlapping hits and picks the one with the lowest 
evalue to represent that cluster.

=head1 CONTACT

Post questions to the Ensembl development list: ensembl-<EMAIL>


=cut


package Bio::EnsEMBL::Analysis::Runnable::BlastRfam;

use strict;
use warnings;

use Bio::EnsEMBL::Analysis::Runnable;
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
use Bio::EnsEMBL::Utils::Argument qw( rearrange );
use Bio::EnsEMBL::Analysis::Runnable::Blast;
use vars qw(@ISA);

@ISA = qw(Bio::EnsEMBL::Analysis::Runnable::Blast);

=head2 run

  Arg [1]   : Bio::EnsEMBL::Analysis::Runnable
  Arg [2]   : string, directory
  Function  : a generic run method. This checks the directory specifed
  to run it, write the query sequence to file, marks the query sequence
  file and results file for deletion, runs the analysis parses the 
  results and deletes any files
  Returntype: 1
  Exceptions: throws if no query sequence is specified
  Example   : 

=cut


sub run{
  my ($self, $dir) = @_;
  my $coverage = 70; # Coverage filter for very common blast hits
  $self->workdir($dir) if($dir);
  throw("Can't run ".$self." without a query sequence") 
    unless($self->query);
  $self->checkdir();
  my $filename = $self->write_seq_file();
  $self->files_to_delete($filename);
  $self->files_to_delete($self->resultsfile);
  $self->run_ncbi_analysis("-W7 -F F -b 1000000 -v 1000000","/ecs2/work2/sw4/RFAM/Partition/low_copy.fasta");  
  $self->parse_results(20);
  $self->{'results_files'} = ();
  $self->database("/ecs2/work2/sw4/RFAM/Partition/high_copy.fasta");
  $self->run_analysis;
  $self->parse_results($coverage);
  $self->delete_files;
  return 1;
}


=head2 parse_results

  Arg [1]   : Bio::EnsEMBL::Analysis::Runnable::Blast
  Function  : override results parsing to allow for 
coverage calculations
  Returntype: none
  Exceptions: none
  Example   : 

=cut


sub parse_results{
  my ($self,$coverage_cutoff) = @_;
  my $results = $self->results_files;
  my @daf_coverage_results;
  my $filtered_output;
  my $bplite = $self->parser->get_parsers($results);
  foreach my $blast (@{$bplite}){
      while( my $subject = $blast->nextSbjct){
	 while (my $hsp = $subject->nextHSP) {
	   my @daf_results;
	   my $hsp_length = $hsp->length."\t";
	   my $subject_length = $subject->{'LENGTH'};
	   $subject_length = 1 unless ($subject_length);
	   my $coverage = $hsp_length/$subject_length*100;
	   $coverage =~ s/\.\d+//;
	   if ($coverage_cutoff){
	     next unless($coverage > $coverage_cutoff);
	   }
	   $subject->name =~ /^(\S+)\/\S+\s+(\w+);\w+/;
	   my $name = $2."-".$1;
	   push @daf_results, $self->parser->split_hsp($hsp,$name);
	   # add coverage into daf score?
	   foreach my $daf(@daf_results){
	     # swaps strands over if blast aligns +strand genomic to -ve strand in 
	     # RFAM file,  RFAM file sequences are the correct orientation
	     # Whereas the genomic can be either, with unspliced DNA it becomes
	     # impossibe to tell automatically
	     if ($daf->hstrand == -1 && $daf->strand == 1){
	       $daf->strand(-1);
	       $daf->hstrand(1);
	     }
#	     $daf->score($coverage);
	     push  @daf_coverage_results, $daf;
	   }
	 }
      }
    }
  return undef unless  ( @daf_coverage_results);
  my $output = $self->cluster(\@daf_coverage_results);
  $self->output($output);
}

# foreach hit, look at all the other hits and push into an array anything that overlaps with it
# then sort the array and take the top scoring hits for heach familly that lives in it.

sub cluster{
  my ($self,$dafs_ref)=@_;
  my @dafs = @$dafs_ref;
  @dafs = sort{$a->p_value <=> $b->p_value} @dafs;
  my $start =0;
  my @representative_sequences;    
 DAFS: foreach my $daf (@dafs){
    $start ++;
    my %family_cluster;
    next DAFS unless($daf);	
    my $RFAM = substr($daf->hseqname,0,7);
    push @{$family_cluster{$RFAM}},$daf;
  MATCHES:  for (my $index = $start; $index <= $#dafs ; $index ++){
      next MATCHES unless ($dafs[$index]);
      if ($daf->overlaps($dafs[$index])){
	my $familly = substr($dafs[$index]->hseqname,0,7);
	push @{$family_cluster{$familly}},$dafs[$index];
	$dafs[$index] = undef;
      }
    }
    foreach my $domain (keys %family_cluster){
      my @temp_array = @{$family_cluster{$domain}};
      # get highest scoring alignment for each family to be representative
      @temp_array = sort{$a->p_value <=> $b->p_value} @temp_array;
      push @representative_sequences,shift @temp_array;
    }
  }
  return \@representative_sequences;
}

#####################################################
# USE THIS TO RUN NCBI BLAST IF YOU ARE SO INCLINED #
#####################################################

{'fraction_non_alphanumeric': 0.1582210242587601, 'fraction_numerical': 0.010512129380053909, 'mean_word_length': 3.8841597191750767, 'pattern_counts': {'":': 0, '<': 10, '<?xml version=': 0, '>': 102, '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': '6882689', 'n_tokens_mistral': 4068, 'n_tokens_neox': 3758, 'n_words': 1038}