Commit 6053497d authored by Andy Yates's avatar Andy Yates
Browse files

FASTA dumper pipeline transfer from personal workspace. Consult the personal...

FASTA dumper pipeline transfer from personal workspace. Consult the personal kt7 workspace for more detailed development notes.
parent d2377aa2
package Bio::EnsEMBL::Pipeline::FASTA::Base;
use strict;
use warnings;
use base qw/Bio::EnsEMBL::Hive::Process/;
use Bio::EnsEMBL::Utils::Exception qw/throw/;
use File::Find;
use File::Spec;
use File::Path qw/mkpath/;
use POSIX qw/strftime/;
sub get_Slices {
my ($self, $type) = @_;
my $sa = $self->get_DBAdaptor($type)->get_SliceAdaptor();
return [ sort { $a->length() <=> $b->length() } @{$sa->fetch_all('toplevel', undef, 1, undef, undef)} ];
}
# Registry is loaded by Hive (see beekeeper_extra_cmdline_options() in conf)
sub get_DBAdaptor {
my ($self, $type) = @_;
my $species = $self->param('species');
$type ||= 'core';
return Bio::EnsEMBL::Registry->get_DBAdaptor($species, $type);
}
sub cleanup_DBAdaptor {
my ($self, $type) = @_;
my $dba = $self->get_DBAdaptor($type);
$dba->clear_caches;
$dba->dbc->disconnect_if_idle;
return;
}
sub get_dir {
my ($self, @extras) = @_;
my $base_dir = $self->param('base_path');
my $dir = File::Spec->catdir($base_dir, @extras);
mkpath($dir);
return $dir;
}
sub fasta_path {
my ( $self, @extras ) = @_;
return $self->get_dir('fasta', $self->param('species'), @extras);
}
sub web_name {
my ($self) = @_;
# my $mc = $self->get_DBAdaptor()->get_MetaContainer();
# my $name = $mc->single_value_by_key('species.url'); # change back
my $name = ucfirst($self->production_name());
return $name;
}
sub assembly {
my ($self) = @_;
my $dba = $self->get_DBAdaptor();
return $dba->get_CoordSystemAdaptor()->fetch_all()->[0]->version();
}
sub production_name {
my ($self, $name) = @_;
my $dba;
if($name) {
$dba = Bio::EnsEMBL::Registry->get_DBAdaptor($name, 'core');
}
else {
$dba = $self->get_DBAdaptor();
}
my $mc = $dba->get_MetaContainer();
my $prod = $mc->get_production_name();
$dba->dbc()->disconnect_if_idle();
return $prod;
}
sub old_path {
my ($self, $species) = @_;
my $base = $self->param('ftp_dir');
my $prod = $self->production_name($species);
my $version = $self->param('previous_version');
my $dir = File::Spec->catdir($base, "release-$version", 'fasta', $prod, 'dna');
}
# Closes file handle, and deletes the file stub if it contains no data
# Returns success type
sub tidy_file_handle {
my ($self, $fh, $path) = @_;
if ($fh->tell() == 0) {
$fh->close;
unlink($path);
return 1;
}
close($fh);
return 0;
}
sub info {
my ($self, $msg, @params) = @_;
if ($self->debug() > 1) {
my $formatted_msg;
if(scalar(@params)) {
$formatted_msg = sprintf($msg, @params);
}
else {
$formatted_msg = $msg;
}
printf STDERR "INFO: %s %s\n", strftime('%c',localtime()), $formatted_msg;
}
return
}
sub fine {
my ($self, $msg, @params) = @_;
if ($self->debug() > 2) {
my $formatted_msg;
if(scalar(@params)) {
$formatted_msg = sprintf($msg, @params);
}
else {
$formatted_msg = $msg;
}
printf STDERR "FINE: %s %s\n", strftime('%c',localtime()), $formatted_msg;
}
return
}
sub find_files {
my ($self, $dir, $boolean_callback) = @_;
$self->throw("Cannot find path $dir") unless -d $dir;
my @files;
find(sub {
my $path = $File::Find::name;
if($boolean_callback->($_)) {
push(@files, $path);
}
}, $dir);
return \@files;
}
1;
=pod
=head1 LICENSE
Copyright (c) 1999-2012 The European Bioinformatics Institute and
Genome Research Limited. All rights reserved.
This software is distributed under a modified Apache license.
For license details, please see
http://www.ensembl.org/info/about/code_licence.html
=head1 CONTACT
Please email comments or questions to the public Ensembl
developers list at <dev@ensembl.org>.
Questions may also be sent to the Ensembl help desk at
<helpdesk@ensembl.org>.
=head1 NAME
Bio::EnsEMBL::Pipeline::FASTA::BlastIndexer
=head1 DESCRIPTION
Creates 2bit file of the given GZipped file. The resulting index
is created under the parameter location I<base_path> in blat. The filename
is prefixed with the port number of the blat server this file should be
run on.
The module also performs filtering of non-reference sequence regions
and can filter the redundant Y chromosome piece for human (as 2bit does
not like repeated sequence region names).
Allowed parameters are:
=over 8
=item file - The file to index
=item program - The location of the faToTwoBit program
=item port_offset - Value to add onto the species_id from the website DB
to name the file correctly
=item base_path - The base of the dumps
=item index - The type of file to index; supported values are empty,
I<dna> or I<dna_rm>. If specified we will look for this
string in the filename surrounded by '.' e.g. .dna.
=back
The registry should also have a DBAdaptor for the website schema
registered under the species B<multi> and the group B<web> for species_id to
Blat port number.
=cut
package Bio::EnsEMBL::Pipeline::FASTA::BlatIndexer;
use strict;
use warnings;
use base qw/Bio::EnsEMBL::Pipeline::FASTA::Indexer/;
use File::Spec;
use File::stat;
use Bio::EnsEMBL::Utils::IO qw/work_with_file/;
use Bio::EnsEMBL::Utils::Exception qw/throw/;
use Bio::EnsEMBL::Registry;
sub param_defaults {
my ($self) = @_;
return {
program => 'faToTwoBit',
port_offset => 30000,
'index' => 'dna', #or dna_rm/''
};
}
sub run {
my ($self) = @_;
if($self->run_indexing()) {
$self->SUPER::run();
}
return;
}
sub run_indexing {
my ($self) = @_;
my $index = $self->param('index');
if($index) {
my $file = $self->param('file');
return (index($file, ".${index}.") > -1) ? 1 : 0;
}
return 1;
}
sub index_file {
my ($self, $file) = @_;
my $target_file = $self->target_file();
my $cmd = sprintf(q{%s %s %s},
$self->param('program'), $file, $target_file);
$self->info('About to run "%s"', $cmd);
system($cmd) and throw "Cannot run program '$cmd'";
unlink $file or throw "Cannot remove the file '$file' from the filesystem: $!";
#Check the file size. If it's 16 bytes then reject as that is an empty file for 2bit
my $filesize = stat($target_file)->size();
if($filesize <= 16) {
unlink $file;
my $msg = sprintf(
'The file %s produced a 2bit file %d byte(s). Lower than 17 bytes therefore empty 2 bit file',
$file, $filesize
);
$self->throw($msg);
}
return;
}
sub decompress {
my ($self) = @_;
#If we have no non-reference seq regions then use normal decompress
if(! $self->has_non_refs()) {
return $self->SUPER::decompress();
}
#Filter for non-refs
my $source = $self->param('file');
my $target_dir = $self->target_dir();
my ($vol, $dir, $file) = File::Spec->splitpath($source);
$file =~ s/.gz$//;
my $target = File::Spec->catdir($target_dir, $file);
my $allowed_regions = $self->allowed_regions();
$self->info('Decompressing %s to %s', $source, $target);
open my $src_fh, '-|', "zcat $source" or throw "Cannot decompress $source to $target";
work_with_file($target, 'w', sub {
my ($trg_fh) = @_;
my $transfer = 0;
while(my $line = <$src_fh>) {
#HEADER
if(index($line, '>') == 0) {
#regex is looking for NNN:NNN:NNN:1:80:1 i.e. the name
my ($name) = $line =~ />.+\s(.+:.+:.+:\d+:\d+:\d+)/;
$transfer = ($allowed_regions->{$name}) ? 1 : 0;
if($transfer) {
$self->info('%s was an allowed Slice', $name);
}
else {
$self->info('%s will be skipped; not a reference Slice', $name);
}
}
print $trg_fh $line if $transfer;
}
});
close($src_fh);
return $target;
}
sub allowed_regions {
my ($self) = @_;
my $prod_name = $self->production_name();
my @slices = grep { $_->is_reference() } @{$self->get_Slices()};
my %hash;
foreach my $slice (@slices) {
if($prod_name eq 'homo_sapiens') {
#Skip the extra unused human region
if($slice->seq_region_name() eq 'Y' && $slice->end() < 2649521 ) {
next;
}
}
$hash{$slice->name()} = 1;
}
return \%hash;
}
#Filename like 30061.Homo_sapiens.GRCh37.2bit
sub target_filename {
my ($self) = @_;
my $port = $self->blat_port();
my $name = $self->web_name();
my $assembly = $self->assembly();
return join(q{.}, $port, $name, $assembly, '2bit');
}
sub target_file {
my ($self) = @_;
my $target_dir = $self->target_dir();
my $target_filename = $self->target_filename();
return File::Spec->catfile($target_dir, $target_filename);
return;
}
sub target_dir {
my ($self) = @_;
return $self->get_dir('blat');
}
sub blat_port {
my ($self) = @_;
my $dba = Bio::EnsEMBL::Registry->get_DBAdaptor('multi', 'web');
my $id = $dba->dbc()->sql_helper()->execute_single_result(
-SQL => 'select species_id from species where name =?',
-PARAMS => [$self->web_name()]
);
return $id + $self->param('port_offset');
}
sub has_non_refs {
my ($self) = @_;
my $sql = <<'SQL';
select count(*)
from attrib_type at
join seq_region_attrib sra using (attrib_type_id)
join seq_region sr using (seq_region_id)
join coord_system cs using (coord_system_id)
where cs.species_id =?
and at.code =?
SQL
my $dba = $self->get_DBAdaptor();
return $dba->dbc()->sql_helper()->execute_single_result(
-SQL => $sql, -PARAMS => [$dba->species_id(), 'non_ref']);
}
1;
=pod
=head1 LICENSE
Copyright (c) 1999-2012 The European Bioinformatics Institute and
Genome Research Limited. All rights reserved.
This software is distributed under a modified Apache license.
For license details, please see
http://www.ensembl.org/info/about/code_licence.html
=head1 CONTACT
Please email comments or questions to the public Ensembl
developers list at <dev@ensembl.org>.
Questions may also be sent to the Ensembl help desk at
<helpdesk@ensembl.org>.
=head1 NAME
Bio::EnsEMBL::Pipeline::FASTA::ChecksumGenerator
=head1 DESCRIPTION
Creates a CHECKSUMS file in the given directory which is produced from running
the sum command over every file in the directory. This excludes the CHECKSUMS
file, parent directory or any hidden files.
Allowed parameters are:
=over 8
=item dir - The directory to generate checksums for
=item gzip - If the resulting file should be gzipped. Defaults to false
=back
=cut
package Bio::EnsEMBL::Pipeline::FASTA::ChecksumGenerator;
use strict;
use warnings;
use base qw/Bio::EnsEMBL::Pipeline::FASTA::Base/;
use File::Spec;
use Bio::EnsEMBL::Utils::IO qw/work_with_file gz_work_with_file/;
sub param_defaults {
my ($self) = @_;
return {
gzip => 0
};
}
sub fetch_input {
my ($self) = @_;
my $dir = $self->param('dir');
$self->throw("No 'dir' parameter specified") unless $dir;
$self->throw("Dir $dir does not exist") unless -d $dir;
return;
}
sub run {
my ($self) = @_;
my @checksums;
my $dir = $self->param('dir');
$self->info('Checksumming directory %s', $dir);
opendir(my $dh, $dir) or die "Cannot open directory $dir";
my @files = sort { $a cmp $b } readdir($dh);
closedir($dh) or die "Cannot close directory $dir";
foreach my $file (@files) {
next if $file =~ /^\./; #hidden file or up/current dir
next if $file =~ /^CHECKSUM/;
my $path = File::Spec->catfile($dir, $file);
my $checksum = $self->checksum($path);
push(@checksums, [$checksum, $file])
}
$self->param('checksums', \@checksums);
return;
}
sub write_output {
my ($self) = @_;
my $dir = $self->param('dir');
my $checksum = File::Spec->catfile($dir, 'CHECKSUMS');
$checksum .= '.gz' if $self->param('gzip');
if(-f $checksum) {
$self->info('Checksum file already exists. Removing');
unlink $checksum;
}
my @checksums = @{$self->param('checksums')};
return unless @checksums;
my $writer = sub {
my ($fh) = @_;
foreach my $entry (@checksums) {
my $line = join(qq{\t}, @{$entry});
print $fh $line;
print $fh "\n";
}
return;
};
my @params = ($checksum, 'w', $writer);
if($self->param('gzip')) {
gz_work_with_file(@params);
}
else {
work_with_file(@params);
}
$self->permissions($checksum);
return;
}
sub checksum {
my ($self, $path) = @_;
my $checksum = `sum $path`;
$checksum =~ s/\s* $path//xms;
chomp($checksum);
return $checksum;
}
sub permissions {
my ($self, $file) = @_;
my $mode = 0666;
chmod($mode, $file) or $self->throw("Cannot perform the chmod to mode $mode for file $file");
return;
}
1;
=pod
=head1 LICENSE
Copyright (c) 1999-2012 The European Bioinformatics Institute and
Genome Research Limited. All rights reserved.
This software is distributed under a modified Apache license.
For license details, please see
http://www.ensembl.org/info/about/code_licence.html
=head1 CONTACT
Please email comments or questions to the public Ensembl
developers list at <dev@ensembl.org>.
Questions may also be sent to the Ensembl help desk at
<helpdesk@ensembl.org>.
=head1 NAME
Bio::EnsEMBL::Pipeline::FASTA::ConcatFiles
=head1 DESCRIPTION
Performs a find in the DNA dumps directory for the given species and then
concats files which match a specified name pattern. We only allow
two types of concats; DNA and RM DNA. The concat file is a series
of cat command calls from all other Gzipped FASTA dumps (allowed under
the GZip specification).
Allowed parameters are:
=over 8
=item version - Needed to build the target path
=item species - Required to indicate which species we are working with
=item data_type - The type of data to work with. Can be I<dna> or T<dna_rm>
=item base_path - The base of the dumps
=back
=cut
package Bio::EnsEMBL::Pipeline::FASTA::ConcatFiles;
use strict;
use warnings;
use base qw/Bio::EnsEMBL::Pipeline::FASTA::Base/;
use Bio::EnsEMBL::Utils::Exception qw/throw/;
use File::Spec;
sub param_defaults {
my ($self) = @_;
return {
dna => {
regex => qr/.+\.dna\..+\.fa\.gz$/,
},
dna_rm => {
regex => qr/.+\.dna_rm\..+\.fa\.gz$/,
},
};
}
sub fetch_input {
my ($self) = @_;
foreach my $key (qw/data_type species version base_path/) {
$self->throw("Cannot find the required parameter $key") unless $self->param($key);
}
return;
}