全基因组比对策略

发表于 | 更新于 | 分类于

两个基因组的比对主要在于大片段序列的比对正确性，这和小片段比对（blast）、多序列比对（mafft）不同。基因组层面的比对方法常用的有 lastz 和 mummer。

lastz 流程

axt文件的处理

在使用 lastz 比对完两个基因组之后，生成了以 axt 结尾的文件（比对时可以分割染色体）。axt 文件即记录了2个基因组的比对位置和对应比对的序列，如下：其中 chr21 19645 19675 Chr01 318732 318762是比对位置。

1
2
3
1 chr21 19645 19675 Chr01 318732 318762 + 2328
CACACAGACATACACACATGCACACACACGC
CACACAAACACACACACACACACACACACAC

我们需要把这一行的信息挑选出来。

1
grep chr xxx.axt |cut -f 2,3,4,5,6,7 -d " " > xxx.pos

得到比对位置之后，我们下一步对结果进行过滤。

片段长度过滤

太短的比对无法真实地体现比对位置的正确性，一般比对长度要大于 1000bp

1
>>>awk '$3-$2>1000{print $0}' xxx.pos >xxx.filter.pos

位置信息过滤

理想情况下，应该是一个基因组的染色体比对到另一个基因组的染色体上，而实际情况会有很多碎片散落在其他染色体上，我们挑选出来的位置文件包含了一条染色体比对到另一个基因组所有染色体的信息，所以我们要进行过滤，看看有没有一一对应比对的染色体位置信息。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
>>> head -5 xxx.filter.pos
	chr11 1016468 1018558 Chr01 1589564 1591669
	chr11 2208239 2209250 Chr01 3172963 3173855
	chr11 4940937 4942025 Chr01 6355782 6356928
	chr11 6016505 6017738 Chr01 8231646 8232851
	chr11 6770249 6771299 Chr01 10322020 10323107
>>> cut -f 4 xxx.filter.pos -d " "|uniq -c
	142 Chr01
    275 Chr02
    318 Chr03
    141 Chr04
    216 Chr05
    227 Chr06
    167 Chr07
    313 Chr08
    245 Chr09
    212 Chr10
    273 Chr11
    223 Chr12
    193 Chr13
    236 Chr14
    303 Chr15
    170 Chr16
    304 Chr17
   5130 Chr18
    240 Chr19
    333 Chr20
    209 Chr21
    289 Chr22
    141 Chr23
    142 Chr24
>>> grep Chr18 xxx.filter.pos > xxx.result

可以看到 chr11 的绝大部分序列比对到了 Chr18（5130），我们由此得知 chr11-Chr18是一对儿比较好的共线性染色体。

作图

我们使用LINKVIEW 进行比对的可视化，输入文件就是上面的 xxx.result。

1
python LINKVIEW.py xxx.result -o xxx # https://github.com/YangJianshun/LINKVIEW

mummer 流程

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
>>> nucmer genome1.fa genome2.fa
>>> delta-filter -i 89 -l 1000 -1 out.delta > out.delta.filter #比对率大于89%，对比长度大于1000
>>> show-coords out.delta.filter > out.coords
>>> grep -P  "chr1\s+" out.delta.coord|awk '{print $12,$13}' |sort |uniq -c
	  6 chr1 Chr01
      6 chr1 Chr02
      7 chr1 Chr03
     12 chr1 Chr04
   6688 chr1 Chr05
      4 chr1 Chr06
     12 chr1 Chr07
      9 chr1 Chr08
      4 chr1 Chr09
     17 chr1 Chr10
      8 chr1 Chr11
      4 chr1 Chr12
     12 chr1 Chr13
      7 chr1 Chr14
     23 chr1 Chr15
      8 chr1 Chr16
      8 chr1 Chr17
      7 chr1 Chr18
     11 chr1 Chr19
     10 chr1 Chr20
      4 chr1 Chr21
      7 chr1 Chr22
      3 chr1 Chr23
      6 chr1 Chr24

可以看到 chr1 的绝大部分序列比对到了 Chr05（6688），我们由此得知 chr1-Chr05是一对儿比较好的共线性染色体。

作图

我们使用LINKVIEW 进行比对的可视化。

1
2
3
grep -P  "chr1\s+Chr05" out.delta.coord|awk '{print $12,$1,$2,$13,$4,$5}' >chr1_chr5.pos
sed -i 's/ /\t/g' chr1_chr5.pos
python LINKVIEW.py chr1_chr5.pos -o xxx # https://github.com/YangJianshun/LINKVIEW