annotate PsiCLASS-1.0.2/samtools-0.1.19/misc/psl2sam.pl @ 0:903fc43d6227 draft default tip

Uploaded
author lsong10
date Fri, 26 Mar 2021 16:52:45 +0000
parents
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
1 #!/usr/bin/perl -w
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
2
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
3 # Author: lh3
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
4
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
5 # This script calculates a score using the BLAST scoring
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
6 # system. However, I am not sure how to count gap opens and gap
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
7 # extensions. It seems to me that column 5-8 are not what I am
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
8 # after. This script counts gaps from the last three columns. It does
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
9 # not generate reference skip (N) in the CIGAR as it is not easy to
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
10 # directly tell which gaps correspond to introns.
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
11
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
12 use strict;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
13 use warnings;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
14 use Getopt::Std;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
15
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
16 my %opts = (a=>1, b=>3, q=>5, r=>2);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
17 getopts('a:b:q:r:', \%opts);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
18 die("Usage: psl2sam.pl [-a $opts{a}] [-b $opts{b}] [-q $opts{q}] [-r $opts{r}] <in.psl>\n") if (@ARGV == 0 && -t STDIN);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
19
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
20 my @stack;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
21 my $last = '';
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
22 my ($a, $b, $q, $r) = ($opts{a}, $opts{b}, $opts{q}, $opts{r});
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
23 while (<>) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
24 next unless (/^\d/);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
25 my @t = split;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
26 my @s;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
27 my $cigar = '';
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
28 if ($t[8] eq '-') {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
29 my $tmp = $t[11];
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
30 $t[11] = $t[10] - $t[12];
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
31 $t[12] = $t[10] - $tmp;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
32 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
33 @s[0..4] = ($t[9], (($t[8] eq '+')? 0 : 16), $t[13], $t[15]+1, 0);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
34 @s[6..10] = ('*', 0, 0, '*', '*');
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
35 $cigar .= $t[11].'H' if ($t[11]); # 5'-end clipping
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
36 my @x = split(',', $t[18]);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
37 my @y = split(',', $t[19]);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
38 my @z = split(',', $t[20]);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
39 my ($y0, $z0) = ($y[0], $z[0]);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
40 my ($gap_open, $gap_ext) = (0, 0, 0);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
41 for (1 .. $t[17]-1) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
42 my $ly = $y[$_] - $y[$_-1] - $x[$_-1];
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
43 my $lz = $z[$_] - $z[$_-1] - $x[$_-1];
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
44 if ($ly < $lz) { # del: the reference gap is longer
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
45 ++$gap_open;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
46 $gap_ext += $lz - $ly;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
47 $cigar .= ($y[$_] - $y0) . 'M';
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
48 $cigar .= ($lz - $ly) . 'D';
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
49 ($y0, $z0) = ($y[$_], $z[$_]);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
50 } elsif ($lz < $ly) { # ins: the query gap is longer
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
51 ++$gap_open;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
52 $gap_ext += $ly - $lz;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
53 $cigar .= ($z[$_] - $z0) . 'M';
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
54 $cigar .= ($ly - $lz) . 'I';
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
55 ($y0, $z0) = ($y[$_], $z[$_]);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
56 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
57 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
58 $cigar .= ($t[12] - $y0) . 'M';
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
59 $cigar .= ($t[10] - $t[12]).'H' if ($t[10] != $t[12]); # 3'-end clipping
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
60 $s[5] = $cigar;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
61 my $score = $a * $t[0] - $b * $t[1] - $q * $gap_open - $r * $gap_ext;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
62 $score = 0 if ($score < 0);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
63 $s[11] = "AS:i:$score";
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
64 print join("\t", @s), "\n";
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
65 }