<commit_before>#!/usr/bin/env perl
# Copyright [1999-2014] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.




=head1 CONTACT

  Please email comments or questions to the public Ensembl
  developers list at <http://lists.ensembl.org/mailman/listinfo/dev>.

  Questions may also be sent to the Ensembl help desk at
  <helpdesk.org>.

=cut
package Bio::EnsEMBL::Variation::Pipeline::Remapping::InitMapping;

use strict;
use warnings;

use FileHandle;
use Bio::DB::Fasta;
use Bio::EnsEMBL::Registry;

use base ('Bio::EnsEMBL::Hive::Process');

sub fetch_input {
    my $self = shift;
    my $registry = 'Bio::EnsEMBL::Registry';
    $registry->load_all($self->param('registry_file'));
    my $cdba = $registry->get_DBAdaptor($self->param('species'), 'core');
    my $vdba = $registry->get_DBAdaptor($self->param('species'), 'variation');

    $self->param('registry', $registry);
    $self->param('cdba', $cdba);
    $self->param('vdba', $vdba);
}

sub run {
    my $self = shift;
    if ($self->param('generate_fasta_files')) {
        $self->dump_features();
        $self->generate_mapping_input();
    }
}

sub write_output {
    my $self = shift;
    # initialise mapping jobs
    my $file_count      = $self->param('file_count');
    my $fasta_files_dir = $self->param('fasta_files_dir');
    my $bam_files_dir   = $self->param('bam_files_dir');
    my @jobs;
    my $i = 1;
    while ($i <= $file_count) {
        push @jobs, {
            'file_number'   => $i,
            'bam_files_dir' => $bam_files_dir,
            'fasta_file'    => "$fasta_files_dir/$i.fa",
            'sam_file'      => "$bam_files_dir/$i.sam",
            'bam_file'      => "$bam_files_dir/$i.bam",
            'err_file'      => "$bam_files_dir/$i.err",
            'out_file'      => "$bam_files_dir/$i.out",
        };
        $i++;
    }
    $self->dataflow_output_id(\@jobs, 2);
}

sub generate_mapping_input {
    my $self = shift;

    my $old_assembly_fasta_file_dir = $self->param('old_assembly_fasta_file_dir');
    my $fasta_db = Bio::DB::Fasta->new($old_assembly_fasta_file_dir, -reindex => 1);
    $self->param('fasta_db', $fasta_db);

    # store end-coordinates for all seq_regions to check that variation_location + flank_seq_length < slice_end
    my $cdba = $self->param('cdba');
    my $sa = $cdba->get_SliceAdaptor;
    # don't include asm exceptions
    my $slices = $sa->fetch_all('toplevel', undef, 0);
    my $seq_regions = {};
    foreach my $slice (@$slices) {
        my $end             = $slice->end;
        my $seq_region_name = $slice->seq_region_name;
        if ($end > $seq_regions->{$seq_region_name}->{end}) {
            $seq_regions->{$seq_region_name}->{end} = $end;
        }
    }
    foreach my $seq_region_name (keys %$seq_regions) {
        my $end = $seq_regions->{$seq_region_name}->{end};
    }

    $self->param('seq_regions', $seq_regions);

    my $dump_features_dir = $self->param('dump_features_dir');
    my $fasta_files_dir   = $self->param('fasta_files_dir');
    my $pipeline_dir      = $self->param('pipeline_dir');
    my $fh_allele_length = FileHandle->new("$pipeline_dir/qc_allele_string_length.txt", 'w');
    my $fh_ref_seq       = FileHandle->new("$pipeline_dir/qc_ref_seq.txt", 'w');

    $self->param('qc_allele_string_length', $fh_allele_length);
    $self->param('qc_ref_seq', $fh_ref_seq);

    my $variants_with_multi_map = {};

    opendir(DIR, $dump_features_dir) or die $!;
    while (my $file = readdir(DIR)) {
        if ($file =~ /^(.+)\.txt$/) {
            my $file_number = $1;
            my $fh = FileHandle->new("$dump_features_dir/$file", 'r');
            my $fh_fasta_file = FileHandle->new("$fasta_files_dir/$file_number.fa", 'w');
            while (<$fh>) {
                chomp;
                my $data = $self->read_line($_);
                my $seq_region_name = $data->{seq_region_name},
                my $feature_id      = $data->{variation_feature_id};
                my $start           = $data->{seq_region_start};
                my $end             = $data->{seq_region_end};
                my $strand          = $data->{seq_region_strand};
                my $allele_string   = $data->{allele_string};
                my $map_weight      = $data->{map_weight};
                my $variation_name  = $data->{variation_name};
                if ($map_weight > 1) {
                    $variants_with_multi_map->{$variation_name}++;
                    next;
                }
                my ($flank_start, $upstream_flank_length, $downstream_flank_length, $flank_end, $variant_length) = @{$self->flank_coordinates($seq_region_name, $start, $end, $strand)};

                # variant surrounded by flank sequences
                my $query_sequence = $self->get_query_sequence($seq_region_name, $flank_start, $flank_end, $strand);

                # replace empty space with underscores <_>
                $allele_string =~ s/\s/_/g;
                $self->qc_alleles($query_sequence, $upstream_flank_length, $variant_length, $allele_string, $variation_name);
#                $allele_string = $self->qc_alleles($config, $query_sequence, $upstream_flank_length, $variant_length, $allele_string, $variation_name);

                my $name = "$seq_region_name:$start:$end:$strand:$allele_string:$variation_name";
                my $id = ">$feature_id-$upstream_flank_length-$variant_length-$downstream_flank_length-$name";
                print $fh_fasta_file "$id\n$query_sequence\n";
            } # end while (read feature file for seq_region)
            $fh->close();
            $fh_fasta_file->close();
        }
    }

    foreach my $file_type (qw/qc_allele_string_length qc_ref_seq/) {
        my $fh = $self->param($file_type);
        $fh->close();
    }

    # store multi map
    my $fh_qc_multi_map = FileHandle->new("$pipeline_dir/multi_map.txt", 'w');

    foreach my $variation_name (keys %$variants_with_multi_map) {
        my $count = $variants_with_multi_map->{$variation_name};
        print $fh_qc_multi_map "$variation_name\t$count\n";
    }
    $fh_qc_multi_map->close();
}<commit_msg>Remove allele_string from query_name in fasta file<commit_after>#!/usr/bin/env perl
# Copyright [1999-2014] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
#
# Licensed under the Apache License, Version 2.0 (the "License");
# you may not use this file except in compliance with the License.
# You may obtain a copy of the License at
#
#      http://www.apache.org/licenses/LICENSE-2.0
#
# Unless required by applicable law or agreed to in writing, software
# distributed under the License is distributed on an "AS IS" BASIS,
# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
# See the License for the specific language governing permissions and
# limitations under the License.




=head1 CONTACT

  Please email comments or questions to the public Ensembl
  developers list at <http://lists.ensembl.org/mailman/listinfo/dev>.

  Questions may also be sent to the Ensembl help desk at
  <helpdesk.org>.

=cut
package Bio::EnsEMBL::Variation::Pipeline::Remapping::InitMapping;

use strict;
use warnings;

use FileHandle;
use Bio::DB::Fasta;
use Bio::EnsEMBL::Registry;

use base ('Bio::EnsEMBL::Hive::Process');

sub fetch_input {
    my $self = shift;
    my $registry = 'Bio::EnsEMBL::Registry';
    $registry->load_all($self->param('registry_file'));
    my $cdba = $registry->get_DBAdaptor($self->param('species'), 'core');
    my $vdba = $registry->get_DBAdaptor($self->param('species'), 'variation');

    $self->param('registry', $registry);
    $self->param('cdba', $cdba);
    $self->param('vdba', $vdba);
}

sub run {
    my $self = shift;
    if ($self->param('generate_fasta_files')) {
        $self->dump_features();
        $self->generate_mapping_input();
    }
}

sub write_output {
    my $self = shift;
    # initialise mapping jobs
    my $file_count      = $self->param('file_count');
    my $fasta_files_dir = $self->param('fasta_files_dir');
    my $bam_files_dir   = $self->param('bam_files_dir');
    my @jobs;
    my $i = 1;
    while ($i <= $file_count) {
        push @jobs, {
            'file_number'   => $i,
            'bam_files_dir' => $bam_files_dir,
            'fasta_file'    => "$fasta_files_dir/$i.fa",
            'sam_file'      => "$bam_files_dir/$i.sam",
            'bam_file'      => "$bam_files_dir/$i.bam",
            'err_file'      => "$bam_files_dir/$i.err",
            'out_file'      => "$bam_files_dir/$i.out",
        };
        $i++;
    }
    $self->dataflow_output_id(\@jobs, 2);
}

sub generate_mapping_input {
    my $self = shift;

    my $old_assembly_fasta_file_dir = $self->param('old_assembly_fasta_file_dir');
    my $fasta_db = Bio::DB::Fasta->new($old_assembly_fasta_file_dir, -reindex => 1);
    $self->param('fasta_db', $fasta_db);

    # store end-coordinates for all seq_regions to check that variation_location + flank_seq_length < slice_end
    my $cdba = $self->param('cdba');
    my $sa = $cdba->get_SliceAdaptor;
    # don't include asm exceptions
    my $slices = $sa->fetch_all('toplevel', undef, 0);
    my $seq_regions = {};
    foreach my $slice (@$slices) {
        my $end             = $slice->end;
        my $seq_region_name = $slice->seq_region_name;
        if ($end > $seq_regions->{$seq_region_name}->{end}) {
            $seq_regions->{$seq_region_name}->{end} = $end;
        }
    }
    foreach my $seq_region_name (keys %$seq_regions) {
        my $end = $seq_regions->{$seq_region_name}->{end};
    }

    $self->param('seq_regions', $seq_regions);

    my $dump_features_dir = $self->param('dump_features_dir');
    my $fasta_files_dir   = $self->param('fasta_files_dir');
    my $pipeline_dir      = $self->param('pipeline_dir');
    my $fh_allele_length = FileHandle->new("$pipeline_dir/qc_allele_string_length.txt", 'w');
    my $fh_ref_seq       = FileHandle->new("$pipeline_dir/qc_ref_seq.txt", 'w');

    $self->param('qc_allele_string_length', $fh_allele_length);
    $self->param('qc_ref_seq', $fh_ref_seq);

    my $variants_with_multi_map = {};

    opendir(DIR, $dump_features_dir) or die $!;
    while (my $file = readdir(DIR)) {
        if ($file =~ /^(.+)\.txt$/) {
            my $file_number = $1;
            my $fh = FileHandle->new("$dump_features_dir/$file", 'r');
            my $fh_fasta_file = FileHandle->new("$fasta_files_dir/$file_number.fa", 'w');
            while (<$fh>) {
                chomp;
                my $data = $self->read_line($_);
                my $seq_region_name = $data->{seq_region_name},
                my $feature_id      = $data->{variation_feature_id};
                my $start           = $data->{seq_region_start};
                my $end             = $data->{seq_region_end};
                my $strand          = $data->{seq_region_strand};
                my $allele_string   = $data->{allele_string};
                my $map_weight      = $data->{map_weight};
                my $variation_name  = $data->{variation_name};
                if ($map_weight > 1) {
                    $variants_with_multi_map->{$variation_name}++;
                    next;
                }
                my ($flank_start, $upstream_flank_length, $downstream_flank_length, $flank_end, $variant_length) = @{$self->flank_coordinates($seq_region_name, $start, $end, $strand)};

                # variant surrounded by flank sequences
                my $query_sequence = $self->get_query_sequence($seq_region_name, $flank_start, $flank_end, $strand);

                # replace empty space with underscores <_>
                $allele_string =~ s/\s/_/g;
                $self->qc_alleles($query_sequence, $upstream_flank_length, $variant_length, $allele_string, $variation_name);
#                $allele_string = $self->qc_alleles($config, $query_sequence, $upstream_flank_length, $variant_length, $allele_string, $variation_name);

                my $name = "$seq_region_name:$start:$end:$strand:$variation_name";
                my $id = ">$feature_id-$upstream_flank_length-$variant_length-$downstream_flank_length-$name";
                print $fh_fasta_file "$id\n$query_sequence\n";
            } # end while (read feature file for seq_region)
            $fh->close();
            $fh_fasta_file->close();
        }
    }

    foreach my $file_type (qw/qc_allele_string_length qc_ref_seq/) {
        my $fh = $self->param($file_type);
        $fh->close();
    }

    # store multi map
    my $fh_qc_multi_map = FileHandle->new("$pipeline_dir/multi_map.txt", 'w');

    foreach my $variation_name (keys %$variants_with_multi_map) {
        my $count = $variants_with_multi_map->{$variation_name};
        print $fh_qc_multi_map "$variation_name\t$count\n";
    }
    $fh_qc_multi_map->close();
}
{'fraction_non_alphanumeric': 0.14554692169754932, 'fraction_numerical': 0.002839210998206814, 'mean_word_length': 2.818830242510699, 'pattern_counts': {'":': 0, '<': 15, '<?xml version=': 0, '>': 155, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, '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': '4554529', 'n_tokens_mistral': 4677, 'n_tokens_neox': 4429, 'n_words': 894}