Bio-ViennaNGS

 view release on metacpan or  search on metacpan

lib/Bio/ViennaNGS/UCSC.pm  view on Meta::CPAN

# -*-CPerl-*-
# Last changed Time-stamp: <2017-06-26 17:28:26 mtw>

package Bio::ViennaNGS::UCSC;

use Bio::ViennaNGS;
use Exporter;
use strict;
use warnings;
use Template;
use Cwd;
use File::Basename;
use File::Find;
use IPC::Cmd qw(can_run run);
use File::Share ':all';
use Path::Class;
use Data::Dumper;
use Carp;
use Bio::ViennaNGS::Util qw(fetch_chrom_sizes bed2bigBed);
use version; our $VERSION = version->declare("$Bio::ViennaNGS::VERSION");

our @ISA = qw(Exporter);

our @EXPORT_OK = qw( make_assembly_hub make_track_hub retrieve_color
make_track make_multi_bigwig_container_track
make_bigwig_container_track retrieve_bigwig_tracks
retrieve_bigbed_url_tracks parse_fasta_header valid_ncbi_accession
modify_fasta_header make_group retrieve_chromosome_size
write_chromosome_size_file convert_tracks);

our @EXPORT = ();

#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#
#^^^^^^^^^^^ Subroutines ^^^^^^^^^^#
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#

sub make_assembly_hub{
  my ($fasta_path, $filesdir, $basedir, $baseURL, $big_wig_ids, $log) = @_;
  my ($basename,$dir,$ext);
  my $this_function = (caller(0))[3];
  #check arguments
  croak ("ERROR [$this_function] \$fasta_path does not exist\n") 
    unless (-e $fasta_path);
  croak ("ERROR [$this_function] \$basedir does not exist\n") 
    unless (-d $basedir);
  croak ("ERROR [$this_function]: no URL (network location for upload to UCSC) provided") 
    unless(defined $baseURL);

  unless ($baseURL =~ /\/$/) { $baseURL .= "/"; }

  my $tmp_path = dist_file('Bio-ViennaNGS', "hub.txt" );
  ($basename,$dir,$ext) = fileparse($tmp_path,qr/\..*/);
  my $template_path = dir($dir,"template-AssemblyHub");

  croak ("ERROR [$this_function] template directory not found\n") 
    unless (-d $template_path);
  my $faToTwoBit = can_run('faToTwoBit') or
    croak ("ERROR [$this_function] faToTwoBit is not installed!");

  # bedfiles path
  my $parsedHeaderRef = parse_fasta_header($fasta_path,$this_function);
  my @parsedHeader = @$parsedHeaderRef;
  my $unchecked_accession = $parsedHeader[0];
  my $scientificName = $parsedHeader[1];
  my $accession = valid_ncbi_accession($unchecked_accession);
  # create assembly hub directory structure
  my $assembly_hub_name = "assemblyHub";
  my $assembly_hub_directory = dir($basedir, $assembly_hub_name);
  my $genome_assembly_name = $accession;
  my $genome_assembly_directory = dir($assembly_hub_directory,$genome_assembly_name);
  mkdir $assembly_hub_directory;
  mkdir $genome_assembly_directory;

  if (defined $log){
    open(LOG, ">>", $log) or croak "$!";
    print LOG "LOG Base directory:          $assembly_hub_directory\n";
    print LOG "LOG Assembly Hub directory:  $genome_assembly_directory\n";
    close(LOG);
  }

  #2-bit fasta file conversion
  my $fa_modified = file($assembly_hub_directory, $accession.".fa");
  modify_fasta_header($fasta_path,$fa_modified,$accession);
  my $twoBit = file($genome_assembly_directory, $accession.".2bit");
  my $fastaToTwobit_cmd = "$faToTwoBit $fa_modified $twoBit";

  if (defined $log){
    open(LOG, ">>", $log) or croak "$!";
    print LOG "LOG [$this_function] $fastaToTwobit_cmd\n";
    close(LOG);
  }

  my( $success, $error_message, $full_buf, $stdout_buf, $stderr_buf ) =
    run( command => $fastaToTwobit_cmd, verbose => 0 );
  if( !$success ) {
    print STDERR "ERROR [$this_function] External command call unsuccessful\n";
    print STDERR "ERROR: this is what the command printed:\n";
    print join "", @$full_buf;
    croak $!;
  }

  #template definition
  my $template = Template->new({
                                INCLUDE_PATH => ["$template_path"],
                                RELATIVE=>1,
  });

lib/Bio/ViennaNGS/UCSC.pm  view on Meta::CPAN


  #construct group.txt
  my $group_txt_path = file($genome_assembly_directory, "groups.txt")->stringify;
  my $group_txt_file = 'groups.txt';
  my $group_txt_vars = 
    {
     groups  => "$groups",
    };
  $template->process($group_txt_file,$group_txt_vars,$group_txt_path) or
    croak "Template process failed: ", $template->error(), "\n";

  #construct big bed
  my $tracksList;
  #Bigbeds are only created from infolder
  unless($filesdir =~ /-/){
    my $chromosome_size = retrieve_chromosome_size($fasta_path);
    my $chromosome_size_filepath = file($genome_assembly_directory,"$accession.chrom.sizes");
    write_chromosome_size_file($chromosome_size_filepath,$accession,$chromosome_size);
    convert_tracks($filesdir, $genome_assembly_directory, $chromosome_size_filepath, $log);
    my @trackfiles = retrieve_tracks($genome_assembly_directory, $baseURL, $assembly_hub_name, $accession);

    foreach my $track (@trackfiles){
      my $trackString = make_track(@$track);
      $tracksList .= $trackString;
    }
  }

  #big wigs
  my $bigwig_tracks_string = "";
  unless($big_wig_ids=~/^-$/){
    $bigwig_tracks_string = retrieve_bigwig_tracks($genome_assembly_directory, $baseURL, $assembly_hub_name, $accession, $big_wig_ids);
  }
  $tracksList .= $bigwig_tracks_string;
  #construct trackDb.txt
  my $trackDb_txt_path = file($genome_assembly_directory, "trackDb.txt")->stringify;
  my $trackDb_txt_file = 'trackDb.txt';
  my $trackDb_txt_vars =
    {
     tracks  => "$tracksList"
    };
  $template->process($trackDb_txt_file,$trackDb_txt_vars,$trackDb_txt_path) or
    croak "Template process failed: ", $template->error(), "\n";

  if (defined $log){
    open(LOG, ">>", $log) or croak "$!";
    print LOG "LOG Assembly Hub created\n";
    close(LOG);
  }
}

sub make_track_hub{
  my ($species, $basedir, $baseURL, $big_bed_ids, $big_wig_ids, $log) = @_;
  my ($basename,$dir,$ext);
  my $this_function = (caller(0))[3];

  #check arguments
  croak ("ERROR [$this_function] \no species provided\n")
    unless ($species);
  croak ("ERROR [$this_function] \$basedir does not exist\n")
    unless (-d $basedir);
  croak ("ERROR [$this_function]: no URL (network location for upload to UCSC) provided")
    unless(defined $baseURL);

  if (defined $log){
    open(LOG, ">>", $log) or croak "$!";
  }

  unless ($baseURL =~ /\/$/) { $baseURL .= "/"; }

  my $tmp_path = dist_file('Bio-ViennaNGS', "hub.txt" );
  ($basename,$dir,$ext) = fileparse($tmp_path,qr/\..*/);
  my $template_path = dir($dir,"template-TrackHub");

  croak ("ERROR [$this_function] template directory not found\n") 
    unless (-d $template_path);
  my $faToTwoBit = can_run('faToTwoBit') or
    croak ("ERROR [$this_function] faToTwoBit is not installed!");

  # create track hub directory structure
  my $track_hub_name = "trackHub";
  my $track_hub_directory = dir($basedir, $track_hub_name);
  my $genome_assembly_name = $species;
  my $genome_assembly_directory = dir($track_hub_directory,$genome_assembly_name);
  mkdir $track_hub_directory;
  mkdir $genome_assembly_directory;
  if (defined $log){
    print LOG "LOG Base directory:          $track_hub_directory\n";
    print LOG "LOG Track Hub directory:  $genome_assembly_directory\n";
  }

  #template definition
  my $template = Template->new({
                                INCLUDE_PATH => ["$template_path"],
                                RELATIVE=>1,
  });

  #construct hub.txt
  my $hubtxt_path = file($track_hub_directory,"hub.txt")->stringify;
  my $hubtxt_file = "hub.txt";
  my $hubtxt_vars =
    {
     hubName => $species,
     shortLabel => $species,
     longLabel => $species,
     genomesFile => "genome.txt",
     email => 'email',
     descriptionURL => "$baseURL" . "description.html"
    };
  $template->process($hubtxt_file,$hubtxt_vars,"$hubtxt_path") ||
    croak "Template process failed: ", $template->error(), "\n";

  #construct genome.txt
  my $genometxt_path = file($track_hub_directory, "genome.txt")->stringify;
  my $genometxt_file = "genome.txt";
  my $genometxt_vars =
    {
     genome => $species,
     trackDb => file($species, "trackDb.txt"),
     groups => file($species, "groups.txt"),
     description => "$species",
     twoBitPath => file($species,$species.".2bit"),

lib/Bio/ViennaNGS/UCSC.pm  view on Meta::CPAN

  return 1;
}

1;
__END__

=head1 NAME

Bio::ViennaNGS::UCSC - Perl extension for easy UCSC Genome Browser
integration.

=head1 SYNOPSIS

  use Bio::ViennaNGS::UCSC;

=head1 DESCRIPTION

L<Bio::ViennaNGS::UCSC> facilitates routine tasks for automatic
visualization of next-generation sequencing data in the UCSC Genome
Browser.

=head2 EXPORT

Routines: 
  make_assembly_hub, make_track_hub

Variables:
  none

=head3 make_assembly_hub()

Build an UCSC Genome Browser Assembly Hub. This function takes 4
parameters:

=over

=item 1 absolute path of input fasta file(e.g. /home/user/input.fa)

=item 2 path to the ouput directory (e.g. /home/user/assemblyhubs/)

=item 3 base URL under which the new Assembly Hub will be available on
the Internet, e.g. http://www.foo.com/folder/. The generated Assembly
Hub files must be copied to the corresponding file system location on
the Web Server and the URL readable for the UCSC Genome Browser server.

=item 4 path for the log file (/home/user/logs/assemblyhubconstructionlog)

=back

=head3 make_track_hub()

Build track hubs for the UCSC genome browser.
This function takes 4 parameters:

=over

=item 1 chromosome id as used in existing ucsc assembly hub (e.g. chr1)

=item 2 path to the ouput directory (e.g. /home/user/assemblyhubs/)

=item 3 base URL where the output folder will be placed for upload to the UCSC genome browser (e.g. http://www.foo.com/folder/)

=item 4 path for the log file (/home/user/logs/assemblyhubconstructionlog)

=back

=head1 SEE ALSO

=over

=item L<Bio::ViennaNGS>

=back

=head1 AUTHORS

=over

=item Michael T. Wolfinger, E<lt>michael@wolfinger.euE<gt>

=item Florian Eggenhofer, E<lt>florian.eggenhofer@univie.ac.atE<gt>

=back

=head1 COPYRIGHT AND LICENSE

Copyright (C) 2014-2017 Michael T. Wolfinger, E<lt>michael@wolfinger.euE<gt>

This library is free software; you can redistribute it and/or modify
it under the same terms as Perl itself, either Perl version 5.10.0 or,
at your option, any later version of Perl 5 you may have available.

This program is distributed in the hope that it will be useful,
but WITHOUT ANY WARRANTY; without even the implied warranty of
MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.


=cut



( run in 3.267 seconds using v1.01-cache-2.11-cpan-b16cb0d3907 )