Mercurial > repos > siyuan > prada

comparison pyPRADA_1.2/tools/samtools-0.1.16/misc/zoom2sam.pl @ 0:acc2ca1a3ba4

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author siyuan
date Thu, 20 Feb 2014 00:44:58 -0500
parents
children
comparison
equal deleted inserted replaced
-1:000000000000 0:acc2ca1a3ba4
1 #!/usr/bin/perl -w
2
3 # Contact: lh3
4 # Version: 0.1.0
5
6 use strict;
7 use warnings;
8 use Getopt::Std;
9
10 &zoom2sam;
11 exit;
12
13 sub mating {
14 my ($s1, $s2) = @_;
15 my $isize = 0;
16 if ($s1->[2] ne '*' && $s1->[2] eq $s2->[2]) { # then calculate $isize
17 my $x1 = ($s1->[1] & 0x10)? $s1->[3] + length($s1->[9]) : $s1->[3];
18 my $x2 = ($s2->[1] & 0x10)? $s2->[3] + length($s2->[9]) : $s2->[3];
19 $isize = $x2 - $x1;
20 }
21 # update mate coordinate
22 if ($s2->[2] ne '*') {
23 @$s1[6..8] = (($s2->[2] eq $s1->[2])? "=" : $s2->[2], $s2->[3], $isize);
24 $s1->[1] |= 0x20 if ($s2->[1] & 0x10);
25 } else {
26 $s1->[1] |= 0x8;
27 }
28 if ($s1->[2] ne '*') {
29 @$s2[6..8] = (($s1->[2] eq $s2->[2])? "=" : $s1->[2], $s1->[3], -$isize);
30 $s2->[1] |= 0x20 if ($s1->[1] & 0x10);
31 } else {
32 $s2->[1] |= 0x8;
33 }
34 }
35
36 sub zoom2sam {
37 my %opts = ();
38 getopts("p", \%opts);
39 die("Usage: zoom2sam.pl [-p] <readLen> <aln.zoom>
40 Warnings: This script only supports the default Illumina outputs.\n") if (@ARGV < 2);
41 my $is_paired = defined($opts{p});
42 my $len = shift(@ARGV);
43 # core loop
44 my @s1 = ();
45 my @s2 = ();
46 my ($s_last, $s_curr) = (\@s1, \@s2);
47 while (<>) {
48 &zoom2sam_aux($_, $s_curr, $is_paired, $len);
49 if (@$s_last != 0 && $s_last->[0] eq $s_curr->[0]) {
50 &mating($s_last, $s_curr);
51 print join("\t", @$s_last), "\n";
52 print join("\t", @$s_curr), "\n";
53 @$s_last = (); @$s_curr = ();
54 } else {
55 print join("\t", @$s_last), "\n" if (@$s_last != 0);
56 my $s = $s_last; $s_last = $s_curr; $s_curr = $s;
57 }
58 }
59 print join("\t", @$s_last), "\n" if (@$s_last != 0);
60 }
61
62 sub zoom2sam_aux {
63 my ($line, $s, $is_paired, $len) = @_;
64 chomp($line);
65 my @t = split("\t", $line);
66 @$s = ();
67 # read name
68 $s->[0] = $t[0];
69 # initial flag (will be updated later)
70 $s->[1] = 0;
71 $s->[1] |= 1 | 1<<6 if ($s->[0] =~ /_F$/);
72 $s->[1] |= 1 | 1<<7 if ($s->[0] =~ /_R$/);
73 $s->[1] |= 2 if ($is_paired);
74 # read & quality
75 $s->[9] = "*"; $s->[10] = "*";
76 # cigar
77 $s->[5] = $len . "M";
78 # coor
79 my @s = split(/\s+/, $t[1]);
80 $s->[2] = $s[0];
81 $t[1] =~ /:(\d+)$/;
82 $s->[3] = $1 + 1;
83 if ($s->[0] =~ /_[FR]$/) {
84 my $u = ($s->[0] =~ /_F$/)? 1 : 0;
85 my $w = ($t[2] eq '+')? 1 : 0;
86 $s->[1] |= 0x10 if ($u ^ $w);
87 $s->[0] =~ s/_[FR]$//;
88 } else {
89 $s->[1] |= 0x10 if ($t[2] eq '-');
90 }
91 # mapQ
92 $s->[4] = 30;
93 # mate coordinate
94 $s->[6] = '*'; $s->[7] = $s->[8] = 0;
95 # aux
96 push(@$s, "NM:i:$t[3]");
97 }