[笔记] sam文件挑取uniq reads

本来呢，samtools本来就可以一句话就能提取，只是没得网络没得下载samtools，只能保留自传脚本。

输入文件为raw.sam

输出为提取uniq reads后的uniq.sam（仍然保留着sam格式）

终端里会打印出简单的报告

use strict;
use warnings;

my ($in, $out) = @ARGV; 
die <DATA> unless (@ARGV == 2);
open IN, "$in";
open OUT, ">$out";

my ($sum, $uniq, $xnd, $other) = (0,0,0,0); 
my $header;
while (<IN>){

    my ($a, $b) = (0, 0);
    chomp;
    if ($_ =~ /^\@/){
        print OUT "$_\n";
        $header ++;
    }
    else {$sum ++};
    $a = 1 if (/AS:i/);
    $b = 1 if (/XS:i/);
    if ($a ==1 ){
        if ($b == 0){
            $uniq ++;
            print OUT "$_\n";
        }
        if ($b == 1){
            $xnd ++;
        }
    }
    elsif ($a == 0){
        $other ++;
        #header number will counted in this list
    }
}
print "input sam is : $in\n";
print "output uniq sam is : $out\n";
print "all reads num: $sum\n";
print "uniq reads num: $uniq\n";
print "2nd reads num: $xnd\n";
print "error reads num: " . ($other-$header) . "\n";
my $ratio = $uniq/$sum;
printf "uniq percentage: %.1f\n", ($ratio*100); 

__DATA__
usage: ./sam_uniq_reads.pl <sam_input.sam> <uniq_reads_output.sam>