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

Uploaded

#!/usr/bin/perl
#V1.0.1 added log, option parameters
use strict;
use warnings;
use Getopt::Long;

my $input_sam_file;
my $output_sam_file;
my $log_file;

my %bitscore_all;
my %bitscore_selected;

GetOptions (
"input_sam_file=s" => \$input_sam_file,
"output_sam_file=s" => \$output_sam_file,
"log_file=s" => \$log_file
) or die("Error in command line arguments\n");

open(IN, $input_sam_file) or die ("Can't open $input_sam_file\n");
while (my $line=<IN>){
	if (($line =~ /^\@SQ/)||($line =~ /^\@PG/)){
		#Header conservation
		print $line;
	}
	else {
		#Optionnal flag verification
		my @fields_all = split (/\s+/,$line);
		my $bit = $fields_all[1];
		if ($bitscore_all{$bit}){
			$bitscore_all{$bit}++;
		}
		else {
			$bitscore_all{$bit}=1;
		}
		if (($line =~ /XT\:A\:U/)&&($line =~ /X0\:i\:1/)&&($line =~ /X1\:i\:0\s/)){
			my @fields_selected = split (/\s+/,$line);
			if (($fields_selected[1]==83)||($fields_selected[1]==163)||($fields_selected[1]==147)||($fields_selected[1]==99)){
				print $line;
				my $bit = $fields_selected[1];
				if ($bitscore_selected{$bit}){
					$bitscore_selected{$bit}++;
				}
				else {
					$bitscore_selected{$bit}=1;
				}
			}
		}
	}
}

close (IN);

open (LF,">$log_file") or die("Can't open $log_file\n");
print LF "\n####\t Sam filtering \n";
print LF "## Before filtering\n";
print LF "bitscore\t:\t";
foreach my $key (sort {$bitscore_all{$b} <=> $bitscore_all{$a}} keys %bitscore_all) {
	print LF $key,"\t*\t";
}
print LF "\n number \t:\t";
foreach my $key (sort {$bitscore_all{$b} <=> $bitscore_all{$a}} keys %bitscore_all) {
	print LF $bitscore_all{$key},"\t*\t";
}
print LF "\n";
print LF "## After filtering\n";
print LF "bitscore\t:\t";
foreach my $key (sort {$bitscore_selected{$b} <=> $bitscore_selected{$a}} keys %bitscore_selected) {
	print LF $key,"\t*\t";
}
print LF "\n number \t:\t";
foreach my $key (sort {$bitscore_selected{$b} <=> $bitscore_selected{$a}} keys %bitscore_selected) {
	print LF $bitscore_selected{$key},"\t*\t";
}
print LF "\n";
close (LF);