annotate PsiCLASS-1.0.2/samtools-0.1.19/misc/blast2sam.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 use strict;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
4 use warnings;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
5 use Getopt::Std;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
6
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
7 &blast2sam;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
8
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
9 sub blast2sam {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
10 my %opts = ();
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
11 getopts('s', \%opts);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
12 die("Usage: blast2sam.pl <in.blastn>\n") if (-t STDIN && @ARGV == 0);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
13 my ($qlen, $slen, $q, $s, $qbeg, $qend, @sam, @cigar, @cmaux, $show_seq);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
14 $show_seq = defined($opts{s});
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
15 @sam = (); @sam[0,4,6..8,10] = ('', 255, '*', 0, 0, '*');
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
16 while (<>) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
17 if (@cigar && (/^Query=/ || /Score =.*bits.*Expect/)) { # print
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
18 &blast_print_sam(\@sam, \@cigar, \@cmaux, $qlen - $qend);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
19 @cigar = ();
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
20 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
21 if (/^Query= (\S+)/) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
22 $sam[0] = $1;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
23 } elsif (/\((\S+)\s+letters\)/) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
24 $qlen = $1; $qlen =~ s/,//g;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
25 } elsif (/^>(\S+)/) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
26 $sam[2] = $1;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
27 } elsif (/Length = (\d+)/) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
28 $slen = $1;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
29 } elsif (/Score =\s+(\S+) bits.+Expect(\(\d+\))? = (\S+)/) { # the start of an alignment block
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
30 my ($as, $ev) = (int($1 + .499), $3);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
31 $ev = "1$ev" if ($ev =~ /^e/);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
32 @sam[1,3,9,11,12] = (0, 0, '', "AS:i:$as", "EV:Z:$ev");
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
33 @cigar = (); $qbeg = 0;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
34 @cmaux = (0, 0, 0, '');
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
35 } elsif (/Strand = (\S+) \/ (\S+)/) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
36 $sam[1] |= 0x10 if ($2 eq 'Minus');
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
37 } elsif (/Query\:\s(\d+)\s*(\S+)\s(\d+)/) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
38 $q = $2;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
39 unless ($qbeg) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
40 $qbeg = $1;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
41 push(@cigar, ($1-1) . "H") if ($1 > 1);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
42 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
43 $qend = $3;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
44 if ($show_seq) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
45 my $x = $q;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
46 $x =~ s/-//g; $sam[9] .= $x;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
47 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
48 } elsif (/Sbjct\:\s(\d+)\s*(\S+)\s(\d+)/) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
49 $s = $2;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
50 if ($sam[1] & 0x10) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
51 $sam[3] = $3;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
52 } else {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
53 $sam[3] = $1 unless ($sam[3]);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
54 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
55 &aln2cm(\@cigar, \$q, \$s, \@cmaux);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
56 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
57 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
58 &blast_print_sam(\@sam, \@cigar, \@cmaux, $qlen - $qend);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
59 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
60
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
61 sub blast_print_sam {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
62 my ($sam, $cigar, $cmaux, $qrest) = @_;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
63 push(@$cigar, $cmaux->[1] . substr("MDI", $cmaux->[0], 1));
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
64 push(@$cigar, $qrest . 'H') if ($qrest);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
65 if ($sam->[1] & 0x10) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
66 @$cigar = reverse(@$cigar);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
67 $sam->[9] = reverse($sam->[9]);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
68 $sam->[9] =~ tr/atgcrymkswATGCRYMKSW/tacgyrkmswTACGYRKMSW/;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
69 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
70 $sam->[9] = '*' if (!$sam->[9]);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
71 $sam->[5] = join('', @$cigar);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
72 print join("\t", @$sam), "\n";
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
73 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
74
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
75 sub aln2cm {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
76 my ($cigar, $q, $s, $cmaux) = @_;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
77 my $l = length($$q);
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
78 for (my $i = 0; $i < $l; ++$i) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
79 my $op;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
80 # set $op
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
81 if (substr($$q, $i, 1) eq '-') { $op = 2; }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
82 elsif (substr($$s, $i, 1) eq '-') { $op = 1; }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
83 else { $op = 0; }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
84 # for CIGAR
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
85 if ($cmaux->[0] == $op) {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
86 ++$cmaux->[1];
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
87 } else {
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
88 push(@$cigar, $cmaux->[1] . substr("MDI", $cmaux->[0], 1));
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
89 $cmaux->[0] = $op; $cmaux->[1] = 1;
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
90 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
91 }
903fc43d6227 Uploaded
lsong10
parents:
diff changeset
92 }