<commit_before>#!/usr/bin/env perl

# Copyright [1999-2015] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
# Copyright [2016] 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 NAME

  prepare_cdnas.pl

=head1 SYNOPSIS

  prepare_cdnas.pl

=head1 DESCRIPTION

  prepare_cdnas.pl makes a file of fasta sequences and also prints out
  an annotation file for use by exonerate's cdna2genome

  The script reads in each cDNA, checks the cDNA against
  the Ensembl kill_list DB if required, and clips it using
  Bio::EnsEMBL::Analysis::Tools::PolyAClipping. It then looks for the
  cDNA in Mole and finds its coordinates.

  If the coordinates are not found or can't be parsed, a message is
  sent to STDERR and the cDNA is not written to the outfile. (Parsing
  problem often caused by the cDNA not having exact start/end position
  data in the EMBL-CDS DB.) If the cDNA is found, we check that it has
  not been clipped into the CDS and then write the cDNA's fasta to
  outfile (/path/to/cdna.clip) and its annotation to annotation_file
  (/path/to/cdna.annot).

  The annotation file contains additional information for the user
  to know what has happened during clipping. In order for this
  annotation file to be used by exonerate's cdna2genome model, the
  annotation file will have to be edited to contain only the first
  4 columns of each line eg.

    less /path/to/cdna.annot | awk '{print $1"\t"$2"\t"$3"\t"$4}'

=head1 OPTIONS

  -dbnames      database names for Mole databases separated by commas without
                any whitespaces, eg. embl_91,refseq_23

  -infile       path to the original downloaded cDNA fasta file

  -outfile      path to the clipped cDNA fasta file to be written

  -min_length   minimum length of cDNA after its polyA tail has been clipped.
                set to be 0 by default.

  -annotation   path to the annotation file

=head1 EXAMPLE

  perl prepare_cdnas.pl -infile all_DL_cDNAs.fa -outfile clipped_annotated_cDNAs.fa\
  -annotation cDNAs.annotate -min_length 60 -dbnames embl_91,refseq_23


=cut

use strict;
use warnings;

use Bio::EnsEMBL::Analysis::Tools::PolyAClipping;
use Bio::SeqIO;
use Bio::EnsEMBL::KillList::KillList;
use Getopt::Long qw(:config no_ignore_case);
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
use Bio::EnsEMBL::ExternalData::Mole::DBXref;
use Bio::EnsEMBL::ExternalData::Mole::Entry;
use Bio::EnsEMBL::ExternalData::Mole::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Pipeline::SeqFetcher::Mfetch;

$| = 1; # disable buffering

my (
        @dbnames,
        $infile,
        $outfile,
        $annotation_file,
);

my $buffer      = 10;
my $window_size = 3;
my $dbuser      = 'genero';
my $dbhost      = 'cbi5d';
my $dbport      = 3306;
my $min_length  = 0;

GetOptions(
        'dbnames|db|D=s'              => \@dbnames,
        'dbuser|user|u=s'               => \$dbuser,
        'dbhost|host|h=s'               => \$dbhost,
        'dbport|port|P=s'               => \$dbport,
        'infile=s'               => \$infile,
        'outfile=s'              => \$outfile,
        'annotation=s'           => \$annotation_file,
        'buffer=s'               => \$buffer,
        'window=s'               => \$window_size,
        'min_length:s'           => \$min_length,
);

# check commandline
if ($annotation_file) {
  if (!defined($dbhost) || !defined($dbuser) || !scalar(@dbnames)) {
    throw ("Please set dbhost (-hdbost), dbport (-dbport), dbnames (-dbnames) and dbuser (-dbuser)");
  }
}

if ( !defined($infile) || !defined($outfile) ) {
  throw( "Please enter a file name to read (-infile), "
       . "a file name to write to (-outfile) "
       . "and a file name for annotation (-annotation). "
       . "The infile must contain a list of accessions "
       . "whose CDS coordinates you want to find" );
}

@dbnames = split(/,/,join(',',@dbnames));

# connect to databases
my @dbs;
foreach my $dbname (@dbnames) {
  my $db = Bio::EnsEMBL::ExternalData::Mole::DBSQL::DBAdaptor->new(
        '-dbname' => $dbname,
        '-host'   => $dbhost,
        '-user'   => $dbuser,
        '-port'   => $dbport,
  );
  push @dbs, $db;
}

# open the downloaded cdna file
my $seqin  = new Bio::SeqIO( -file   => "<$infile",
                             -format => "Fasta",
                           );

my $seqout = new Bio::SeqIO( -file   => ">$outfile",
                             -format => "Fasta"
                           );

# Various counter to keep track of cDNAs lost/discarded in the process:

my $total_cDNAs          = 0;
my $X_count              = 0;
my $killed_count         = 0;
my $parse_problem        = 0;
my $partial_cds          = 0;
my $pe_low               = 0;
my $not_in_mole          = 0;
my $too_short_after_clip = 0;
my $AT_only_count        = 0;
my $cdna_written         = 0;
my $zero_length          = 0;

my $kill_list_object = Bio::EnsEMBL::KillList::KillList->new( -TYPE => "cDNA" );
my %kill_list = %{ $kill_list_object->get_kill_list() };

if ($annotation_file) { open( ANNOTATION, ">>$annotation_file" ) or die "Cannot open $annotation_file\n"; }

# create hash containing all cdnas from the input file
my %cdnas;
while (my $raw_cdna = $seqin->next_seq) {
  $total_cDNAs++;
  if (!$raw_cdna->seq) {
    print STDERR $raw_cdna->display_id . " has 0 length seq, skipping it.\n";
    $zero_length++;<commit_msg>Updated code so it works with the latest ensembl-killlist<commit_after>#!/usr/bin/env perl

# Copyright [1999-2015] Wellcome Trust Sanger Institute and the EMBL-European Bioinformatics Institute
# Copyright [2016] 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 NAME

  prepare_cdnas.pl

=head1 SYNOPSIS

  prepare_cdnas.pl

=head1 DESCRIPTION

  prepare_cdnas.pl makes a file of fasta sequences and also prints out
  an annotation file for use by exonerate's cdna2genome

  The script reads in each cDNA, checks the cDNA against
  the Ensembl kill_list DB if required, and clips it using
  Bio::EnsEMBL::Analysis::Tools::PolyAClipping. It then looks for the
  cDNA in Mole and finds its coordinates.

  If the coordinates are not found or can't be parsed, a message is
  sent to STDERR and the cDNA is not written to the outfile. (Parsing
  problem often caused by the cDNA not having exact start/end position
  data in the EMBL-CDS DB.) If the cDNA is found, we check that it has
  not been clipped into the CDS and then write the cDNA's fasta to
  outfile (/path/to/cdna.clip) and its annotation to annotation_file
  (/path/to/cdna.annot).

  The annotation file contains additional information for the user
  to know what has happened during clipping. In order for this
  annotation file to be used by exonerate's cdna2genome model, the
  annotation file will have to be edited to contain only the first
  4 columns of each line eg.

    less /path/to/cdna.annot | awk '{print $1"\t"$2"\t"$3"\t"$4}'

=head1 OPTIONS

  -dbnames      database names for Mole databases separated by commas without
                any whitespaces, eg. embl_91,refseq_23

  -infile       path to the original downloaded cDNA fasta file

  -outfile      path to the clipped cDNA fasta file to be written

  -min_length   minimum length of cDNA after its polyA tail has been clipped.
                set to be 0 by default.

  -annotation   path to the annotation file

=head1 EXAMPLE

  perl prepare_cdnas.pl -infile all_DL_cDNAs.fa -outfile clipped_annotated_cDNAs.fa\
  -annotation cDNAs.annotate -min_length 60 -dbnames embl_91,refseq_23


=cut

use strict;
use warnings;

use Bio::EnsEMBL::Analysis::Tools::PolyAClipping;
use Bio::SeqIO;
use Bio::EnsEMBL::KillList::KillList;
use Getopt::Long qw(:config no_ignore_case);
use Bio::EnsEMBL::Utils::Exception qw(throw warning);
use Bio::EnsEMBL::ExternalData::Mole::DBXref;
use Bio::EnsEMBL::ExternalData::Mole::Entry;
use Bio::EnsEMBL::ExternalData::Mole::DBSQL::DBAdaptor;
use Bio::EnsEMBL::Pipeline::SeqFetcher::Mfetch;

$| = 1; # disable buffering

my (
        @dbnames,
        $infile,
        $outfile,
        $annotation_file,
);

my $buffer      = 10;
my $window_size = 3;
my $dbuser      = 'genero';
my $dbhost      = 'cbi5d';
my $dbport      = 3306;
my $killist_dbname = 'gb_kill_list';
my $killist_dbuser;
my $killist_dbhost;
my $killist_dbport = 3306;
my $min_length  = 0;
my $tax_id;

GetOptions(
        'dbnames|db|D=s'              => \@dbnames,
        'dbuser|user|u=s'               => \$dbuser,
        'dbhost|host|h=s'               => \$dbhost,
        'dbport|port|P=s'               => \$dbport,
        'killdbnames|kdb=s'              => \$killist_dbname,
        'killdbuser|kuser=s'               => \$killist_dbuser,
        'killdbhost|khost=s'               => \$killist_dbhost,
        'killdbport|kport=s'               => \$killist_dbport,
        'infile=s'               => \$infile,
        'outfile=s'              => \$outfile,
        'annotation=s'           => \$annotation_file,
        'buffer=s'               => \$buffer,
        'window=s'               => \$window_size,
        'min_length:s'           => \$min_length,
        'tax_id:s'           => \$tax_id,
);

# check commandline
if ($annotation_file) {
  if (!defined($dbhost) || !defined($dbuser) || !scalar(@dbnames)) {
    throw ("Please set dbhost (-hdbost), dbport (-dbport), dbnames (-dbnames) and dbuser (-dbuser)");
  }
}

if ( !defined($infile) || !defined($outfile) ) {
  throw( "Please enter a file name to read (-infile), "
       . "a file name to write to (-outfile) "
       . "and a file name for annotation (-annotation). "
       . "The infile must contain a list of accessions "
       . "whose CDS coordinates you want to find" );
}

@dbnames = split(/,/,join(',',@dbnames));

# connect to databases
my @dbs;
foreach my $dbname (@dbnames) {
  my $db = Bio::EnsEMBL::ExternalData::Mole::DBSQL::DBAdaptor->new(
        '-dbname' => $dbname,
        '-host'   => $dbhost,
        '-user'   => $dbuser,
        '-port'   => $dbport,
  );
  push @dbs, $db;
}

# open the downloaded cdna file
my $seqin  = new Bio::SeqIO( -file   => "<$infile",
                             -format => "Fasta",
                           );

my $seqout = new Bio::SeqIO( -file   => ">$outfile",
                             -format => "Fasta"
                           );

# Various counter to keep track of cDNAs lost/discarded in the process:

my $total_cDNAs          = 0;
my $X_count              = 0;
my $killed_count         = 0;
my $parse_problem        = 0;
my $partial_cds          = 0;
my $pe_low               = 0;
my $not_in_mole          = 0;
my $too_short_after_clip = 0;
my $AT_only_count        = 0;
my $cdna_written         = 0;
my $zero_length          = 0;

my $kill_list_object = Bio::EnsEMBL::KillList::KillList->new( -TYPE => "cDNA",
    -KILL_LIST_DB => {
        '-dbname' => $killist_dbname,
        '-host'   => $killist_dbhost,
        '-user'   => $killist_dbuser,
        '-port'   => $killist_dbport,
    });
if ($tax_id) {
  $kill_list_object->FILTER_PARAMS->{-for_species} = [$tax_id];
}
my %kill_list = %{ $kill_list_object->get_kill_list() };

if ($annotation_file) { open( ANNOTATION, ">>$annotation_file" ) or die "Cannot open $annotation_file\n"; }

# create hash containing all cdnas from the input file
my %cdnas;
while (my $raw_cdna = $seqin->next_seq) {
  $total_cDNAs++;
  if (!$raw_cdna->seq) {
    print STDERR $raw_cdna->display_id . " has 0 length seq, skipping it.\n";
    $zero_length++;
{'fraction_non_alphanumeric': 0.11578359399936164, 'fraction_numerical': 0.01005426109160549, 'mean_word_length': 2.7445473558410516, 'pattern_counts': {'":': 0, '<': 5, '<?xml version=': 0, '>': 71, 'https://': 0, 'lorem ipsum': 0, 'www.': 2, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'dir': 'git-commits-cleaned', 'id': '2251378', 'n_tokens_mistral': 4332, 'n_tokens_neox': 4043, 'n_words': 1298}