phylop 中遇到的maf sort 与 split 问题

在运行phylop时遇到"ERROR maf_read_subset: invalid idx_offset -2" 的报错, 作者是这样解释的

image.png

1、 maf文件的第一个物种在每一个比对块中必须是一样的，作为参考基因组
2、 maf文件比对块的顺序必须是被排序的
3、 maf文件的参考基因组文件，应该只包含一条染色体

上述描述翻译成人话就是，maf文件需要按照参考基因组的染色体分割以及按照坐标由低到高排序，但是他分享的这些maf_sort，maf-sort， mafSplit.py 都不太好用，花半天时间下载可能还有一堆报错，所以还不如自己写一个简单的split 与 sort 的脚本。

参考maf格式如下

e12a1c1d1f3c41c7b4b79c92cc875a6.png

自动在maf输入文件下生成分割并排序好的maf，这个脚本可能会比较占内存，maf过大的或者服务器内存比较小的慎用，脚本如下

#!/usr/bin/envs perl -w

use strict;

my $usage=<<USAGE;

perl $0 maf

script will create each splited maf by chromosomes of reference at the path of maf input.

USAGE

die "$usage\n" unless @ARGV == 1;

open MF, "$ARGV[0]" or die "";
$/="\na\n";
my ($header, %hash);
while (<MF>)
{
        chomp;
        if (m/^#/)
        {
                $header="$_\n";
        }else{
                my @line=split/\s+/, $_;
                $hash{$line[1]}{$line[2]}=$_;
        }
}

close MF;

$/="\n";

foreach my $chr (sort keys %hash){
        next if ($chr =~ m/ChrM/ || $chr =~m/ChrC/);
        #next unless $chr =~ m/A_tha_pep/;
        my $out=split/\./, $chr;
        open OUT, ">$chr.maf" or die "";
        print OUT "$header\n";
        foreach my $pos (sort {$a<=>$b} keys %{$hash{$chr}}){
                print OUT "a\n$hash{$chr}{$pos}";
        }
        close OUT;
        print "chr has been done! \n ";
}