rf-correlate

#!/usr/bin/env perl

##
# RF Correlate
# RNA Framework [http://www.rnaframework.com]
#
# Author:  Danny Incarnato (dincarnato[at]rnaframework.com)
# Summary: Calculates correlation between datasets
#
# This program is free software, and can be redistribute  and/or modified
# under the terms of the GNU General Public License as published by
# the Free Software Foundation, either version 3 of the License, or
# any later version.
#
# Please see <http://www.gnu.org/licenses/> for more informations.
##

use strict;
use Config;
use File::Basename;
use FindBin qw($Bin);
use Getopt::Long qw(:config no_ignore_case);
use threads;
use threads::shared;

use lib $Bin . "/lib";

use Core::Mathematics qw(:all);
use Core::Statistics;
use Core::Utils;
use Data::Sequence::Utils;
use RF::Data::IO::RC;
use RF::Data::IO::XML;
use Term::Constants qw(:screen);
use Term::Table;

die "\n  [!] Error: This program requires ithreads." .
    "\n             Please recompile Perl with ithreads and try again\n\n" unless(exists $Config{useithreads});

$|++;

my ($output, $overwrite, $help, $minvalues,
    $threads, $table, $skipoverall, $spearman,
    $ignoreseq, $singlefile, $isRC, $keepbases,
    $coverage, $ratio, $minCov, $maxReact,
    $capReact, $medianCov, @pool, @files, 
    @overall, %files);

my @ids : shared;
my @results : shared;
my @allreactivities1 : shared;
my @allreactivities2 : shared;
my %results : shared;
%results = ( diffseq     => 0,
             nominvalues => 0,
             failed      => 0,
             correlated  => 0,
             missing     => 0 );

do {

    local $SIG{__WARN__} = sub { };

    GetOptions( "h|help"               => \$help,
                "o|output=s"           => \$output,
                "ow|overwrite"         => \$overwrite,
                "m|min-values=s"       => \$minvalues,
                "p|processors=i"       => \$threads,
                "s|skip-overall"       => \$skipoverall,
                "S|spearman"           => \$spearman,
                "kb|keep-bases=s"      => \$keepbases,
                "i|ignore-sequence"    => \$ignoreseq,
                "c|coverage"          => \$coverage,
                "r|ratio"              => \$ratio,
                "ec|median-coverage=s" => \$medianCov,
                "mc|min-coverage=i"    => \$minCov,
		        "mr|max-react=s"       => \$maxReact,
                "cr|cap-react=s"       => \$capReact ) or help(1);

    @files = uniq(@ARGV);

};

help() if ($help);

# Default
$output ||= "rf_correlate.txt";
$threads ||= 1;
$keepbases ||= "ACGT";
$minCov ||= 1;
$capReact ||= 1e9;

##
# Input validation
##

die "\n  [!] Error: No output directory specified\n\n" unless(defined $output);
die "\n  [!] Error: 2 XML directories/files required\n\n" if (@files < 2);
die "\n  [!] Error: Number of processors must be an integer greater than 0\n\n" if ($threads < 1);
die "\n  [!] Error: Miminum number of values must be greater than 0\n\n" if (defined $minvalues &&
                                                                             (!isnumeric($minvalues) ||
                                                                              $minvalues == 0));
die "\n  [!] Error: Parameters -r and -c are mutually exclusive\n\n" if ($ratio && $coverage);
die "\n  [!] Error: Invalid bases to keep\n\n" if (!isiupac($keepbases) && $keepbases !~ m/^all$/i);
die "\n  [!] Error: Maximum reactivity must be positive\n\n" if (defined $maxReact && !ispositive($maxReact));
die "\n  [!] Error: Median coverage must be positive\n\n" if (defined $medianCov && !ispositive($medianCov));
die "\n  [!] Error: Reactivity cap must be positive\n\n" if (!ispositive($capReact));

$minvalues = round($minvalues) if ($minvalues > 1);
$singlefile = 1 if (!-d $files[0] && !-d $files[1]);
$keepbases = $keepbases =~ m/^all$/i ? "ACGT" : join("", sort(uniq(split("", join("", iupac2nt(rna2dna(uc($keepbases))))))));

print "\n[+] Importing input XML directories\/files...";

for (@files) {

    die "\n\n  [!] Error: Provided XML directory\/file doesn't exist\n\n" if (!-e $_);

    if (-d $_) {

        $_ =~ s/\/?$/\//;

        opendir(my $dh, $_) or die "\n\n  [!] Error: Unable to read from input directory \"" . $_ . "\" (" . $! . ")\n\n";
        while(my $file = readdir($dh)) {

            next if ($file !~ m/\.xml$/);

            $file =~ s/\.xml//;
            $files{$file}++;

        }
        close($dh);

    }
    else {

        $files{(fileparse($_, qw(.rc .xml)))[0]}++;
        $isRC++ if ($_ =~ m/\.rc$/);

    }

}

# Only one of the two files is in RC format
die "\n\n  [!] Error: Mixed file types provided\n\n" if ($isRC == 1);
die "\n\n  [!] Error: Provided RC files have unequal sizes\n\n" if ($isRC && uniq(map { -s $_ } @files) != 1);

# Delete all transcripts that are not common to all XML directories
if (!$singlefile) {

    for (keys %files) { delete($files{$_}) if ($files{$_} != @files); }

    print " " . (keys %files) . " common transcripts.";

    die "\n\n  [!] Error: No common transcript ID found between XML directories\/files\n\n" unless(keys %files);

}

if (-e $output) {

    if ($overwrite) { unlink($output) or die  "\n\n  [!] Error: Unable to overwrite output file (" . $! . ")\n\n"; }
    else { die "\n\n  [!] Error: Output file already exists." .
               "\n      Please use -ow (or --overwrite) to overwrite output file\n\n"; }

}

print "\n[+] Calculating correlations [Last: none]";

if ($isRC) {

    my $io = RF::Data::IO::RC->new(file => $files[0]);
    @ids = $io->ids();

}
else { @ids = $singlefile ? (keys %files)[0] : keys %files; }

@pool = map{ threads->create(\&correlate) } 1 .. $threads;
$_->join() for(@pool);

if ($singlefile && !$isRC) {

    if (@results) { print "\n[+] Correlation: " . sprintf("%.3f", $results[0]->[1]) . " (p-value: " .  sprintf("%.2e", $results[0]->[2]) . ")"; }
    else { print "\n\n  [!] Error: Correlation calculation failed"; }

}
else {

    if (!$skipoverall && @allreactivities1) {

        @overall = _correlate(\@allreactivities1, \@allreactivities2);

        print "\n[+] Overall correlation (over " . scalar(@allreactivities1) . " bases): " . sprintf("%.3f", $overall[0]) . " (p-value: " . sprintf("%.2e", $overall[1]) . ")";

    }

    if (@results) {

        my $top = min(scalar(@results), 10);
        @results = sort { $b->[1] <=> $a->[1] } @results;

        print "\n[+] Top $top correlated transcripts:\n\n";

        $table = Term::Table->new(indent => 2);
        $table->head("Transcript", "Correlation", "p-value");
        $table->row($_->[0], sprintf("%.3f", $_->[1]), sprintf("%.2e", $_->[2])) for (@results[0 .. $top - 1]);
        $table->print();

        open(my $wh, ">", $output) or die "\n\n  [!] Error: Unable to write output file (" . $! . ")\n\n";
        select((select($wh), $|=1)[0]);
        print $wh "#Transcript\tCorrelation\tp-value\n";
        print $wh join("\t", @$_) . "\n" for (@results);
        close($wh);

    }
    else { print "\n\n  [!] Error: Correlation calculation failed for all transcripts"; }

}

print "\n\n[+] Correlation statistics:\n" .
      "\n  [*] Combined transcripts:  " . $results{correlated} .
      "\n  [*] Discarded transcripts: " . $results{failed} . " total" .
      "\n                             " . ($results{failed} - $results{diffseq} - $results{nominvalues}) . " parsing failed" .
      "\n                             " . $results{diffseq} . " mismatch between transcript sequences" .
      "\n                             " . $results{nominvalues} . " not enough values for correlation calculation";
print "\n                             " . $results{missing} . " transcripts absent in second RC file" if ($isRC);

print "\n\n[+] All done.\n\n";

sub correlate {

    my ($io1, $io2);

    if ($isRC) {

        $io1 = RF::Data::IO::RC->new(file => $files[0]);
        $io2 = RF::Data::IO::RC->new(file => $files[1]);

    }

    while (1) {

        my ($id, $entry1, $entry2, $nbases,
            $sequence, @reactivity1, @reactivity2,
            @commonidx, @pearson);

        { lock(@ids);
          $id = shift(@ids) if (@ids); }

        last unless($id);

        eval {

            no warnings;

            if (!$isRC) {

                $entry1 = RF::Data::IO::XML->new(file => -d $files[0] ? $files[0] . $id . ".xml" : $files[0]);
                $entry2 = RF::Data::IO::XML->new(file => -d $files[1] ? $files[1] . $id . ".xml" : $files[1]);

            }

        };

        if ($@) { # Exception from eval

            lock(%results);
            $results{failed}++;

            undef($@);

            next;

        }

        if ($isRC) {

            $entry1 = $io1->read($id);
            $entry2 = $io2->read($id);

            if (!$entry2) {

                lock(%results);
                $results{missing}++;
                $results{failed}++;

                undef($@);

                next;

            }

            if ($medianCov && 
                ($entry1->mediancoverage() < $medianCov || $entry2->mediancoverage() < $medianCov)) {

                lock(%results);
                $results{nominvalues}++;
                $results{failed}++;

                undef($@);

                next;

            }

        }

        if ((!$ignoreseq && $entry1->sequence() ne $entry2->sequence()) ||
            ($ignoreseq && $entry1->length() != $entry2->length())) {

            lock(%results);
            $results{diffseq}++;
            $results{failed}++;

            undef($@);

            next;

        }

        $keepbases ||= $entry1->reactive() if (!$isRC);
        $sequence = $entry1->sequence();
        ($nbases) = $sequence =~ s/([$keepbases])/$1/ge;

        if ($isRC) {

            my (@counts1, @counts2, @cov1, @cov2);
            @counts1 = $entry1->counts();
            @counts2 = $entry2->counts();
            @cov1 = $entry1->coverage();
            @cov2 = $entry2->coverage();

            if ($ratio) {

		        my @i = grep { $cov1[$_] } 0 .. $#cov1;
                @i = grep { $cov2[$_] } @i;
		        @counts1 = @counts1[@i];
                @counts2 = @counts2[@i];
                @cov1 = @cov1[@i];
                @cov2 = @cov2[@i];

	        }

            @reactivity1 = $coverage ? @cov1 : ($ratio ? map { $counts1[$_] / $cov1[$_] } 0 .. $#cov1 : @counts1);
            @reactivity2 = $coverage ? @cov2 : ($ratio ? map { $counts2[$_] / $cov2[$_] } 0 .. $#cov2 : @counts2);

            @commonidx = grep { substr($sequence, $_, 1) =~ m/^[$keepbases]$/ &&
                                $cov1[$_] >= $minCov && $cov2[$_] >= $minCov && sum($cov1[$_], $cov2[$_]) > 0 } 0 .. $#cov1;

        }
        else {

            @reactivity1 = map { isnan($_) ? "NaN" : $_ } $entry1->reactivity();
            @reactivity2 = map { isnan($_) ? "NaN" : $_ } $entry2->reactivity();

            @commonidx = grep { isnumeric($reactivity1[$_]) &&
                                isnumeric($reactivity2[$_]) } 0 .. $#reactivity1;

        }

        if ($maxReact && (($isRC && $ratio) || !$isRC)) {

            @commonidx = grep { $reactivity1[$_] <= $maxReact &&
                                $reactivity2[$_] <= $maxReact } @commonidx;

        }

        if (defined $minvalues || @commonidx <= 2) {

            if (($minvalues < 1 && @commonidx / $nbases < $minvalues) ||
                ($minvalues >= 1 && @commonidx < $minvalues) ||
                @commonidx <= 2) {

                lock(%results);
                $results{nominvalues}++;
                $results{failed}++;

                undef($@);

                next;

            }

        }

        @reactivity1 = @reactivity1[@commonidx];
        @reactivity2 = @reactivity2[@commonidx];
        @pearson = _correlate(\@reactivity1, \@reactivity2);

        { lock(%results);
          lock(@results);
          $results{correlated}++;
          push(@results, shared_clone([$id, @pearson]));

          if (!$skipoverall) {

              push(@allreactivities1, @reactivity1);
              push(@allreactivities2, @reactivity2);

          } }

        print CLRRET . "[+] Calculating correlations [Last: $id]";

    }

    threads->exit();

}

sub _correlate {

    local $SIG{__WARN__} = sub { }; # Suppresses warnings from Pearson and Spearman when stdev = 0

    return($spearman ? spearman(@_, { cap => $capReact }) : pearson(@_, { cap => $capReact }));

}

sub help {

    print "\n  [!] Error: Invalid option. Please check the help\n" if ($_[0]);

    die <<HELP;

 RF Correlate (v$Core::Utils::VERSION)
 RNA Framework [http://www.rnaframework.com]

 Author:  Danny Incarnato (dincarnato[at]rnaframework.com)
 Summary: Calculates pairwise correlations of structure probing experiments

 Usage:   rf-correlate [Options] XML_folder_rep1/ XML_folder_rep2/           # Whole transcriptome
          rf-correlate [Options] file_rep1.rc file_rep2.rc                   # Whole transcriptome
          rf-correlate [Options] file_rep1.xml file_rep2.xml                 # Single transcript

 Options                                  Description
 -p  or --processors       <int>          Number of processors to use (Default: 1)
 -o  or --output           <string>       Output TSV file (Default: rf_correlate.txt)
 -ow or --overwrite                       Overwrites output file (if the specified file already exists)
 -m  or --min-values       <float>        Minimum number of values to calculate correlation (Default: off)
                                          Note: if a value between 0 and 1 is provided, this is interpreted as a
                                                fraction of the transcript's length
 -cr or --cap-react        <float>        Maximum reactivity value to cap reactivities to (>0, Default: 1e9)
                                          Note: if processing RC files, this parameter only applies to ratios (-r)
 -mr or --max-react        <float>        Reactivity values above this threshold will be excluded from correlation
                                          calculation (>0, Default: none)
                                          Note: if processing RC files, this parameter only applies to ratios (-r)
 -s  or --skip-overall                    Skips overall experiment correlation calculation (faster)
 -i  or --ignore-sequences                Ignores sequence differences (e.g. SNVs) between the compared transcripts
 -S  or --spearman                        Uses Spearman instead of Pearson to calculate correlation

 RC file-specific options
 -kb or --keep-bases       <string>       Bases on which correlation should be calculated (Default: all)
                                          Note: this option has effect only on RC files. For XML files, reactive
                                                bases are automatically identified from the "reactive" attribute
 -mc or --min-coverage     <int>          Restricts the correlation analysis to bases exceeding this coverage
 -ec or --median-coverage  <float>        Restricts the correlation analysis to transcripts with median coverage
                                          above this threshold (>=0, Default: 0)
 -c  or --coverage                        Correlation is calculated on the coverage, rather than on the raw RT
                                          stop/mutation counts
 -r  or --ratio                           Correlation is calculated on the ratio between, the RT stop/mutation
                                          counts and the coverage, rather than on the raw RT stop/mutation counts

HELP

}