使用PHAST分析保守性

发表于 更新于 分类于

简介

全基因组多序列比对的主要目的是揭示基因组中进化上保守的和非保守的区段以及它们的分布。但是，在碱基层次观察多个全基因组的多序列比对结果和序列保守性非常不方便；其次，许多功能区段的保守性介乎高度保守和完全不保守之间；再者，保守性要能予以方便的显示。上述三点使UCSC基因组浏览器使用两个软件PhastCons和PhyloP把由MULTIZ产生的多序列比对结果转换成两种单一的保守性计分和显示，这两个方法都基于已知的种系树结构和利用一个称作phylo-HMM的隐马尔可夫模型种系分析方法。

PhyloP只考虑比对的当前列而PhastCons同时也考虑比对列的相邻列，这使得PhastCons对于保守区段的出现更敏感而PhyloP对于保守区段的界定更有效。PhyloP和PhastCons之间的另一个区别是，PhyloP能够识别加速进化和保守进化的位点，它们分别产生正的和负的计分，而PhastCons的计分总是一个0至1之间的正值。

安装PHAST

PHAST（Phylogenetic Analysis with Space/Time models）官网：http://compgen.cshl.edu/phast/

PHAST依赖CLAPACK，需要先编译CLAPACK，然后编译PHAST时指定CLAPACK的路径

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
# 安装CLAPACK
cd ~/software/
wget http://www.netlib.org/clapack/clapack.tgz
tar zxf clapack.tgz
cd CLAPACK-3.2.1/
cp make.inc.example make.inc && make f2clib && make blaslib && make lib

# 安装PHAST
cd ~/software/
wget https://github.com/CshlSiepelLab/phast/archive/v1.5.tar.gz
tar zxf v1.5.tar.gz
cd phast-1.5/src/
make CLAPACKPATH=/home/chenwen/software/CLAPACK-3.2.1

# 将PHAST添加到环境变量
echo "export PATH=$HOME/software/phast-1.5/bin:\$PATH" >> ~/.bashrc
source ~/.bashrc

# 删除安装包
rm ~/software/clapack.tgz
rm -rf ~/software/CLAPACK-3.2.1
rm ~/software/v1.5.tar.gz

安装UCSC工具

本文除了使用PHAST外，还需要使用UCSC工具中的mafSplit、wigToBigWig和bigWigAverageOverBed。

1
2
3
4
5
6
mkdir -p ~/software/kent/
cd ~/software/kent/
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/mafSplit
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/wigToBigWig
wget https://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bigWigAverageOverBed
chmod +x *

运行phastCons

我们这里以斑马鱼（Danio rerio）、鲤鱼（Cyprinus carpio）、鲫鱼（Carassius auratus）基因组为例，斑马鱼的基因组作为参考基因组。假设我们已经通过多基因组比对，得到了Danio_rerio_mulitway.maf。

运行phyloFit生成MOD文件

1
phyloFit -i MAF Danio_rerio_mulitway.maf

会在当前文件夹生成phyloFit.mod文件。

运行phastCons计算保守性分数

phastCons一次只能处理一条染色体，所以我们需要使用mafSplit先分割MAF文件。

1
2
mkdir maf
mafSplit _.bed maf/ Danio_rerio_mulitway.maf -byTarget -useFullSequenceName

使用-byTarget参数分割MAF，第一个参数_.bed会被忽略掉，但是必须要有这个参数占位。-useFullSequenceName参数的作用是生成的文件名是染色体名。

接下来我们使用phastCons计算每个碱基的phastCons分数和保守的区块。

1
2
3
4
5
mkdir wig
mkdir bed
for i in maf/*.maf; do x=`basename $i .maf`; phastCons $i phyloFit.mod --target-coverage 0.25 --expected-length 12 --rho 0.4 --msa-format MAF --seqname $x --most-conserved bed/$x_most-cons.bed > wig/$x.wig; done
cat wig/*.wig >> Danio_rerio_mulitway.wig
cat bed/*_most-cons.bed >> Danio_rerio_most-cons.bed

Danio_rerio_mulitway.wig文件有斑马鱼全基因组单碱基phastCons分数，Danio_rerio_most-cons.bed文件有斑马鱼保守元件位置。

进一步分析

  1. 使用wigToBigWig将wig文件转换为二进制的bigWig文件后，就能在Genome Browser中展示了。 

    1
    wigToBigWig in.wig chrom.sizes out.bw

  2. 假设给定一条序列比如lncRNA，我们想计算它的平均phastCons分数，可以使用bigWigAverageOverBed实现。

    1
    bigWigAverageOverBed in.bw in.bed out.tab

  3. 保守元件去掉编码区的之后，就得到了保守的非编码元件（Conserved noncoding DNA elements，CNEs）。