maq 0.7.1-5 / usr / share / maq / paf_utils.pl

#!/usr/bin/perl -w

# Contact: lh3

use strict;
use warnings;
use Getopt::Std;

# PAF: Pairwise Alignment Format:
#   qName qBegin strand tName tBegin score CIGAR qSeq/. qQual/. meta
#
# The `score' field can be any non-negative integer but ideally is a
# Phred-scaled error probability score.

my $version = '0.3.2';
my %cmd_hash = (eland2paf=>\&eland2paf, mapview=>\&mapview,
				soap2paf=>\&soap2paf, pafeval=>\&pafeval, maq2paf=>\&maq2paf,
				rmap2paf=>\&rmap2paf, mapview2map=>\&mapview2map,
				shrimp2paf=>\&shrimp2paf, novo2paf=>\&novo2paf, help=>\&help);
die(qq/
Program: paf_utils.pl (Utilities for PAF - the Pairwise Alignment Format)
Version: $version
Contact: Heng Li <lh3\@sanger.ac.uk>

Usage:   paf_utils.pl <command> <arguments>

Command: eland2paf      convert eland alignment format to PAF
         soap2paf       convert soap alignment format to PAF
         rmap2paf       convert rmap alignment format to PAF
         maq2paf        convert maq alignment format to PAF
         shrimp2paf     convert SHRiMP alignment to PAF
         pafeval        evaluate alignment from MAQ simulation
         mapview        view the maq alignment (.map)
         mapview2map    convert mapview to binary .map format
         help           about PAF format
\n/) if (@ARGV == 0);

my $cmd = shift(@ARGV);
die(qq/Unknown command "$cmd".\n/) unless (defined($cmd_hash{$cmd}));
&{$cmd_hash{$cmd}};
exit;

sub help {
print qq{
PAF stands for Pairwise Alignment Format. I design this format to
facilitate the comparison between different alignment programs.  PAF is
TAB delimited. It consists of nine predefined fields and flexible number
of additional fields specific to the alignment program.

In PAF, each line must contain the following 9 fields:

  qName qBegin strand tName tBegin score CIGAR qSeq/. qQual/.

In this line, `q' stands for query, `t' stands for target, `strand' is
either + or -, and `score' is an integer which is ideally the Phred-like
quality of the alignment reliability. Both `qBegin' and `tBegin' are the
leftmost position of the alignment on the query and target,
respectively. The rightmost position of the alignment can be calculated
in combination with `CIGAR'.  `qSeq' is the full (not padded) sequence
of the query. It is presented on the same strand as the target sequence,
no matter whether `strand' is + or -. `qQual' is on the same strand as
`qSeq'. Either can be a dot if sequence or quality is not available.

The 10th column is the name of alignment program. The information of
additional fields is determined by the alignment program:

  eland count0 count1 count2 [mismatch [...]]
  maq mateDist flag snglMapQ altMapQ nMis sumErr seedCount0 seedCount1
  soap a/b flag [mismatch [...]]
  rmap nMismatches

For eland, `score' is estimated by assuming uniform base qualities.
Scores are not available for both soap and rmap.

};
}

sub pafeval {
  my %opts;
  getopts('p', \%opts);
  die("Usage: paf_utils.pl pafeval [-p] <in.paf>\n") if (@ARGV == 0 && -t STDIN);
  my (@c0, @c1);
  my $max_q = 0;
  while (<>) {
	my @t = split;
	my $q = int($t[5]/10);
	my $is_correct = 1;
	my $l = length($t[7]); # is the read length
	if ($t[7] eq '.') { # no read sequence
	  my $tmp = $_;
	  $l = $t[1]-1; $_ = $t[6];
	  s/(\d+)[MI]/$l+=$1,'x'/eg;
	  $_ = $tmp;
	}
	$max_q = $q if ($q > $max_q);
	++$c0[$q];
	if ($t[0] =~ /^([^_\s]+)_(\d+)_(\d+)/) {
	  if ($1 ne $t[3]) {
		$is_correct = 0;
	  } elsif ($2 < $3) {
		if ($t[2] eq '+') {
		  $is_correct = 0 if (abs($t[4]-$t[1]+1 - $2) > 1);
		} else {
		  $is_correct = 0 if (abs($t[4]-$t[1]+1+$l-1 - $3) > 3); # why not zero? because of the indels
		}
	  } else {
		if ($t[2] eq '+') {
		  $is_correct = 0 if (abs($t[4]-$t[1]+1 - $3) > 1);
		} else {
		  $is_correct = 0 if (abs($t[4]-$t[1]+1+$l-1 - $2) > 3);
		}
	  }
	} else {
	  die("[pafeval] the reads seem not to be generated by MAQ.\n");
	}
	++$c1[$q] unless ($is_correct);
	print $_ if (defined($opts{p}) && !$is_correct);
  }
  my ($cc0, $cc1) = (0, 0);
  for (my $i = $max_q; $i >= 0; --$i) {
	$c0[$i] = 0 unless (defined $c0[$i]);
	$c1[$i] = 0 unless (defined $c1[$i]);
	$cc0 += $c0[$i]; $cc1 += $c1[$i];
	printf("%.2dx %12d / %-12d  %12d  %.3e\n", $i, $c1[$i], $c0[$i], $cc0, $cc1/$cc0);
  }
}

sub eland2paf {
  my %opts = (q=>25);
  getopts('sq:', \%opts);
  die("Usage: paf_utils.pl eland2paf [-s] [-q 25] <in.eland>\n") if (@ARGV == 0 && -t STDIN);
  my $q = $opts{q}; # average base quality
  my @logn;
  $logn[$_] = int(3.434 * log($_) + 0.5) for (1 .. 255);
  my $fh = \*STDOUT;
  open($fh, "| sort -k4,4 -k5,5n -k3,3") if (defined($opts{s}));
  while (<>) {
	my @t = split;
	next unless ($t[2] =~ /^U/);
	my $size = length($t[1]);
	$t[0] =~ s/^>//;
	$t[6] =~ s/\.fa((sta)?)$//i; # in most cases, we do not want to see the file extension
	$t[8] =~ tr/FR/+-/;
	if ($t[8] eq '-') { # reverse
	  $t[1] = reverse($t[1]);
	  $t[1] =~ tr/acgtACGTswmryk/tgcaTGCAswkyrs/;
	}
	my $score = 0;
	if ($t[3] == 1) {
	  if ($t[4]) {
		$score = $q - $logn[$t[4]];
	  } elsif ($t[5]) {
		$score = $q * 2 - $logn[$t[5]];
	  } else {
		$score = $q * 3;
	  }
	} elsif ($t[4] == 1) {
	  $score = $t[5]? $q - $logn[$t[5]] : 2 * $q;
	} else { # $t[5] == 1
	  $score = $q;
	}
	print $fh join("\t", $t[0], 1, @t[8,6,7], $score, $size."M", $t[1], '.', 'eland', @t[3,4,5]);
	print "\t", join("\t", @t[10..$#t]) if ($t[3] != 1);
	print "\n";
  }
  close($fh) if (defined($opts{s}));
}

sub novo2paf {
  my %opts;
  getopts('s', \%opts);
  die("Usage: paf_utils.pl soap2paf [-s] <in.soap>\n") if (@ARGV == 0 && -t STDIN);
  my $fh = \*STDOUT;
  open($fh, "| sort -k4,4 -k5,5n -k3,3") if (defined($opts{s}));
  while (<>) {
	next if (/^#/);
	my @t = split("\t");
	next if ($t[4] ne 'U');
	$t[0] = substr($t[0], 1); # trim @/>
	$t[7] = substr($t[7], 1); # trim >
	$t[9] = ($t[9] =~ /F/)? '+' : '-';
	$t[3] = "." unless ($t[3]);
	my $cigar = length($t[2]) . "M";
	if ($t[13] && $t[13] =~ /-/) {
	  $cigar .= '1D';
	}
	print join("\t", $t[0], 1, @t[9,7,8,6], $cigar, @t[2,3]), "\n";
  }
  close($fh) if (defined($opts{s}));
}
sub soap2paf {
  my %opts;
  getopts('s', \%opts);
  die("Usage: paf_utils.pl soap2paf [-s] <in.soap>\n") if (@ARGV == 0 && -t STDIN);
  my $fh = \*STDOUT;
  open($fh, "| sort -k4,4 -k5,5n -k3,3") if (defined($opts{s}));
  while (<>) {
	my @t = split;
	my $size = length($t[1]);
	my $cigar = $size . "M";
	if ($t[9] > 100 && $t[9] <= 200) {
	  my $i = $t[9] - 100;
	  my $j = $size - $t[10];
	  $cigar = "$t[10]M$i"."D$j"."M";
	} elsif ($t[9] > 200) {
	  my $i = $t[9] - 200;
	  my $j = $size - $i - $t[10];
	  $cigar = "$t[10]M$i"."I$j"."M";
	}
	print $fh join("\t", $t[0], 1, @t[6..8], 0, $cigar, @t[1,2], 'soap', @t[4,9]);
	print "\t", join("\t", @t[10..$#t]) if ($t[9]);
	print "\n";
  }
  close($fh) if (defined($opts{s}));
}

sub rmap2paf {
  my %opts;
  getopts('s', \%opts);
  die("Usage: paf_utils.pl rmap2paf [-s] <in.rmap>\n") if (@ARGV == 0 && -t STDIN);
  my $fh = \*STDOUT;
  open($fh, "| sort -k4,4 -k5,5n -k3,3") if (defined($opts{s}));
  while (<>) {
	my @t = split;
	++$t[1]; # change to 1-based coordinate
	my $cigar = ($t[2] - $t[1] + 1)."M";
	print $fh join("\t", $t[3], 1, @t[5,0,1], 0, $cigar, '.', '.', 'rmap', $t[4]), "\n";
  }
  close($fh) if (defined $opts{s});
}

sub shrimp2paf {
  my %opts = (q=>25);
  getopts('sq:', \%opts);
  die("Usage: paf_utils.pl shrimp2paf [-s] [-q 25] <in.shrimp>\n") if (@ARGV == 0 && -t STDIN);
  my ($fh2, $fh);
  $fh = \*STDOUT;
  my $tmp = join(" ", @ARGV);
  open($fh, "| sort -k4,4 -k5,5n -k3,3") if (defined($opts{s}));
  open($fh2, qq{grep ^'>' $tmp | sed 's/[a-z_]*=//g' | sed 's/"//g' | sort -k2,2 -k5,5nr |}) || die;
  my (@s, @score, @count, @logn);
  my $q = $opts{q};
  $logn[$_] = int(3.434 * log($_) + 0.5) for (1 .. 255);
  $s[0] = '';
  while (<$fh2>) {
	if (/^>/) {
	  my @t = split;
	  if ($t[1] ne $s[0]) { # a new read
		if ($s[0]) {
		  my $s = 0;
		  if ($count[0] > 1) {
			$s = 0;
		  } elsif (@score == 1) {
			$score[0] = 3 if ($score[0] > 3);
			$s = $q * (3 - $score[0]);
		  } else {
			$s = ($score[1] - $score[0]) * $q - $logn[$count[1]];
		  }
		  $s[5] = $s;
		  print $fh join("\t", @s), "\n";
		}
		my $cigar = ($t[10]+$t[11])."M";
		$cigar .= "$t[12]D" if ($t[12]);
		$cigar .= "$t[13]I" if ($t[13]);
		@s = (@t[1,7,3,2,5], $t[4], $cigar, ".", ".", "shrimp", @t[4,11]);
		push(@s, $t[15]) if (defined $t[15]);
		@score = (); @count = ();
		$score[0] = $t[11] + $t[12] + $t[13];
		$count[0] = 1;
	  } else {
		my $s = $t[11] + $t[12] + $t[13];
		$score[@score] = $s if ($s != $score[$#score]);
		++$count[$#score];
	  }
	}
  }
  my $s = 0;
  if ($count[0] > 1) {
	$s = 0;
  } elsif (@score == 1) {
	$score[0] = 3 if ($score[0] > 3);
	$s = $q * (3 - $score[0]);
  } else {
	$s = ($score[1] - $score[0]) * $q - $logn[$count[1]];
  }
  $s[5] = $s;
  print $fh join("\t", @s), "\n";
  close($fh) if (defined $opts{s});
  close($fh2);
}

sub maq2paf {
  die("Usage: paf_utils.pl maq2paf <in.map>\n") if (@ARGV == 0);
  my ($fh, $buf, @ref);
  open($fh, "gzip -dc $ARGV[0] |") || die("[maq] Cannot decompress $ARGV[0]");
  binmode($fh);
  # read header
  $buf = '';
  read($fh, $buf, 8);
  my ($format, $n_ref) = unpack("i2", $buf);
  for (1..$n_ref) {
	read($fh, $buf, 4);
	my $len = unpack("i", $buf);
	read($fh, $buf, $len);
	($_) = unpack("Z$len", $buf);
	push(@ref, $_);
  }
  read($fh, $buf, 8);
  # read alignment
  my @ACGT = ('A', 'C', 'G', 'T');
  while (read($fh, $buf, 120) == 120) {
	my ($sq, $qse, $size, $q, $i1, $i2, $c0, $c1, $f, $aq, $sid, $x, $d, $rn)
	  = unpack("a63cC8I2iZ36", $buf);
	my $cigar = $size."M";
	if ($f >= 128 && $q != 0) {
	  $cigar = ($qse > 0)? "${q}M${qse}I".($size-$qse-$q)."M" : "${q}M".(-$qse)."D".($size-$q)."M";
	}
	$q = $aq if ($f >= 128);
	print "$rn\t1\t", ($x&1)?'-':'+', "\t$ref[$sid]\t", ($x>>1)+1, "\t$q\t$cigar";
	my @t = unpack("C$size", $sq);
	my $seq = '';
	foreach my $p (@t) {
	  if ($p == 0) { $seq .= 'N'; }
	  else { $seq .= $ACGT[$p>>6]; }
	  $p = chr(33 + ($p&0x3f));
	}
	print "\t$seq\t", @t, "\tmaq\t", join("\t", $d, $f, $qse, $aq, $i1&0xf, $i2, $c0, $c1), "\n";
  }
  close($fh);
}

sub mapview {
  my %opts;
  getopts('bq:', \%opts);
  die("Usage: paf_utils.pl mapview [-b] <in.map>\n") if (@ARGV == 0);
  my ($fh, $buf, @ref);
  open($fh, "gzip -dc $ARGV[0] |") || die("[maq] Cannot decompress $ARGV[0]");
  binmode($fh);
  # read header
  $buf = '';
  read($fh, $buf, 8);
  my ($format, $n_ref) = unpack("i2", $buf);
  for (1..$n_ref) {
	read($fh, $buf, 4);
	my $len = unpack("i", $buf);
	read($fh, $buf, $len);
	($_) = unpack("Z$len", $buf);
	push(@ref, $_);
  }
  read($fh, $buf, 8); # what is read is the number of reads, but frequently it is zero.
  # read alignment
  my $is_verbose = defined($opts{b})? 0 : 1;
  my @ACGT = ('A', 'C', 'G', 'T');
  while (read($fh, $buf, 120) == 120) {
	my ($sq, $qse, $size, $q, $i1, $i2, $c0, $c1, $f, $aq, $sid, $x, $d, $rn)
	  = unpack("a63cC8I2iZ36", $buf);
	print "$rn\t$ref[$sid]\t", ($x>>1)+1, "\t", ($x&1)?'-':'+', "\t$d\t$f\t$q\t$qse\t$aq\t",
	  $i1&0xf, "\t$i2\t$c0\t$c1\t$size";
	if ($is_verbose) { # print sequence and quality
	  my @t = unpack("C$size", $sq);
	  my $seq = '';
	  foreach my $p (@t) {
		if ($p == 0) { $seq .= 'N'; }
		else { $seq .= $ACGT[$p>>6]; }
		$p = chr(33 + ($p&0x3f));
	  }
	  print "\t$seq\t", @t;
	}
	print "\n";
  }
  close($fh);
}

sub mapview2map {
  die(qq/
Usage: paf_utils.pl mapview2map <in.list> <in.mapview> > <out.map>

Notes: This `mapview2map' command works on a stream. It begins to write
       the output before it reads to the end of the input mapview
       file. Unfortunately, this command has to write the reference
       names in the header of <out.map> and the user must provide
       correct information to make it work properly.

       The input MUST satisfy the following conditions:

       1) <in.list> is a file that consists of the names of all
          reference sequences, one name per line.

       2) <in.mapview> must be sorted according to the read position.

       3) The names in <in.list> MUST appear in the identical order of
          sequence names in <in.mapview>.
\n/) if (@ARGV == 0);

  my ($fh, @refs, %chr2id);
  my @ACGT; # array should be faster than hash
  $ACGT[$_] = 0 for (0..255);
  $ACGT[ord('A')] = $ACGT[ord('a')] = 0; $ACGT[ord('C')] = $ACGT[ord('c')] = 1;
  $ACGT[ord('G')] = $ACGT[ord('g')] = 2; $ACGT[ord('T')] = $ACGT[ord('t')] = 3;
  # read the name of the reference sequences
  $_ = shift(@ARGV);
  open($fh, $_) || die("[mapview2map] fail to open $_\n");
  while (<$fh>) {
	push(@refs, $1) if (/^(\S+)/);
  }
  close($fh);
  open($fh, "| gzip -f") || die;
  # write header
  syswrite($fh, pack("i2", -1, scalar(@refs)), 8);
  for my $i (0 .. $#refs) {
	my $chr = $refs[$i];
	$chr2id{$chr} = $i;
	my $l = length($chr) + 1;
	syswrite($fh, pack("iZ$l", $l, $chr), 4 + $l);
  }
  syswrite($fh, pack("I2", 0, 0), 8);
  # write alignment
  my ($seqid, $pos) = (0, 0);
  while (<>) {
	my @t = split;
	$t[0] = substr($t[0], length($t[0])-35, 35) if (length($t[0]) > 35);
	my $seq = '';
	my ($l, $s, $q) = @t[13..15];
	for my $p (0 .. $l-1) {
	  my $q = ord(substr $q, $p, 1)-33;
	  $q = 63 if ($q > 63);
	  $seq .= chr($ACGT[ord(substr $s, $p, 1)]<<6 | $q);
	}
	my $sid = $chr2id{$t[1]};
	die(qq/[mapview2map] undefined sequence name '$t[1]'\n/) unless (defined $sid);
	$t[1] = $sid;
	die(qq/[mapview2map] inconsistent order of reference sequences\n/) if ($sid < $seqid);
	if ($sid > $seqid) {
	  $seqid = $sid;
	  $pos = 0;
	}
	die(qq/[mapview2map] unsorted mapview input\n/) if ($t[2] < $pos);
	$t[2] = ($t[2]-1)<<1 | (($t[3] eq '+')? 0:1);
	syswrite($fh, pack("a63cC8I2iZ36", $seq, @t[7,13,6,9..12,5,8,1,2,4,0]), 120);
  }
  close($fh);
}