两基因组比对

发表于 更新于 分类于

本教程根据UCSC Genome Browser提供的教程略作修改。原教程地址http://genomewiki.ucsc.edu/index.php/Whole_genome_alignment_howto。

软件安装

安装UCSC工具

Jim Kent创建UCSC Genome Browser，开发了一系列工具，这套工具也叫Kent Utils。我们这里需要使用到有faSize、faToTwoBit、faSplit、axtChain、chainMergeSort、chainPreNet、chainNet、netSyntenic、netToAxt、axtSort、axtToMaf，直接下载相应的二进制文件即可。

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
mkdir -p ~/software/kent
cd ~/software/kent/
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faSize
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faToTwoBit
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/faSplit
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/axtChain
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainMergeSort
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainPreNet
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/chainNet
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/netSyntenic
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/netToAxt
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/axtSort
wget https://hgdownload.soe.ucsc.edu/admin/exe/linux.x86_64/axtToMaf
chmod +x *
echo "export PATH=$HOME/software/kent:\$PATH" >> ~/.bashrc
source ~/.bashrc

安装LASTZ

原教程中使用的是BLASTZ，BLASTZ已经被LASTZ替代了，而LASTZ也停止更新了。这里我们使用Github上的可以被现代编译器编译的版本。

1
2
3
4
5
6
7
8
cd ~/software/
wget https://github.com/lastz/lastz/archive/refs/tags/1.04.15.tar.gz
tar zxf 1.04.15.tar.gz
cd lastz-1.04.15/
make
cp src/lastz ~/software/kent/
rm -rf ~/software/lastz-1.04.15
rm ~/software/1.04.15.tar.gz

第0步：准备

基因组比对首先需要屏蔽基因组上的重复序列。我们以斑马鱼（Danio rerio）和鲤鱼（Cyprinus carpio）的基因组为例，从Ensembl下载软屏蔽的基因组文件。如果没有屏蔽了重复序列的基因组，可以使用RepeatModeler和RepeatMasker自己屏蔽。

1
2
3
4
5
6
7
8
9
mkdir -p ~/Danio_rerio_Cyprinus_carpio
cd ~/Danio_rerio_Cyprinus_carpio/
wget http://ftp.ensembl.org/pub/release-104/fasta/danio_rerio/dna/Danio_rerio.GRCz11.dna_sm.primary_assembly.fa.gz
gunzip Danio_rerio.GRCz11.dna_sm.primary_assembly.fa.gz
mv Danio_rerio.GRCz11.dna_sm.primary_assembly.fa Danio_rerio.fa

wget http://ftp.ensembl.org/pub/release-104/fasta/cyprinus_carpio/dna/Cyprinus_carpio.common_carp_genome.dna_sm.toplevel.fa.gz
gunzip Cyprinus_carpio.common_carp_genome.dna_sm.toplevel.fa.gz
mv Cyprinus_carpio.common_carp_genome.dna_sm.toplevel.fa Cyprinus_carpio.fa

将fa文件转换生成2bit文件，计算每条染色体的长度，用于后续的分析。

1
2
3
4
faToTwoBit Danio_rerio.fa Danio_rerio.2bit
faToTwoBit Cyprinus_carpio.fa Cyprinus_carpio.2bit
faSize Danio_rerio.fa -detailed > Danio_rerio.sizes
faSize Cyprinus_carpio.fa -detailed > Cyprinus_carpio.sizes

lastz本身不支持并行，所以我们将斑马鱼基因组按照染色体切分，手动并行。

1
2
mkdir fa
faSplit byName Danio_rerio.fa fa

第1步：使用lastz比对

1
2
mkdir axt
for i in fa/*.fa; do prefix=$(basename $i .fa); lastz $i Cyprinus_carpio.fa O=400 E=30 K=3000 L=3000 H=2200 T=1 --format=axt --ambiguous=n --ambiguous=iupac > axt/${prefix}.axt; done

上面用了shell代码对斑马鱼的每条染色体进行比对，这里贴一条看得更加清楚的代码。

1
lastz fa/chr1.fa Cyprinus_carpio.fa O=400 E=30 K=3000 L=3000 H=2200 T=1 --format=axt --ambiguous=n --ambiguous=iupac > axt/chr1.axt

第2步：Chaining

从axt文件得到chain文件

1
2
mkdir chain 
for i in axt/*.axt; do prefix=$(basename $i .axt); axtChain $i Danio_rerio.2bit Cyprinus_carpio.2bit chain/${prefix}.chain -minScore=3000 -linearGap=medium

上面用了shell代码获得斑马鱼的每条染色体的chain文件，这里贴一条看得更加清楚的代码。

1
axtChain axt/chr1.axt Danio_rerio.2bit Cyprinus_carpio.2bit chain/chr1.chain -minScore=3000 -linearGap=medium

合并chain文件，并过滤

1
2
chainMergeSort chain/*.chain > all.chain
chainPreNet all.chain Danio_rerio.sizes Cyprinus_carpio.sizes all.pre.chain

第3步：Netting

1
chainNet all.pre.chain -minSpace=1 Danio_rerio.sizes Cyprinus_carpio.sizes stdout /dev/null | netSyntenic stdin noClass.net

第4步：Maffing

1
2
netToAxt noClass.net all.pre.chain Danio_rerio.2bit Cyprinus_carpio.2bit stdout | axtSort stdin Danio_rerio.Cyprinus_carpio.axt
axtToMaf Danio_rerio.Cyprinus_carpio.axt Danio_rerio.sizes Cyprinus_carpio.sizes Danio_rerio.Cyprinus_carpio.maf -tPrefix=Danio_rerio. -qPrefix=Cyprinus_carpio.

总结

经过这4步，就得到了基因组两两比对的maf文件Danio_rerio.Cyprinus_carpio.maf。原教程还有第5步计算Phastcons，不过通常的做法是根据多基因组比对的结果计算Phastcons，所以这部分我们另外讲。我写了一个Python脚本，让大家方便地运行以上步骤。