Mercurial > repos > mcharles > rapsosnp

comparison rapsodyn/filtersam_mapped_and_unique.pl @ 7:3f7b0788a1c4 draft

Find changesets by keywords (author, files, the commit message), revision number or hash, or revset expression.
Uploaded
author mcharles
date Tue, 07 Oct 2014 10:34:34 -0400
parents 442a7c88b886
children 0a6c1cfe4dc8
comparison
equal deleted inserted replaced
6:1776b8ddd87e 7:3f7b0788a1c4
1 #!/usr/bin/perl 1 #!/usr/bin/perl
2 #V1.0.1 added log, option parameters
2 use strict; 3 use strict;
3 use warnings; 4 use warnings;
5 use Getopt::Long;
4 6
5 open(IN, $ARGV[0]) or die ("Can't open $ARGV[0]\n"); 7 my $input_sam_file;
8 my $output_sam_file;
9 my $log_file;
10
11 my %bitscore_all;
12 my %bitscore_selected;
13
14 GetOptions (
15 "input_sam_file=s" => \$input_sam_file,
16 "output_sam_file=s" => \$output_sam_file,
17 "log_file=s" => \$log_file
18 ) or die("Error in command line arguments\n");
19
20 open(IN, $input_sam_file) or die ("Can't open $input_sam_file\n");
6 while (my $line=<IN>){ 21 while (my $line=<IN>){
7 if ($line =~ /^\@/){ 22 if (($line =~ /^\@SQ/)||($line =~ /^\@PG/)){
8 #Header conservation 23 #Header conservation
9 print $line; 24 print $line;
10 } 25 }
11 else { 26 else {
12 #Optionnal flag verification 27 #Optionnal flag verification
28 my @fields_all = split (/\s+/,$line);
29 my $bit = $fields_all[1];
30 if ($bitscore_all{$bit}){
31 $bitscore_all{$bit}++;
32 }
33 else {
34 $bitscore_all{$bit}=1;
35 }
13 if (($line =~ /XT\:A\:U/)&&($line =~ /X0\:i\:1/)&&($line =~ /X1\:i\:0\s/)){ 36 if (($line =~ /XT\:A\:U/)&&($line =~ /X0\:i\:1/)&&($line =~ /X1\:i\:0\s/)){
14 my @fields = split (/\s+/,$line); 37 my @fields_selected = split (/\s+/,$line);
15 if (($fields[1]==83)||($fields[1]==163)||($fields[1]==147)||($fields[1]==99)){ 38 if (($fields_selected[1]==83)||($fields_selected[1]==163)||($fields_selected[1]==147)||($fields_selected[1]==99)){
16 print $line; 39 print $line;
40 my $bit = $fields_selected[1];
41 if ($bitscore_selected{$bit}){
42 $bitscore_selected{$bit}++;
43 }
44 else {
45 $bitscore_selected{$bit}=1;
46 }
17 } 47 }
18 } 48 }
19 } 49 }
20 } 50 }
21 51
52 close (IN);
22 53
23 close (IN); 54 open (LF,">$log_file") or die("Can't open $log_file\n");
55 print LF "\n####\t Sam filtering \n";
56 print LF "## Before filtering\n";
57 print LF "bitscore\t:\t";
58 foreach my $key (sort {$bitscore_all{$b} <=> $bitscore_all{$a}} keys %bitscore_all) {
59 print LF $key,"\t*\t";
60 }
61 print LF "\n number \t:\t";
62 foreach my $key (sort {$bitscore_all{$b} <=> $bitscore_all{$a}} keys %bitscore_all) {
63 print LF $bitscore_all{$key},"\t*\t";
64 }
65 print LF "\n";
66 print LF "## After filtering\n";
67 print LF "bitscore\t:\t";
68 foreach my $key (sort {$bitscore_selected{$b} <=> $bitscore_selected{$a}} keys %bitscore_selected) {
69 print LF $key,"\t*\t";
70 }
71 print LF "\n number \t:\t";
72 foreach my $key (sort {$bitscore_selected{$b} <=> $bitscore_selected{$a}} keys %bitscore_selected) {
73 print LF $bitscore_selected{$key},"\t*\t";
74 }
75 print LF "\n";
76 close (LF);
77
78
79
80
81
82
83