sam2gff3

最新推荐文章于 2025-04-24 15:54:36 发布

原创最新推荐文章于 2025-04-24 15:54:36 发布 · 843 阅读

1 ·

CC 4.0 BY-SA版权

文章标签：

#sam #gff3 #trans

文本转化同时被 2 个专栏收录

1 篇文章

订阅专栏

perl

1 篇文章

订阅专栏

博客讨论了在将sam格式的序列比对文件转换为gff3格式时可能出现的重复问题。由于一条序列可能比对到基因组多个位置，转换后会导致gff3文件中的重复。通常，这不影响大部分研究，但可以通过筛选sam文件来解决此问题。

摘要生成于 C知道，由 DeepSeek-R1 满血版支持，前往体验 >

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

#!/usr/bin/perl -w
use 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}};
    }
}