By liyupeng, 29 December, 2025
Forums
Researches

seqkit倒转染色体方向

介绍

seqkit 能够通过反向互补操作来倒转染色体的方向。具体来说，它首先对DNA序列中的每个碱基进行互补替换：A变为T，T变为A，C变为G，G变为C。然后，它将整个序列的顺序倒置，即将末端的碱基移到起始位置，形成反向互补的序列。这一过程常用于基因组学中的序列方向修正，尤其是在处理基因注释时，帮助确保正链和反链的正确表示和分析。通过这种方式，seqkit使得原本正向的DNA序列转换为反向序列，以适应不同的分析需求。

使用方法

  1. 先找一个方向调整好的基因组，和目标基因组进行全基因组比对
    1. 这里用的基因组都是杜鹃的基因组，染色体数量也一样，方向的话Rsimsii这个是调整过的，所以只要通过全基因组比对两个基因组中染色体的连线情况，就能够判断目标基因组染色体（ctg_chr_scaffolds.fa）的方向情况
nohup /home/soft/NGenomeSyn/bin/GetTwoGenomeSyn.pl \
    -InGenomeA ../Rsimsii.chr.fa \
    -InGenomeB ctg_chr_scaffolds.fa \
    -OutPrefix RS_vs_ctg \
    -BinDir /home/soft/mummer-4.0.1 \
    -MinAlnLen 10000 &
  1. 全基因组比对结果图
  2. 老实说，全基因组比对大概是无法判断方向情况，个别倒是能够看得出来，但是连染色体编号也无法知道的情况下，那就牡蛎了
  3. 但是，GetTwoGenomeSyn.pl生成的图片虽然用不上，但它基于比对的数据结果（.link）的这个文件，只要加以解读就能够知道哪个染色体对应到哪个染色体上

 

  1. 统计每一对染色体的比对结果的数量

awk '{key=$1"\t"$4; count[key]++} END{for (k in count) print k"\t"count[k]}' RS_vs_ctg.link \
| sort -k1,1 -k2,2 > chr_scaffold_link.count.txt
  1. 筛选最有染色体比对
    1. 两个基因组各有13条染色体，比对的组合一共有91中。不过，吾辈只需要最优的那个比对。
    2. 通过计算，保留比对数量最多的，这就是最后的13个比对结果统计，也没有重复

awk '
{
  chr=$1; scaffold=$2; cnt=$3
  if (!(chr in max) || cnt > max[chr]) {
    max[chr]=cnt
    best[chr]=scaffold
  }
}
END{
  for (chr in max)
    print chr "\t" best[chr] "\t" max[chr]
}' chr_scaffold_link.count.txt \
| sort -k1,1 > chr_best_scaffold.txt
  1. 拆分染色体，按照对应的染色体重新进行比对，查看方向
    1. 拆分Rsimsii基因组
awk '
/^>/ {
  if (out) close(out)
  name=substr($0,2)
  out=name".fa"
}
{ print > out }
' Rsimsii.chr.fa
  1. 拆分ctg基因组
awk '
/^>/ {
  if (out) close(out)
  name=substr($0,2)
  out=name".fa"
}
{
  print > out
}
' ctg_chr_scaffolds.fa
  1. 单个染色体进行比对
    1. 相较于一整个基因组的比对，单个染色体之间的比对能够看得更加清楚一些
nohup /home/soft/NGenomeSyn/bin/GetTwoGenomeSyn.pl \
    -InGenomeA chr01.fa \
    -InGenomeB scaffold_1.fa \
    -OutPrefix 01 \
    -BinDir /home/soft/mummer-4.0.1 \
    -MinAlnLen 10000 &
  1. 这是染色体方向相同的

  1. 这是染色体方向相反的

  1. seqkit染色体到位
    1. 根据图片和统计表来决定哪些染色体序列需要倒位
    2. 参数
      1. r，对每个序列执行反向互补操作，即生成序列的互补链，并反转序列的方向。
      2. p，此选项指示 seqkit 对输入文件中的所有序列都应用 -r 参数（即对所有序列进行反向互补处理）。
      3. t，指定输入的序列类型为 DNA（默认），seqkit 将根据此信息进行相应的处理（如计算互补链时，基于 DNA 的碱基对规则 A ↔ T 和 C ↔ G）。
seqkit seq -r -p -t dna scaffold_7.fa  > scaffold_7.revcomp.fa
seqkit seq -r -p -t dna scaffold_11.fa > scaffold_11.revcomp.fa
seqkit seq -r -p -t dna scaffold_12.fa > scaffold_12.revcomp.fa
seqkit seq -r -p -t dna scaffold_13.fa > scaffold_13.revcomp.fa
seqkit seq -r -p -t dna scaffold_2.fa  > scaffold_2.revcomp.fa
  1. 合并到位后的染色体序列
cat \
  scaffold_1.fa \
  scaffold_2.revcomp.fa \
  scaffold_3.fa \
  scaffold_4.fa \
  scaffold_5.fa \
  scaffold_6.fa \
  scaffold_7.revcomp.fa \
  scaffold_8.fa \
  scaffold_9.fa \
  scaffold_10.fa \
  scaffold_11.revcomp.fa \
  scaffold_12.revcomp.fa \
  scaffold_13.revcomp.fa \
> all_scaffolds_merged.fa
  1. 根据id重命名
awk '
NR==FNR {
    if ($1 ~ /^chr/) map[$2] = "Chr" substr($1,4);
    next
}
{
    if ($0 ~ /^>/) {
        id = substr($0,2);
        if (id in map) print ">" map[id];
        else print $0;
    } else {
        print
    }
}
' chr_best_scaffold.txt all_scaffolds_merged.fa > all_scaffolds_merged_renamed.fa
  1. 重新全基因组比对一遍
    1. 将方向调整后的目标基因组和参考基因组比对，结果图能够清楚地看到两个基因组上所有染色体方向一致
    2. 不过，相较于这样一个一个拆分进行比对的有些肝的方法，如果能够在全基因组比对出图那里，获得更多的情报，比如染色体编号、更清晰的比对情况，就不需要拆分染色体这种法子了，写在这里作为参考使用
/home/soft/NGenomeSyn/bin/GetTwoGenomeSyn.pl \
    -InGenomeA Rsimsii.chr.fa \
    -InGenomeB all_scaffolds_merged_renamed.fa \
    -OutPrefix rs_ctg_last \
    -BinDir /home/soft/mummer-4.0.1 \
    -MinAlnLen 10000