sam2gff3

需要注意的是，有时候一条序列可以比对到基因组多个地方，转换成gff3格式之后就会出现重复。这就要评估你的研究目的了。一般情况不会影响。当然也可以将sam文件中这种情况进行筛选。

#!/usr/bin/perl -wuse strict;open IN,"2.sam"; #输入文件open OUT,">output.gff3";#输出文件my %gtf;my %pos;while (<IN>){    chomp;    next if /^@/;    my @array=split/\t/;    my @num=split/\D+/,$array[5];    my @base=split/\d+/,$array[5];    shift @base;    my @pos;    my $start=$array[3];    my $end=$array[3];    for my $i(0..$#num){    if ($base[$i] eq "H" or $base[$i] eq "S"){        if ($start <=> $end){        push @pos,($start,$end-1);        }        $end+=$num[$i];        $start=$end;    }elsif ($base[$i] eq "N"){        if ($start <=> $end){        push @pos,($start,$end-1);        }        $end+=$num[$i];        $start=$end;    }elsif($base[$i] eq "M"){        $end+=$num[$i];    }elsif($base[$i] eq "D"){        $end+=$num[$i];    }elsif($base[$i] eq "I"){        next;    }else{        print "line$.,$base[$i]\n";    }    }    if ($start <=> $end){    push @pos,($start,$end-1);    }    my $direct;    if ($array[1]==16){    $direct="-"    }else{    $direct="+";    }my $nam= "$array[2]\tSAMtrans\tgene\t$pos[0]\t$pos[-1]\t.\t$direct\t.\tID=$array[0];Name=".(split(/\./,$array[0]))[0]."\n";    push @{$gtf{$array[2]}{$array[0]}},$nam;    for (my $i=0;$i<=$#pos;$i+=2){        my $mlj="$array[2]\tSAMtrans\texon\t$pos[$i]\t$pos[$i+1]\t.\t$direct\t.\tParent=$array[0];Name=".(split(/\./,$array[0]))[0]."\n";    push @{$gtf{$array[2]}{$array[0]}},$mlj}    $pos{$array[2]}{$array[0]}=$pos[0];}foreach my $chr(sort keys %pos){    foreach my $gene(sort {$pos{$chr}{$a}<=>$pos{$chr}{$b}} keys %{$pos{$chr}}){    print OUT "$_" foreach @{$gtf{$chr}{$gene}};    }}