annotate bin/ppp.pm @ 13:98b4a5ec8e63 draft

Uploaded

package ppp;

use strict;
use warnings;
use FindBin;
use lib $FindBin::Bin;
use Rcall qw ( histogram );
use Math::CDF;

use Exporter;
our @ISA = qw( Exporter );
our @EXPORT_OK = qw( &ping_pong_partners );

sub ping_pong_partners
{
  my ( $TE_fai, $sam, $dir, $max ) = @_;
  my ( $hashRef, $dupRef ) = count_mapped ( $TE_fai, $sam );
  my ( %num_per_overlap_size, $overlap_number, $reverseR, $begRev, $endRev, $sensR, $begSens, $endSens, $snum, $rnum, $overlap );
  my ( $SP, $AP, $SN, $AN, $txt); 
  my $flag = 0;
  my @distri_overlap = (); my @overlaps_names = ();

  open my $ppp_f, '>', $dir."ppp.txt" || die "cannot create ppp.txt $!\n";   
  foreach my $k ( sort keys %{$hashRef} )
  {
    my $v = $hashRef->{$k};
    my $TE_dir = $dir.$k.'/';

    %num_per_overlap_size = (); $overlap_number = 0;      
    $flag = 0;
    for ( my $i = 0; $i <= $#{$v->[1]} ; $i++ )
    {
      $reverseR = ${$v->[1]}[$i] ;
      $begRev = $reverseR->[0];
      $endRev = $begRev + length($reverseR->[1]) - 1;

      my $revR = reverse($reverseR->[1]);
      $revR =~ tr/atgcuATGCU/tacgaTACGA/;

      for ( my $j = 0; $j <= $#{$v->[0]}; $j++ )
      {
        $sensR = ${$v->[0]}[$j];
        $begSens = $sensR->[0];

        $endSens = $begSens + length($sensR->[1]) - 1;

        if ( $begSens <= $endRev && $endSens > $endRev )
        {
          $flag = 1;
          
          mkdir $TE_dir;
          open  $SP, '>>', $TE_dir."sensPPP.txt" || die "cannot create sensPPP\n";
          open  $AP, '>>', $TE_dir."antisensPPP.txt" || die "cannot create antisensPPP\n";      
          open  $SN, '>>', $TE_dir."sens.txt"  || die "cannot create sens\n";
          open  $AN, '>>', $TE_dir."antisens.txt"  || die "cannot create antisens\n";
          open  $txt, '>', $TE_dir.'overlap_size.txt' || die "cannot open repartition\n";      
          
          $overlap = $endRev - $begSens + 1;
          $snum =  $dupRef->{$sensR->[0].$sensR->[1].$sensR->[2].$sensR->[3]};
          $rnum = $dupRef->{$reverseR->[0].$reverseR->[1].$reverseR->[2].$reverseR->[3]};

          if ( $overlap == 10 )
          {
           print $SP ">$sensR->[0]|$sensR->[2]|$sensR->[3]|$snum\n$sensR->[1]\n";
           print $AP ">$reverseR->[0]|$reverseR->[2]|$reverseR->[3]|$rnum\n$revR\n";
          }
          else
          {
            print $SN ">$sensR->[0]|$sensR->[2]|$sensR->[3]|$snum\n$sensR->[1]\n";
            print $AN ">$reverseR->[0]|$reverseR->[2]|$reverseR->[3]|$rnum\n$revR\n";
          }
          next if $overlap > $max;
          if ( $snum < $rnum )
          {
            $num_per_overlap_size{$overlap} += $snum;
            $overlap_number += $snum;
          }
          else
          {
            $num_per_overlap_size{$overlap} += $rnum ;
            $overlap_number += $rnum ;
          }
        }
      }
    }

    if ( $max != 0 )
    {
      my @overlaps = (); 
      push @overlaps_names, $k;  
      for my $i (1..$max)
      { 
        $num_per_overlap_size{$i} = 0 unless exists( $num_per_overlap_size{$i} );
        push @overlaps, $num_per_overlap_size{$i};
      }
      push @distri_overlap, \@overlaps;
    } 

    if ( $flag == 1 )
    {
      my $histo_png = $TE_dir.'histogram.png';
      histogram( \%num_per_overlap_size, $histo_png, $overlap_number );
      print $txt "size\tnumber\tpercentage of the total overlap number\n";
      foreach my $k ( sort {$a <=> $b} keys %num_per_overlap_size )
      {
        my $percentage = 0;
        $percentage = $num_per_overlap_size{$k} * 100 / $overlap_number unless $overlap_number == 0;
        print $txt "$k\t$num_per_overlap_size{$k}\t"; printf $txt "%.2f\n",$percentage;   
          
      }
      close $txt;
    }

  }

  foreach my $tabP (  @distri_overlap )
  {
    my $sum = sum($tabP);
    my $ten = $tabP->[9];
    my $mean = mean($tabP);
    my $std = standard_deviation($tabP, $mean);
    my $zsc = z_significance($ten, $mean, $std);
    my $name = shift @overlaps_names;
    my $prob = 'NA';
    $prob =  1 - &Math::CDF::pnorm( $zsc ) if $zsc ne 'NA';
    print $ppp_f (join ("\t", $name, $sum, $ten, $mean, $std, $zsc, $prob ),"\n" );
  }
  close $ppp_f; 
} 

sub count_mapped
{
  my ( $fai, $in_file ) = @_;
  my ( %mapped, %dup );

  open my $f, '<', $fai || die "cannot open $fai $! \n";
  while(<$f>)
  {
    if ($_ =~ /(.*)\t(\d+)\n/)
    {
      $mapped{$1} = [];
      $mapped{$1}->[0] = []; $mapped{$1}->[1] = [];        
    }
  }
  close $f;

  open my $infile, "samtools view  $in_file |"|| die "cannot open input file $! \n";
  while(<$infile>)
  {
    unless ($_ =~ /^\@[A-Za-z][A-Za-z](\t[A-Za-z][A-Za-z0-9]:[ -~]+)+$/ || $_ =~ /^\@CO\t.*/ )
    {    
      my @line = split (/\t/,$_);
      if ($line[1] == 0)
      {
        push @{$mapped{$line[2]}->[0]} , [$line[3], $line[9], $line[1],  $line[2]] unless exists ($dup{$line[3].$line[9].$line[1].$line[2]});
        $dup{$line[3].$line[9].$line[1].$line[2]}+=1;
      }
      elsif ($line[1] == 16)
      {
        push @{$mapped{$line[2]}->[1]} , [$line[3], $line[9], $line[1],  $line[2]] unless exists ($dup{$line[3].$line[9].$line[1].$line[2]});
        $dup{$line[3].$line[9].$line[1].$line[2]}+=1;
      }
    }
  }
  close $infile;
  return (\%mapped, \%dup );
}

sub sum
{
  my $arrayref = shift;
  my $result = 0;
  foreach (@$arrayref) {$result += $_}
  return $result;
}

sub mean
{
  my $arrayref = shift;
  my $result;
  foreach (@$arrayref) {$result += $_}
  return $result / scalar(@$arrayref);
}

sub standard_deviation
{
  my ($arrayref, $mean) =  @_;
  return sqrt ( mean ( [map $_**2 , @$arrayref ]) - ($mean**2));  
}

sub z_significance
{
  my ($ten, $mean, $std) = @_;
  my $z = 'NA';
  $z = (($ten - $mean) / $std) if $std != 0;
  return $z;
}

1;