annotate pyPRADA_1.2/tools/samtools-0.1.16/misc/interpolate_sam.pl @ 3:f17965495ec9 draft default tip

Uploaded
author siyuan
date Tue, 11 Mar 2014 12:14:01 -0400
parents acc2ca1a3ba4
children
Ignore whitespace changes - Everywhere: Within whitespace: At end of lines:
rev   line source
0
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
1 #!/usr/bin/perl
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
2 use strict;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
3
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
4 ###Builds interpolated pileup from SAM file
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
5 ##@description counts bases between paired ends and piles up single end reads.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
6 ##@output, uses a #header for the RNAME and then the number of reads per base
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
7 ##@author sm8@sanger.ac.uk, Stephen B. Montgomery
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
8
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
9 ##@caveats
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
10 ##Requires RNAME to have format as per example
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
11 ## chromosome:NCBI36:18:1:76117153:1
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
12 ## supercontig::NT_113883:1:137703:1
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
13 ## clone::AC138827.3:1:149397:1
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
14 ##Expects simple CIGAR characters, M, I and D
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
15 ##Expects SAM file to be sorted.
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
16 ##Expects 0x0010 to mark second read in PE file (as has been the observed case from MAQ output) (important for line 77)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
17
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
18 ##Verify and read in SAM file
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
19 my $sam_file = $ARGV[0];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
20 if(!defined($sam_file)) { die("No sam file defined on arg 1"); }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
21 unless(-f $sam_file) { die("Sam file does not exist: $sam_file"); }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
22 open(SAM, $sam_file) || die("Cannot open sam file");
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
23
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
24 ##Globals
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
25 my $current_location = ""; ##Current RNAME being processed
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
26 my $current_size = 0; ##Size of sequence region being processed
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
27 my $current_position = 1; ##Current base being processed
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
28 my $open = 0; ##Number of open reads (PE reads that have not been closed)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
29 my %close = (); ##Hash of closing positions, when the current_position gets to this position it subtracts the
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
30 ##contained value from those open and deletes the indexed position from the hash
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
31
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
32 while (my $line = <SAM>) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
33 my @tokens = split /\t/, $line;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
34
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
35 if ($current_location ne $tokens[2]) { ##Start a new sequence region
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
36 for (my $i = $current_position; $i <= $current_size; $i++) { ##Close the previous sequence region
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
37 if (defined($close{$i})) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
38 $open = $open - $close{$i};
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
39 delete $close{$i};
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
40 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
41 print $open . "\n";
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
42 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
43 if ($current_location ne "") {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
44 print "\n";
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
45 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
46
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
47 ##Initiate a new sequence region
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
48 my @location_tokens = split /:/, $tokens[2];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
49 $current_position = 1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
50 $current_location = $tokens[2];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
51 $current_size = $location_tokens[4];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
52 $open = 0;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
53 %close = ();
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
54 print "#" . $tokens[2] . "\n";
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
55
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
56 ##Print pileup to just before the first read (will be 0)
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
57 for (my $current_position = 1; $current_position < $tokens[3]; $current_position++) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
58 print $open . "\n";
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
59 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
60 $current_position = $tokens[3];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
61
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
62 } else { ##Sequence region already open
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
63 if ($tokens[3] > $current_position) { ##If the new read's position is greater than the current position
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
64 ##cycle through to catch up to the current position
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
65 for (my $i = $current_position; $i < $tokens[3]; $i++) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
66 if (defined($close{$i})) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
67 $open = $open - $close{$i};
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
68 delete $close{$i};
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
69 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
70 print $open . "\n";
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
71 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
72 $current_position = $tokens[3];
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
73 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
74 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
75 $open++; ##Increment the number of open reads
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
76
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
77 if (($tokens[1] & 0x0080 || $tokens[1] & 0x0040) && $tokens[1] & 0x0010 && $tokens[1] & 0x0002) { ##if second read of mate pair, add close condition
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
78 $open--;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
79 my $parsed_cig = &parseCigar($tokens[5]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
80 my $seq_region_end = $tokens[3] + $parsed_cig->{'M'} + $parsed_cig->{'D'} - 1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
81 if (!defined($close{$seq_region_end + 1})) { $close{$seq_region_end + 1} = 0; }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
82 $close{$seq_region_end + 1} = $close{$seq_region_end + 1} + 1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
83 } elsif (!($tokens[1] & 0x0001) || !($tokens[1] & 0x0002)) { ##if unpaired, add close condition
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
84 my $parsed_cig = &parseCigar($tokens[5]);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
85 my $seq_region_end = $tokens[3] + $parsed_cig->{'M'} + $parsed_cig->{'D'} - 1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
86 if (!defined($close{$seq_region_end + 1})) { $close{$seq_region_end + 1} = 0; }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
87 $close{$seq_region_end + 1} = $close{$seq_region_end + 1} + 1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
88 } else {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
89 #do nothing
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
90 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
91 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
92 for (my $i = $current_position; $i <= $current_size; $i++) { ##Finish up the last sequence region
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
93 if (defined($close{$i})) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
94 $open = $open - $close{$i};
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
95 delete $close{$i};
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
96 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
97 print $open . "\n";
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
98 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
99 print "\n";
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
100 close(SAM);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
101 exit(0);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
102
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
103 ##reads and tokenizes simple cigarline
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
104 sub parseCigar() {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
105 my $cigar_line = shift;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
106 $cigar_line =~ s/([0-9]*[A-Z]{1})/$1\t/g;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
107 my @cigar_tokens = split /\t/, $cigar_line;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
108 my %parsed = ('M' => 0,
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
109 'I' => 0,
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
110 'D' => 0);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
111 my @events = ();
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
112 for(my $i = 0; $i < scalar(@cigar_tokens); $i++) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
113 if ($cigar_tokens[$i] =~ /([0-9]+)([A-Z]{1})/g) {
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
114 if (!defined($parsed{$2})) { $parsed{$2} = 0; }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
115 my $nt = $2;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
116 if ($nt ne "M" && $nt ne "D" && $nt ne "I") { $nt = "M"; }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
117 $parsed{$nt} += $1;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
118 my %event_el = ("t" => $nt,
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
119 "n" => $1);
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
120 push @events, \%event_el;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
121 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
122 }
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
123 $parsed{'events'} = \@events;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
124 return \%parsed;
acc2ca1a3ba4 Uploaded
siyuan
parents:
diff changeset
125 }