Bio-ViennaNGS
view release on metacpan or search on metacpan
lib/Bio/ViennaNGS/AnnoC.pm view on Meta::CPAN
lazy => 1,
);
before 'featstat' => sub {
my $self = shift;
$self->_get_nr_of_features();
};
sub _set_featstat {
my $self = shift;
my $this_function = (caller(0))[3];
my %fs = ();
confess "ERROR [$this_function] \$self->features not available"
unless ($self->has_features);
$fs{total} = 0;
$fs{origin} = "$this_function ".$VERSION;
$fs{count} = $self->nr_features;
foreach my $uid ( keys %{$self->features} ){
my $gbkey = ${$self->features}{$uid}->{gbkey};
$fs{total} += 1;
unless (exists $fs{$gbkey}){
$fs{$gbkey} = 0;
}
$fs{$gbkey} += 1;
}
return \%fs;
}
sub _get_nr_of_features {
my $self = shift;
my $this_function = (caller(0))[3];
confess "ERROR [$this_function] \$self->features not available"
unless ($self->has_features);
return (keys %{$self->features});
}
sub parse_gff {
my ($self,$in_file) = @_;
my ($i,$gffio,$header,$f,$gbkey);
my $this_function = (caller(0))[3];
$gffio = Bio::Tools::GFF->new(-file => $in_file,
-gff_version => 3,
);
$gffio->ignore_sequence(1);
if ($header = $gffio->next_segment() ){
$self->accession( $header->display_id() );
}
else{ carp "ERROR [$this_function] Cannot parse GFF header\n" }
lib/Bio/ViennaNGS/AnnoC.pm view on Meta::CPAN
}
}
$gffio->close();
}
sub features2bed {
my ($self,$gbkey,$dest,$bn,$log) = @_;
my ($chrom,$chrom_start,$chrom_end,$name,$score,$strand,$thick_start);
my ($thick_end,$reserved,$block_count,$block_sizes,$block_starts);
my @ft = ();
my $this_function = (caller(0))[3];
croak "ERROR [$this_function] $self->features not available"
unless ($self->has_features);
croak "ERROT [$this_function] $self->featstat not available"
unless ($self->has_featstat);
croak "ERROR [$this_function] $dest does not exist"
unless (-d $dest);
if (defined $log){open(LOG, ">>", $log) or croak $!;}
if (defined $gbkey){ # dump features of just one genbank key
lib/Bio/ViennaNGS/AnnoC.pm view on Meta::CPAN
sortbed($bedname_u,".",$bedname,1,undef); # sort bed file
} # end foreach
if (defined $log){close(LOG)};
}
sub feature_summary {
my ($self, $dest) = @_;
my ($fn,$fh);
my $this_function = (caller(0))[3];
croak "ERROR [$this_function] $dest does not exist\n"
unless (-d $dest);
croak "ERROR [$this_function] $self->accession not available\n"
unless ($self->has_accession);
$fn = dir($dest,$self->accession.".summary.txt");
open $fh, ">", $fn or croak $!;
print $fh "Accession\t ".$self->accession."\n";
lib/Bio/ViennaNGS/Bam.pm view on Meta::CPAN
sub split_bam {
my %data = ();
my @NHval = ();
my @processed_files = ();
my ($verbose,$nh_warning_issued) = (0)x2;
my ($bamfile,$reverse,$want_uniq,$want_bed,$dest_dir,$log) = @_;
my ($bam,$sam,$bn,$path,$ext,$header,$flag,$NH,$eff_strand,$tmp);
my ($bam_pos,$bam_neg,$tmp_bam_pos,$tmp_bam_neg,$bamname_pos,$bamname_neg);
my ($bed_pos,$bed_neg,$bedname_pos,$bedname_neg);
my ($seq_id,$start,$stop,$strand,$target_names,$id,$score);
my $this_function = (caller(0))[3];
my %count_entries = (
total => 0,
uniq => 0,
pos => 0,
neg => 0,
skip => 0,
cur => 0,
mult => 0,
se_alis => 0,
pe_alis => 0,
lib/Bio/ViennaNGS/Bam.pm view on Meta::CPAN
}
sub uniquify_bam {
my %data = ();
my ($bamfile,$dest,$log) = @_;
my ($bam, $bn,$path,$ext,$read,$header);
my ($tmp_uniq,$tmp_mult,$fn_uniq,$fn_mult,$bam_uniq,$bam_mult,$NH);
my ($count_all,$count_uniq,$count_mult,$unmapped,$nh_warning_issued) = (0)x5;
my @processed_files = ();
my @NHval = ();
my $this_function = (caller(0))[3];
my %count_entries = (
total => 0,
uniq => 0,
mult => 0,
unmapped => 0,
nonhtag => 0,
);
$data{count} = \%count_entries;
$data{nh_issues} = 0; # will be set to 1 if NH attribute is missing
lib/Bio/ViennaNGS/Bam.pm view on Meta::CPAN
my ($bamfile,$dest,$log) = @_;
my ($bam, $bn,$path,$ext,$read,$header,$ali);
my ($tmp_uniq,$tmp_mult,$fn_uniq,$fn_mult,$bam_uniq,$bam_mult,$NH);
my ($count_all,$count_uniq,$count_mult,$unmapped,$nh_warning_issued) = (0)x5;
my $pushme = 0;
my $allgomulti = 0;
my @processed_files = ();
my @NHval = ();
my $lastqname = undef;
my @band = ();
my $this_function = (caller(0))[3];
my %count_entries = (
total => 0,
uniq => 0,
mult => 0,
unmapped => 0,
nonhtag => 0,
);
$data{count} = \%count_entries;
$data{nh_issues} = 0; # will be set to 1 if NH attribute is missing
lib/Bio/ViennaNGS/Bed.pm view on Meta::CPAN
isa => 'Int',
builder => '_build_length',
lazy => 1,
predicate => 'has_length',
);
sub _build_length {
my ($self) = @_;
my ($i,$len) = (0)x2;
my $this_function = (caller(0))[3];
my @blockSizes = split (/,/ ,$self->blockSizes);
my $bc = scalar @blockSizes;
croak "ERROR [$this_function] invalid blockSount in BED12 line"
unless ($bc == $self->blockCount);
for ($i=0;$i<$bc;$i++){
$len += $blockSizes[$i];
}
$self->length($len);
}
sub as_bed_line {
my ($self,$n) = @_;
my $this_function = (caller(0))[3];
croak "ERROR [$this_function] no argument provided"
unless (defined $n);
croak "ERROR [$this_function] argument of as_bed_line() must be 6 or 12"
unless ( ($n == 6) | ($n == 12) );
my $bed6= join ("\t",
$self->chromosome,
$self->start,
$self->end,
$self->name,
$self->score,
lib/Bio/ViennaNGS/Expression.pm view on Meta::CPAN
has 'nr_features' => (
is => 'rw',
isa => 'Int',
predicate => 'has_features',
);
sub parse_readcounts_bed12 {
my ($self,$file) = @_;
my @mcData = ();
my ($i,$n) = (0)x2;
my $this_function = (caller(0))[3];
croak "ERROR [$this_function] readcount / multicov file $self->readcountfile not available\n"
unless (-e $file);
$self->readcountfile($file);
open (RC_IN, "< $file") or croak $!;
while (<RC_IN>){
$n++;
chomp;
# 0:chr|1:start|2:end|3:name|4:score|5:strand
lib/Bio/ViennaNGS/Expression.pm view on Meta::CPAN
# print Dumper(${$self->data}[$i]);
}
}
$self->nr_features($n);
close(RC_IN);
}
sub write_expression_bed12 {
my ($self,$item,$dest,$base_name) = @_;
my ($bedname,$bedname_u,$outfile,$feat,$i);
my $this_function = (caller(0))[3];
croak "ERROR [$this_function]: $dest does not exist\n"
unless (-d $dest);
$bedname = $base_name.".".$item.".multicov.bed12";
$bedname_u = $base_name.".".$item.".multicov.u.bed12";
$outfile = file($dest,$bedname_u);
croak "ERROR [$this_function]: $self->conds not available\n"
lib/Bio/ViennaNGS/Fasta.pm view on Meta::CPAN
has 'primaryseqH' => (
is => 'rw',
isa => 'HashRef',
predicate => 'has_primaryseq',
init_arg => undef,
);
sub BUILD {
my $self = shift;
my $this_function = (caller(0))[3];
$self->fastaids([$self->fasta->ids]);
confess "ERROR [$this_function] \$self->fastsids not available"
unless ($self->has_ids);
my %ps = ();
foreach my $id (@{$self->fastaids}){
$ps{$id} = $self->fasta->get_Seq_by_id($id);
}
$self->primaryseqH(\%ps);
}
lib/Bio/ViennaNGS/Fasta.pm view on Meta::CPAN
Returns : 1 if ID C<$id> was found, 0 else.
=back
=cut
sub stranded_subsequence {
my ($self,$id,$start,$end,$strand) = @_;
my ($this_function,$seq,$rc,$p,$obj);
$this_function = (caller(0))[3];
confess "ERROR [$this_function] start corrdinate must be <= end coordinate"
unless ($start <= $end);
confess "ERROR [$this_function] Id $id not found in input Fasta file"
unless ($self->has_sequid($id));
$seq = ${$self->primaryseqH}{$id}->subseq($start => $end);
if ($strand eq '-1' || $strand eq '-') {
$rc = revcom($seq);
$seq = $rc->seq();
}
#print "id:$id\nstart:$start\nend:$end\n";
lib/Bio/ViennaNGS/FeatureChain.pm view on Meta::CPAN
predicate => 'nr_entries',
init_arg => undef, # make this unsettable via constructor
builder => 'count_entries',
lazy => 1,
);
with 'Bio::ViennaNGS::FeatureBase';
sub BUILD {
my $self = shift;
my $this_function = (caller(0))[3];
confess "ERROR [$this_function] \$self->chain not available"
unless ($self->has_chain);
$self->count_entries();
}
sub count_entries {
my $self = shift;
my $cnt = scalar @{$self->chain};
$self->_entries($cnt);
}
before 'as_bed12_line' => sub {
my $self = shift;
my $this_function = (caller(0))[3];
$self->sip( sub { $_[0]->start <=> $_[1]->start} )
};
sub sort_chain_ascending {
my $self = shift;
my $this_function = (caller(0))[3];
$self->sip( sub { $_[0]->start <=> $_[1]->start} )
}
sub as_bed12_line{
my ($self,$name,$score,$strand) = @_;
my ($i,$chr,$start,$end,$feat,$bed12,$bsizes,$bstarts);
my $count=0;
my @blockSizes = ();
my @blockStarts = ();
# TODO check whether all features have the same chromosome id
lib/Bio/ViennaNGS/FeatureIO.pm view on Meta::CPAN
predicate => 'nr_entries',
init_arg => undef, # make this unsettable via constructor
builder => 'count_entries',
lazy => 1,
);
with 'Bio::ViennaNGS::FeatureBase';
sub BUILD { # call a parser method, depending on $self->instanceOf
my $self = shift;
my $this_function = (caller(0))[3];
my $type;
confess "ERROR [$this_function] \$self->file not available"
unless ($self->has_file);
confess "ERROR [$this_function] \$self->filetype not available"
unless ($self->has_filetype);
confess "ERROR [$this_function] \$self->instanceOf not available"
unless ($self->has_instanceOf);
if ($self->filetype eq "BedGraph"){
lib/Bio/ViennaNGS/FeatureIO.pm view on Meta::CPAN
sub count_entries {
my $self = shift;
my $cnt = scalar @{$self->data};
$self->_entries($cnt);
}
# TODO ensure FeatureBase is handled correctly
sub parse_bedgraph_file{
my ($self,$filename) = @_;
my $this_function = (caller(0))[3];
my ($file,$line,$entry,$chr,$start,$end,$val);
$file = read_file( $filename, array_ref => 1, chomp =>1 ) ;
foreach $line (@$file){
croak "ERROR [$this_function] cannot parse bedGraph input from $filename"
unless {$line =~ m/^([a-zA-Z0-9._]+)\t(\d+)\t(\d+)\t(-?\d+\.?\d*)$/};
$chr = $1; $start = $2; $end = $3, $val = $4;
#print "++ \$chr $chr ++ \$start $start ++ \$end $end ++ \$val $val\n";
$entry = Bio::ViennaNGS::BedGraphEntry->new(chromosome => $chr,
start => $start,
end => $end,
dataValue => $val);
push @{$self->data}, $entry;
}
}
sub parse_bed6_file{
my ($self,$file,$typ) = @_;
my $this_function = (caller(0))[3];
my ($line,$feat,$fc);
$file = read_file( $file, array_ref => 1, chomp =>1 );
if ($typ == 2){ # initialize an empty FeatureChain object
$fc = Bio::ViennaNGS::FeatureChain->new(type => "feature",
base => $self->base);
}
# print "********** in parse_bed6: typ= $typ ************\n";
foreach $line (@$file){
lib/Bio/ViennaNGS/FeatureIO.pm view on Meta::CPAN
else{
croak "ERROR [$this_function] don't know how to handle type $typ";
}
} #end foreach
if ($typ == 2) { push @{$self->data}, $fc; }
# print Dumper($self);
}
sub parse_bed12_file{
my ($self,$file,$typ) = @_;
my $this_function = (caller(0))[3];
my ($line,$feat,$fc);
$file = read_file( $file, array_ref => 1, chomp =>1 );
foreach $line (@$file){
my @mcData = split /\t/,$line;
if ($typ == 0){ # ArrayRef of Bio::ViennaNGS::Bed objects
my $bo = Bio::ViennaNGS::Bed->new(chromosome => $mcData[0],
start => $mcData[1],
end => $mcData[2],
name => $mcData[3],
score => $mcData[4],
lib/Bio/ViennaNGS/FeatureInterval.pm view on Meta::CPAN
is => 'rw',
isa => 'Int',
predicate => 'length',
init_arg => undef, # make this unsettable via constructor
);
with 'Bio::ViennaNGS::FeatureBase';
sub BUILD { # call a parser method, depending on $self->instanceOf
my $self = shift;
my $this_function = (caller(0))[3];
my $len = undef;
if ($self->base == 0){
confess "ERROR [$this_function] \$self->end must be > than \$self->start for 0-based start coordinates [==> start ".eval($self->start)." end ".eval($self->end)." <==]"
unless ($self->end > $self->start);
$len = eval($self->end)-eval($self->start)-1;
$self->_length($len);
}
else {
confess "ERROR [$this_function] \$self->end must be >= than \$self->start for 0-based start coordinates [==> start ".eval($self->start)." end ".eval($self->end)." <==]"
lib/Bio/ViennaNGS/Peak.pm view on Meta::CPAN
);
has 'threshold' => (
is => 'ro',
isa => 'Value',
predicate => 'has_threshold',
);
sub populate_data {
my ($self,$filep,$filen) = @_;
my $this_function = (caller(0))[3];
my ($i,$element,$chr,$start,$end,$val,$lastend);
my $have_lastend = 0;
# parse [+] strand
foreach $element (@{$filep->data}){
unless (exists ${$self->region}{$element->chromosome}{pos}{start}){
${$self->region}{$element->chromosome}{pos}{start} = $element->start;
}
if($have_lastend == 1 && $element->start>$lastend){ # fill the gap with 0s
for($i=$lastend;$i<$element->start;$i++){
lib/Bio/ViennaNGS/Peak.pm view on Meta::CPAN
}
return;
}
sub raw_peaks {
my ($self,$dest,$prefix,$log) = @_;
my ($fn,$tfn,$tfh,$maxidx,$have_pstart,$pstart,$pend,$chr);
my ($winstart, $winend, $winsum, $mean, $lastmax, $from, $to);
my ($index_of_maxwin);
my $suffix = "rawpeaks.bed";
my $this_function = (caller(0))[3];
my $flex = 10;
croak "ERROR [$this_function]: $dest does not exist\n"
unless (-d $dest);
if (defined $log){open(LOG, ">", $log) or croak "$!"}
else {croak "ERROR [$this_function] \$log not defined";}
$fn = $prefix.".".$suffix; # filename of final, sorted RAWPEAKS bed file
(undef,$tfn) = tempfile('RAWPEAKS_XXXXXXX',UNLINK=>0);
lib/Bio/ViennaNGS/Peak.pm view on Meta::CPAN
close($tfh);
close(LOG);
sortbed($tfn,$dest,$fn,0,$log);
unlink($tfn);
}
sub final_peaks {
my ($self,$dest,$prefix,$log) = @_;
my ($fn,$tfn,$tfh,$strand,$max,$idx,$val,$peak,$position,$pleft,$pright,$str);
my $suffix = "candidatepeaks.bed";
my $this_function = (caller(0))[3];
croak "ERROR [$this_function]: $dest does not exist\n"
unless (-d $dest);
if (defined $log){open(LOG, ">>", $log) or croak "$!";}
else {croak "ERROR [$this_function] \$log not defined";}
$fn = $prefix.".".$suffix;
(undef,$tfn) = tempfile('CANDIDATEPEAKS_XXXXXXX',UNLINK=>0);
lib/Bio/ViennaNGS/SpliceJunc.pm view on Meta::CPAN
# Extracts splice junctions from BED12 annotation.
#
# Writes a BED6 file for each transcript found in the BED12, listing
# all splice sites of this transcript, optionally flanking it with a
# window of +/-$window nt.
sub bed6_ss_from_bed12{
my ($bed12,$dest,$window,$can,$fastaobjR) = @_;
my ($i,$tr_name,$pos5,$pos3);
my ($splicesites,$c,$totalsj,$cansj) = (0)x4;
my @bedline = ();
my $this_function = (caller(0))[3];
croak "ERROR [$this_function] $bed12 does not exists\n"
unless (-e $bed12);
croak "ERROR [$this_function] $dest does not exist"
unless (-d $dest);
open(BED12IN, "< $bed12") or croak $!;
while(<BED12IN>){
chomp;
my ($chr,$chromStart,$chromEnd,$name,$score,$strand,$thickStart,
lib/Bio/ViennaNGS/SpliceJunc.pm view on Meta::CPAN
# Extracts splice junctions from mapped RNA-seq data
#
# Writes a BED6 file for each splice junction present in the input,
# optionally flanking it with a window of +/-$window nt. Only splice
# junctions supported by at least $mcov reads are considered.
sub bed6_ss_from_rnaseq{
my ($bed_in,$dest,$window,$mcov,$can,$fastaobjR) = @_;
my ($reads,$proper,$passed,$pos5,$pos3);
my $c = 0;
my @bedline = ();
my $this_function = (caller(0))[3];
croak "ERROR [$this_function] $bed_in does not exist\n"
unless (-e $bed_in);
croak "ERROR [$this_function] $dest does not exist"
unless (-d $dest);
open(INBED, "< $bed_in") or croak $!;
while(<INBED>){
chomp;
my ($chr, $start, $end, $info, $score, $strand) = split("\t");
lib/Bio/ViennaNGS/SpliceJunc.pm view on Meta::CPAN
# bed6_ss_to_bed12 ($bed_in,$dest,$window,$mcov)
#
# Produce BED12 from BED6 file holdig splice junctions from mapped
# RNA-seq data
sub bed6_ss_to_bed12{
my ($bed_in,$dest,$window,$mcov,$circ) = @_;
my ($reads,$proper,$passed,$pos5,$pos3,$basename,$fn,$bed12_fn);
my @result = ();
my @bedline = ();
my $this_function = (caller(0))[3];
croak "ERROR [$this_function] $bed_in does not exist\n"
unless (-e $bed_in);
croak "ERROR [$this_function] $dest does not exist"
unless (-d $dest);
$basename = fileparse($bed_in,qr/\.[^.]*/);
$fn = $basename.".bed12";
$bed12_fn = file($dest,$fn);
open(BED12OUT, "> $bed12_fn");
lib/Bio/ViennaNGS/SpliceJunc.pm view on Meta::CPAN
# Writes BED6 files for existing and novel splice junctions to $dest
# and reurns an array with the absolute path to the two resulting BED
# files
sub intersect_sj{
my ($p_annot,$p_mapped,$dest,$prefix,$window,$mil) = @_;
my ($dha,$dhm);
my $processed = 0;
my @junctions = ();
my @transcript_beds = ();
my %asj = (); # annotated splice junctions hash
my $this_function = (caller(0))[3];
my $bedtools = can_run('bedtools') or
croak "ERROR [$this_function] bedtools not found";
my $sortBed = can_run('sortBed') or
croak "ERROR [$this_function] sortBed not found";
croak "ERROR [$this_function] $p_annot does not exist\n"
unless (-d $p_annot);
croak "ERROR [$this_function] $p_mapped does not exist\n"
unless (-d $p_mapped);
lib/Bio/ViennaNGS/SpliceJunc.pm view on Meta::CPAN
# ------------------>
# 5'===]GT..........AG[====3'
#
# <-----------------
# 3'===]GA..........TG[====5'
#
sub ss_isCanonical{
my ($chr,$p5,$p3,$fo) = @_;
my ($seqL,$seqR,$pattern);
my $ss_motif_length = 2;
my $this_function = (caller(0))[3];
my $c = -1;
$seqL = $fo->stranded_subsequence($chr,$p5+1,$p5+$ss_motif_length,"+");
$seqR = $fo->stranded_subsequence($chr,$p3-$ss_motif_length+1,$p3,"+");
$pattern = sprintf("%s|%s",$seqL,$seqR);
#print STDERR "[$this_function] p5->p3 ($p5 -- $p3) $seqL -- $seqR\n";
if ($pattern =~ /^GT|AG$/) { $c = 1000;}
elsif ($pattern =~ /^CT|AC$/) { $c = 1000; }
lib/Bio/ViennaNGS/UCSC.pm view on Meta::CPAN
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 .= "/"; }
lib/Bio/ViennaNGS/UCSC.pm view on Meta::CPAN
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){
lib/Bio/ViennaNGS/UCSC.pm view on Meta::CPAN
croak "Template process failed: ", $template->error(), "\n";
if (defined $log){
print LOG "LOG Track Hub created\n";
close(LOG);
}
}
sub convert_tracks{
my ($filesdir,$genome_assembly_directory,$chromosome_size_filepath,$log) = @_;
my $this_function = (caller(0))[3];
my @bedfiles = ();
find( sub {
return unless -f;
return unless /\.bed$/;
push @bedfiles, $File::Find::name;
} , $filesdir);
if (defined $log){
open(LOG, ">>", $log) or croak $!;
lib/Bio/ViennaNGS/UCSC.pm view on Meta::CPAN
elsif ($acc =~ /^(N[CSTZ]\_[A-Z]*\d{6})$/){
return $1;
}
else {
return 0;
}
}
sub parse_fasta_header{
my $filepath = shift;
my $this_function = (caller(0))[3];
open my $file, '<', "$filepath" or die $!;
my $fastaheader = <$file>;
chomp $fastaheader;
close $file;
my @ids = ();
#>gi|556503834|ref|NC_000913.3| Escherichia coli str. K-12 substr. MG1655
if($fastaheader=~/^>gi/){
my @headerfields = split(/\|/, $fastaheader);
my $accession = $headerfields[3];
my $scientificName = $headerfields[4];
lib/Bio/ViennaNGS/Util.pm view on Meta::CPAN
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#
#^^^^^^^^^^^ Subroutines ^^^^^^^^^^#
#^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^#
sub bed_or_bam2bw {
my ($type,$infile,$chromsizes,$strand,$dest,$want_norm,$size,$scale,$log) = @_;
my ($fn_bg_tmp,$fn_bg,$fn_bw,$fn,$bn,$path,$ext,$cmd);
my @processed_files = ();
my $factor = 1.;
my $this_function = (caller(0))[3];
croak "ERROR [$this_function] \$type is '$type', however it is expected to be either 'bam' or 'bed'\n"
unless ($type eq "bam") || ($type eq "bed");
my $genomeCoverageBed = can_run('genomeCoverageBed') or
croak "ERROR [$this_function] genomeCoverageBed utility not found";
my $bedGraphToBigWig = can_run('bedGraphToBigWig') or
croak "ERROR [$this_function] bedGraphToBigWig utility not found";
my $awk = can_run('awk') or
croak "ERROR [$this_function] awk utility not found";
lib/Bio/ViennaNGS/Util.pm view on Meta::CPAN
if (defined $log){ close(LOG); }
unlink ($fn_bg_tmp);
unlink ($fn_bg);
return $fn_bw;
}
sub bed2bigBed {
my ($infile,$chromsizes,$dest,$log) = @_;
my ($bn,$path,$ext,$cmd,$outfile);
my $this_function = (caller(0))[3];
my $bed2bigBed = can_run('bedToBigBed') or
croak "ERROR [$this_function] bedToBigBed utility not found";
if (defined $log){
open(LOG, ">>", $log) or croak $!;
print LOG "LOG [$this_function] \$infile: $infile -- \$chromsizes: $chromsizes --\$dest: $dest\n";
}
croak "ERROR [$this_function] Cannot find $infile"
unless (-e $infile);
lib/Bio/ViennaNGS/Util.pm view on Meta::CPAN
}
if (defined $log){ close(LOG); }
return $outfile;
}
sub sortbed {
my ($infile,$dest,$outfile,$rm_orig,$log) = @_;
my ($cmd,$out);
my $this_function = (caller(0))[3];
my $bedtools = can_run('bedtools') or
croak "ERROR [$this_function] bedtools utility not found";
croak "ERROR [$this_function] Cannot find $infile"
unless (-e $infile);
croak "ERROR [$this_function] $dest does not exist"
unless (-d $dest);
if (defined $log){open(LOG, ">>", $log) or croak $!;}
$out = file($dest,$outfile);
lib/Bio/ViennaNGS/Util.pm view on Meta::CPAN
$kstring = "";
}
}
return( \%km );
}
sub fetch_chrom_sizes{
my $species = shift;
my %sizes;
my @chromsize;
my $this_function = (caller(0))[3];
my $test_fetchChromSizes = can_run('fetchChromSizes') or
croak "ERROR [$this_function] fetchChromSizes utility not found";
my $cmd = "fetchChromSizes $species";
my( $success, $error_message, $full_buf, $stdout_buf, $stderr_buf ) =
run(command => $cmd, verbose => 0);
if ($success){
@chromsize = @{$stdout_buf};
}
lib/Bio/ViennaNGS/Util.pm view on Meta::CPAN
my ($chr,$size)=split (/\t/,$_);
$sizes{$chr}=$size;
}
}
return(\%sizes);
}
sub mkdircheck {
my @dirstocreate=();
my $this_function = (caller(0))[3];
while(@_){
push @dirstocreate, shift(@_);
}
foreach (@dirstocreate){
my @total = split(/[\/\\]/,$_);
if ($total[0] eq '~'){
$total[0]=$HOME;
carp "WARNING [$this_function] replacing \'~\' in path with \'$HOME\'";
}
my $dir;
lib/Bio/ViennaNGS/Util.pm view on Meta::CPAN
}
my $check = $dir->stringify();
return if (-d $check);
make_path($dir,{ verbose => 0 }) or
croak "ERROR [$this_function] Cannot create directory $check";
}
}
sub rmdircheck {
my @dirstorm=();
my $this_function = (caller(0))[3];
while(@_){
push @dirstorm, shift(@_);
}
foreach (@dirstorm){
my @total = split(/[\/\\]/,$_);
my $dir;
while(@total){
$dir = dir(shift(@total)) unless (defined $dir);
$dir = dir($dir,shift(@total));
}
( run in 0.306 second using v1.01-cache-2.11-cpan-a9ef4e587e4 )