File Coverage

blib/lib/CracTools/Utils.pm

Criterion	Covered	Total	%
statement	155	239	64.8
branch	33	84	39.2
condition	12	27	44.4
subroutine	24	37	64.8
pod	24	24	100.0
total	248	411	60.3

package CracTools::Utils;
{
  $CracTools::Utils::DIST = 'CracTools';
}
# ABSTRACT: A set of useful functions
$CracTools::Utils::VERSION = '1.22';
use strict;
use warnings;

use Carp;
use Fcntl qw( SEEK_SET );


sub reverseComplement($) {
  my $dna = shift;

  # reverse the DNA sequence
  my $revcomp = reverse $dna;

  # complement the reversed DNA sequence
  $revcomp =~ tr/ACGTacgt/TGCAtgca/;
  return $revcomp;
}


sub reverse_tab($) {
  my $string = shift;
  my @tab = split(/,/,$string);
  my $newString;
  if(@tab > 0) {
    for (my $i=$#tab ; $i > 0 ; $i--){
      $newString .= $tab[$i];
      $newString .= ",";
    }
    $newString .= $tab[0];
  }
  return $newString;
}


sub isVersionGreaterOrEqual($$) {
  my ($v1,$v2) = @_;
  my @v1_nums = split(/\./,$v1);
  my @v2_nums = split(/\./,$v2);
  for(my $i = 0; $i < @v1_nums; $i++) {
    if($v1_nums[$i] >= $v2_nums[$i]) {
      return 1;
    }else {
      return 0;
    }
  }
  if(scalar @v2_nums > @v1_nums) {
    return 0;
  } else {
    return 1;
  }
}


my %conversion_hash = ( '+' => 1, '-' => '-1', '1' => '+', '-1' => '-');
sub convertStrand($) {
  my $strand = shift;
  return $conversion_hash{$strand};
}

sub removeChrPrefix($) {
  my $string = shift;
  $string =~ s/^chr//;
  return $string;
}

sub addChrPrefix($) {
  my $string = shift;
  return "chr".removeChrPrefix($string);
}


sub seqFileIterator {
  my ($file,$format) = @_;

  croak "Missing file in argument of seqFileIterator" if !defined $file;

  # Get file handle for $file
  my $fh = getReadingFileHandle($file);

  # Automatic file extension detection
  if(!defined $format) {
    if($file =~ /\.(fasta|fa)(\.|$)/) {
      $format = 'fasta';
    } elsif($file =~ /\.(fastq|fq)(\.|$)/) {
      $format = 'fastq';
    } else {
      croak "Undefined file extension";
    }
  } else {
    $format = lc $format;
  }

  # FASTA ITERATOR
  if ($format eq 'fasta') {
    # Read prev line for FASTA because we dont know the number
    # of line used for the sequence
    my $prev_line = <$fh>;
    chomp $prev_line;
    return sub {
      my ($name,$seq,$qual);
      if(defined $prev_line) {
        ($name) = $prev_line =~ />(.*)$/;
        $prev_line = <$fh>;
        # Until we find a new sequence identifier ">", we
        # concatenate the lines corresponding to the sequence
        while(defined $prev_line && $prev_line !~ /^>/) {
          chomp $prev_line;
          $seq .= $prev_line;
          $prev_line = <$fh>;
        }
        return {name => $name, seq => $seq, qual => $qual};
      } else {
        return undef;
      }
    };
  # FASTQ ITERATOR
  } elsif ($format eq 'fastq') {
    return sub {
      my ($name,$seq,$qual);
      my $line = <$fh>;
      ($name) = $line =~ /@(.*)$/ if defined $line;
      if(defined $name) {
        $seq = <$fh>;
        chomp $seq;
        <$fh>; # skip second seq name (useless line)
        $qual = <$fh>;
        chomp $qual;
        return {name => $name, seq => $seq, qual => $qual};
      } else {
        return undef;
      }
    };
  } else {
    croak "Undefined file format";
  }
}


sub pairedEndSeqFileIterator {
  my($file1,$file2,$format) = @_;

  my $it1 = seqFileIterator($file1,$format);
  my $it2 = seqFileIterator($file2,$format);

  return sub {
    my $entry1 = $it1->();
    my $entry2 = $it2->();
    if(defined $entry1 && defined $entry2) {
      return { read1 => $entry1, read2 => $entry2 };
    } else {
      return undef;
    }
  };
}


sub writeSeq {
  my ($fh,$format,$name,$seq,$qual) = @_;
  if($format eq 'fasta') {
    print $fh ">$name\n$seq\n";
  } elsif($format eq 'fastq') {
    print $fh '@'."$name\n$seq\n+\n$qual\n";
  } else {
    die "Unknown file format. Must be either a FASTA or a FASTQ";
  }
}


sub bedFileIterator {
  return getFileIterator(file => shift, type => 'bed');
}


sub gffFileIterator {
  my $file = shift;
  my $type = shift;
  return getFileIterator(file => $file, type => $type);
}


sub vcfFileIterator {
  my $file = shift;
  return getFileIterator(file => $file, type => 'vcf');
}


sub chimCTFileIterator {
  return getFileIterator(file => shift, type => 'chimCT');
}


sub bamFileIterator {
  my ($file,$region) = @_;
  $region = "" if !defined $region;

  my $fh;
  if($file =~ /\.bam$/) {
    open($fh, "-|", "samtools view $file $region" )or die "Cannot open $file, check if samtools are installed.";
  } else {
    $fh = getReadingFileHandle($file);
  }

  return sub {
    return <$fh>;
  }

}


sub getSeqFromIndexedRef {
  my ($ref_file,$chr,$pos,$length,$format) = @_;
  $format = 'fasta' if !defined $format;
  $pos++; # +1 because samtools faidx is 1-based
  # First we try without the "chr" prefix
  my $region = removeChrPrefix($chr).":$pos-".($pos+$length-1);
  my $fasta_seq = `samtools faidx $ref_file "$region"`;
  # If it has not worked, we try without
  if($fasta_seq !~ /^>\S+\n\S+/) {
    $region = addChrPrefix($chr).":$pos-".($pos+$length-1);
    $fasta_seq = `samtools faidx $ref_file "$region"`;
  }
  if($format eq 'raw') {
    my ($seq) = $fasta_seq =~ /^>\S+\n(\S+)/;
    return $seq;
  }
  return $fasta_seq;
}


sub parseBedLine {
  my $line = shift;
  my %args = @_;
  my($chr,$start,$end,$name,$score,$strand,$thick_start,$thick_end,$rgb,$block_count,$block_size,$block_starts) = split("\t",$line);
  my @blocks;

  # Manage blocks if we have a bed12
  if(defined $block_starts) {
    my @block_size = split(",",$block_size);
    my @block_starts = split(",",$block_starts);
    my $cumulated_block_size = 0;
    for(my $i = 0; $i < $block_count; $i++) {
      push(@blocks,{size        => $block_size[$i],
                   start        => $block_starts[$i],
                   end          => $block_starts[$i] + $block_size[$i],
                   block_start  => $cumulated_block_size,
                   block_end    => $cumulated_block_size + $block_size[$i],
                   ref_start    => $block_starts[$i] + $start,
                   ref_end      => $block_starts[$i] + $start + $block_size[$i],
                 });
      $cumulated_block_size += $block_size[$i];
    }
  }

  return { chr        => $chr,
    start       => $start,
    end         => $end,
    name        => $name,
    score       => $score,
    strand      => $strand,
    thick_start => $thick_start,
    thick_end   => $thick_end,
    rgb         => $rgb,
    blocks      => \@blocks,
  };
}


sub parseGFFLine {
  my $line = shift;
  my $type = shift;
  my $attribute_split;
  if($type =~ /gff3/i) {
    $attribute_split = sub {my $attr = shift; return $attr =~ /(\S+)=(.*)/;};
  } elsif ($type eq 'gtf' || $type eq 'gff2') {
    $attribute_split = sub {my $attr = shift; return $attr  =~ /(\S+)\s+"(.*)"/;};
  } else {
    croak "Undefined GFF format (must be either gff3,gtf, or gff2)";
  }
  my($chr,$source,$feature,$start,$end,$score,$strand,$frame,$attributes) = split("\t",$line);
  my @attributes_tab = split(";",$attributes);
  my %attributes_hash;
  foreach my $attr (@attributes_tab) {
    my ($k,$v) = $attribute_split->($attr);
    $attributes_hash{$k} = $v;
  }
  return { chr        => $chr,
    source     => $source,
    feature    => $feature,
    start      => $start,
    end        => $end,
    score      => $score,
    strand     => $strand,
    frame      => $frame,
    attributes => \%attributes_hash,
  };
}


sub parseVCFLine {
  my $line = shift;
  my($chr,$pos,$id,$ref,$alt,$qual,$filter,$info) = split("\t",$line);
  my @alts = split(",",$alt);
  my %infos;
  my @infos = split(";",$info);
  foreach (@infos) {
    my($k,$v) = split("=",$_);
    $infos{$k} = $v;
  }

  return { chr => $chr,
    pos     => $pos,
    id      => $id,
    ref     => $ref,
    alt     => \@alts,
    qual    => $qual,
    filter  => $filter,
    info    => \%infos,
  };
}


sub parseChimCTLine {
  my $line = shift;
  my($id,$name,$chr1,$pos1,$strand1,$chr2,$pos2,$strand2,$chim_value,$spanning_junction,$spanning_PE,$class,$comments,@others) = split("\t",$line);
  my($sample,$chim_key) = split(":",$id);
  my @comments = split(",",$comments);
  my %comments;
  foreach my $com (@comments){
      my ($key,$val) = split("=",$com);
      $comments{$key} = $val;
  }
  my %extend_fields;
  foreach my $field (@others) {
    my ($key,$val) = split("=",$field);
    if(defined $key && defined $val) {
      $extend_fields{$key} = $val;
    }
  }
  return {
    sample            => $sample,
    chim_key          => $chim_key,
    chr1              => $chr1,
    pos1              => $pos1,
    strand1           => $strand1,
    chr2              => $chr2,
    pos2              => $pos2,
    strand2           => $strand2,
    chim_value        => $chim_value,
    spanning_junction => $spanning_junction,
    spanning_PE       => $spanning_PE,
    class             => $class,
    comments          => \%comments,
    extended_fields     => \%extend_fields,
  };
}


sub parseSAMLineLite {
  my $line = shift;
  my ($qname,$flag,$rname,$pos,$mapq,$cigar,$rnext,$pnext,$tlen,$seq,$qual,@others) = split("\t",$line);
  my @cigar_hash = map { { op => substr($_,-1), nb => substr($_,0,length($_)-1)} } $cigar =~ /(\d+\D)/g;
  return {
    qname => $qname,
    flag => $flag,
    rname => $rname,
    pos => $pos,
    mapq => $mapq,
    cigar => \@cigar_hash,
    original_cigar => $cigar,
    rnext => $rnext,
    pnext => $pnext,
    tlen => $tlen,
    seq => $seq,
    qual => $qual,
    extended_fields => \@others,
  };
}



sub getFileIterator {
  my %args = @_;
  my $file = $args{file};
  my $type = $args{type};
  my $skip = $args{skip};
  my $header_regex = $args{header_regex};
  my $parsing_method = $args{parsing_method};
  my @parsing_arguments = ();


  croak "Missing arguments in getFileIterator" if !defined $file;

  if(!defined $parsing_method && defined $type) {
    if($type =~ /gff3/i || $type eq 'gtf' || $type eq 'gff2') {
      $header_regex = '^#';
      $parsing_method = sub { return parseGFFLine(@_,$type) };
      push (@parsing_arguments,$type);
    } elsif($type =~ /bed/i) {
      $header_regex = '^track\s';
      $parsing_method = \&parseBedLine;
    } elsif ($type =~ /vcf/i) {
      $header_regex = '^#';
      $parsing_method = \&parseVCFLine;
    } elsif ($type =~ /chimCT/i) {
      $header_regex = '^#';
      $parsing_method = \&parseChimCTLine;
    } elsif ($type =~ /SAM/i || $type =~ /BAM/i) {
      $header_regex = '^@';
      $parsing_method = \&parseSAMLineLite;
    } else {
      croak "Undefined format type";
    }
  }

  # set defaults
  #$header_regex = '^#' if !defined $header_regex;
  $parsing_method = sub{ return { line => shift } } if !defined $parsing_method;

  # We get a filehandle on the file
  my $fh = getReadingFileHandle($file);

  # $curr_pos will hold the SEEK_POS of the current_line
  my $curr_pos = tell($fh);
  # $line will hold the content of the current_line
  my $line = <$fh>;


  # if we need to skip lines
  if(defined $skip) {
    for(my $i = 0; $i < $skip; $i++) {
      $curr_pos = tell($fh);
      $line = <$fh>;
    }
  }

  # Skip line that match a specific regex
  if(defined $header_regex) {
    while($line && $line =~ /$header_regex/) {
      $curr_pos = tell($fh);
      $line = <$fh>;
    }
  }

  # The iterator itself (an unnamed subroutine
  return sub {
    if (defined $line) {
      chomp $line;
      # We parse the line with the appropriate methd
      my $parsed_line = $parsing_method->($line,@parsing_arguments);
      # Add seek pos to the output hashref
      $parsed_line->{seek_pos} = $curr_pos;
      $curr_pos = tell($fh);
      $line = <$fh>; # Get next line
      return $parsed_line;
    } else {
      return undef;
    }
  };
}


sub getReadingFileHandle {
  my $file = shift;
  my $fh;
  if($file =~ /\.gz$/) {
    open($fh,"gunzip -c $file |") or die ("Cannot open $file");
  } else {
    open($fh,"< $file") or die ("Cannot open $file");
  }
  return $fh;
}


sub getWritingFileHandle {
  my $file = shift;
  my $fh;
  if($file =~ /\.gz$/) {
    open($fh,"| gzip > $file") or die ("Cannot open $file");
  } else {
    open($fh,"> $file") or die ("Cannot open $file");
  }
  return $fh;
}


sub getLineFromSeekPos {
  my($fh,$seek_pos) = @_;
  seek($fh,$seek_pos,SEEK_SET);
  my $line = <$fh>;
  chomp($line);
  return $line;
}

1;

__END__