Python+Perl联合解决基因组染色体拼接错误问题

原创

已于 2024-03-13 16:40:09 修改 · 1.1k 阅读

28 ·

CC 4.0 BY-SA版权

文章标签：

#生物信息学 #python #perl #linux

于 2023-12-19 21:32:53 首次发布

文章讲述了如何使用Perl脚本根据染色体上的gap打断序列，然后手动调整基因注释顺序，最后使用Python脚本来合并修复好的染色体序列和基因注释，确保基因坐标正确。

基因组拼接完成甚至注释完之后，由于有个染色体因为Hi-C信号不是很丰富，导致着丝粒前的序列和着丝粒后的序列装反了，直到用jcvi开始共线性分析才发现，所以想办法直接对现有的两个文件，该染色体序列和该染色体对应的基因结构注释进行一个顺序变换。不过本操作有两个前提： 1、染色体仍是存在gap的；2、跨gap的基因应当被去除。 根据现在的情况，我要将下图右边的块移动到最左边。
正好反了，把后面一大段移到最前面就能解决
由于跨gap的基因应该被去除，所以我只要将这条染色体根据gap再打断成contig（实际上不是传统意义上的contig，这里专指被N分开的若干段序列），重排contig之后就能获得新染色体序列，而每个contig上的基因相对于contig是不会变化的，只需要计算contig移动之后基因的新坐标就可以获得新染色体基因注释。

1、使用Perl脚本根据gap打断染色体

#!/usr/bin/perl
use strict;
use warnings;

my $author = "Alter X";
my $usage = "$0 <input_chr.fa> <output_broken_contig.list>";
die($usage) unless (@ARGV == 2);

my $fasta_input_file = $ARGV[0];
my $broken_contig_file = $ARGV[1];

my $seq_id;
my $seq = "";
my @seq_data;

open(my $fasta_input,"<",$fasta_input_file) or die "$!";
while(<$fasta_input>){
   
   
    chomp;
    $seq_id = $1 and next if (m/^>(.*)/);
    $seq .= $_;
}
close($fasta_input);

#最精华的一条代码，实现根据gap打断染色体，同时保留gap的序列，分别从存入数组seq_data中
@seq_data = $seq =~ /([ACGT]+|N+)/g;

my $num = 1;
my $length;
my $start = 0;
my $end = 0;

open(my $broken_contig,">",$broken_contig_file);
for my $sub_seq (@seq_data){
   
   
    $start = $end + 1;
    $length = length($sub_seq);
    $end = $start + $length - 1;
    if ($sub_seq =~ /