基因家族进化分析

1 首先要阅读大量的文献，了解背景知识，选择你关心的，与基因家族功能相关的物种。

物种选择完之后，可以把这些物种提交到NCBI 的Taxonomy 数据库中，会得到这些物种的进化关系，即物种树。

2 然后从NCBI 的refseq 数据库中下载所需物种的DNA 和protein 序列，这时需要用到NCBI 的批量下载工具efetch。

######################################################

#This code is used for downloadind sequences from NCBI

######################################################

from Bio import Entrez

Entrez.email = "邮箱" # Always tell NCBI who you are

out_handle = open("result.txt", "w")

search_handle = Entrez.esearch(db="要搜索的数据库",term="关键字 ",usehistory="y")

search_results = Entrez.read(search_handle)

search_handle.close()

gi_list = search_results["IdList"]

count = int(search_results["Count"])

webenv = search_results["WebEnv"]

query_key = search_results["QueryKey"]

batch_size = 1

for start in range(0,count,batch_size):

end = min(count, start+batch_size)

print "Going to download record %i to %i" % (start+1, end)

fetch_handle = Entrez.efetch(db="检索的数据库", rettype="数据类型", retmode="文件类型",

retstart=start, retmax=batch_size,

webenv=webenv, query_key=query_key)

data = fetch_handle.read()

fetch_handle.close()

out_handle.write(data)

out_handle.close()

##################################################################################

3 序列处理，因为下载到的DNA 序列其实是mRNA 的序列，并不是编码区序列。所以要对应着protein 序列，在DNA 序列中找到所对应的conding region。

##########################################################################

#This code is used for turning the protein sequences back to DNA sequences

#####################################################################

4 序列比对，需要注意的是，要以对应密码子的形式进行比对，这里用到的是Prank 软件，默认参数

即可。直接用clustalw 比对，不准确。

prank -d=Exon.fa -o=序列文件 -codon –F

clustalw2 -INFILE=Exon.fa -OUTPUT=FASTA -OUTFILE= alignment.fa -SEQNOS=ON

#####################################################################

5 构建进化树，把序列比对的结果输入到做树的软件中，就能得到进化树。可用的软件有mega，

phylip，phyml，RAxML，fasttree

感觉mega 和RAxML 比较可信，而phylip 要求序列名字在10 个字符以

内，不喜欢这个，phyml 同样需要phylip 格式的序列比对文件。Fasttree 的特点是，速度非常快，很适合

物种多的时候使用，但是它的缺陷在于，只是近似的最大似然法，感觉准确性不高。

FastTree -gtr -nt <_alignment.fa> fasttree.nwk

-gtr 是模型，-nt 是DNA 序列，如果是蛋白序列，省略-nt 即可

raxmlHPC-MPI -s _prank.fa.1.fas -f a -x 1234 -N 1000 -p 1234 -m GTRGAMMA –n nwk -o

Oreochromis_niloticus,Danio_rerio

-s 是输入，-n 是输出的后缀，-o 是选择的外群，-N 是bootstroop 值

phyml -i alignment.phy -d aa -b 1000

-i 是输入，-d 是序列类型（aa 是氨基酸的意思），-b 是bootstroop 值

6 选择压力检测，用PAML 软件，可选branch model，site model， branch site model。

1> Site model：这个模型假设不同的位点可以受到不同的选择压力，而系统树上所有的枝受到的选择压

力相同，codeml.ctr 参数设置，model=0，Nsite=0,1,2,3,7,8。模型M0/M3，M1a/M2a，M7/M8 通过两

两比较LRT，看差异是否显著，确定是否有正选择事件发生。

2> Branch model：模型假设所有位点受到的选择压力相同，不同的枝受到的选择压力可以不同，参数设

置为，model=2，Nsite=0。与零假设模型（model=2,Niste=0）通过LRT 比较，看差异是否显著，确定

正选择是否发生

3> Branch site model：模型假设某一特定分支(foreground)与其他所有分支(background)上的所有位点具有

不同的ω值，即该模型假设正选择仅仅发生在特定枝上的少数位点上，参数设置model=2，Nsite=2.

零假设模型（fix omega=1）与正选择模型（fixomega=0）通过LRT 进行比较

################################################################

#biopython for paml

################################################################

from Bio.Phylo.PAML import codeml

cml = codeml.Codeml()

cml.alignment = "alignment.fa"

cml.tree = "species.nwk"

cml.out_file = "out"

cml.working_dir = "./"

cml.set_options(seqtype=1,

verbose=0,

noisy=0,

RateAncestor=0,

model=2,

NSsites=[2],

CodonFreq=2,

cleandata=1,

fix_alpha=1,

kappa=4.54006)

results = cml.run()

更多精彩欢迎关注