序列间间隔的计较
在此,我们把序列间间隔界说为每个位点产生核酸代替的期望值。同时,我们界说p间隔为序列间差异位点的比例。由于p间隔并不是进化时间的线性方程,所以p间隔一般只用于阐明高度相似的序列,即p<5%。
很容易想到,核酸代替进程中每个位点的下一状态只和当前状态有关,而和前一状态无关。这是典范的马尔可夫性质,所以持续型马尔可夫链被引入作为统计学模子利用。每个位点的代替被看作是马尔可夫链,A、T、C、G四个核酸被看作是个中的四个元素。
下一个需要办理的问题就是位点代替率矩阵的界说,这也是诸多序列间间隔计较要领的主要区别。最经典的模子一般有JC69、K80、F81、HKY85等,此处我们利用最简朴的JC69模子。
JC69模子由Jukes和Cantor于1969年提出,他们假设每种核酸都用沟通的概率λ被其它三种核酸所代替。
为了更好地研究马尔可夫链,我们界说跃迁概率pij为t时间后核酸i变为核酸j的概率,跃迁概率矩阵如下所示:
从p0(t)、p1(t)的形式可以很显然地看出一些重要的性质。首席,矩阵P(t)的每一行的和都是1,也就是说序列中的每个核酸必然是四种核酸中的一种。其次,P(0)=I,这是没有核酸位点的改变,也就是进化没有产生。第三,只思量λt的值,也是说只要λt牢靠,所有环境对我们来讲没有区别。第四,pij(t)趋近于1/4,该马尔可夫模子的极限漫衍是(1/4, 1/4, 1/4, 1/4),也就是说此时方针序列每个位点都以1/4的概率随机漫衍四种核酸。从前面的论述,我们已经知道每种核酸的总代替率是3λ,那么假如两个基因在时间t之前分隔,则这两个序列之间的间隔d=3λt。简朴的讲,我们假定四种碱基之间的突变概率是一致的,但实际环境虽然不是,因为这个跟化学布局的差别有干系。这个环境下,我们可以找出两个物种基因内部呈现差此外碱基,再调核对应的三个碱基所编码的氨基酸是否是一样的。因为三个碱基有64种大概编码,而个中一些差异的编码大概对应沟通的20种中的一种氨基酸,这时我们界说为一个同义变异Ks,同理,假如碱基变异后所编码的氨基酸呈现了差异,这时我们就界说为异义变异Ka。等概率变异环境,则计较在一个基因中相应的个数,而且计较Ka/Ks的值作为基因进化的指标。一般认为,同义突变不受自然选择,而非同义突变则受到自然选择浸染。假如Ka/Ks>1,则认为有正选择效应。假如Ka/Ks=1,则认为存在中性选择。假如Ka/Ks<1,则认为有纯化选择浸染,即很有大概突变造成了卵白质的失活而影响了机体的保留最终被进化浸染所裁减。
对样本1600多个的序列比对,我们要研究,进化速率高的基因是否会相应的有更高的,可能更低的表达量。
数据是从Human Genome数据库中随机选取的约莫1600个基因的编码区在三个物种人类、大猩猩和猴子中的同源序列的比对功效。数据库中包括了45个物种所有的编码序列,这里选取这约莫1600个基因,主要是因为冷泉港尝试室已经对其表达量举办丈量,后头可以将基因表达量作为一个影响因素举办阐明。这里需要指出的是冷泉港尝试室对三个物种差异性此外6个样本举办了丈量,性别间的区别也会成为较量内容。
文件中序列的名目如下:
>NM_203342_hg19_1_17 54 0 0 chr1:29320001-29320054+
ATGCACTGCAAGGTTTCTTTGTTGGATGACACAGTTTATGAATGTGTTGTGGAG
个中“NM_203342”是Reference ID,代表基因的名字。hg19则代表了物种,hg19是human,panTro2是大猩猩,rheMac2是猴子。1代表了17个外显子中的第一个,后头的chr1:29….等信息表白的是这一段序列在基因组上的坐标,最后为该基因的序列。
界说物种间基因在序列程度的演化速度为Ka/Ks,而在表达的差别为(log2(X)-log2(Y))/(log2(X)+log2(Y)),X指基因在物种A的表达量,Y指基因在物种B的表达量。
这时候首先选取输出功效中的基因编号一栏,与表达量数据中的基因编号栏举办match,而且再反match回kaks数据中,这样可以或许找到在两者中都具有的基因,以用来做相关性阐明。获得的功效约莫为900个数据对。此时别离做散点图,而且ka/ks>0,因此调解横轴坐标。同时也可以利用cor(x,y)计较两个变量间的相干系数。
下图是人类和猴子的进化速率与表达差别,横轴为进化速率,纵轴为表达差别。
大猩猩和猴子的进化速率与表达差别,横轴为进化速率,纵轴为表达差别如下:
人类和大猩猩的进化速率与表达差别,横轴为进化速率,纵轴为表达差别
二基因进化速率与成果
既然基因的进化速率和表达差别不存在明明差别,下面我们就进化速率与成果稍作接头。操作方才生成的KaKs矩阵,计较Ka/Ks并举办排序,之后操作GenBank所提供的API接口,获得进化速率最快的30个基因的成果。
ID | Vs species | Ka/Ks | Function |
NM_000798 | rheMac2 | 54.99143243 | D5 subtype of the dopamine receptor |
NM_002341 | rheMac2 | 14.08734465 | TNF superfamily, member 3 |
NM_152308 | panTro2 | 13.59065657 | chromosome 16 open reading frame 75 |
NM_015326 | rheMac2 | 13.04769249 | SLIT-ROBO Rho GTPase activating protein 2 |
NM_000190 | rheMac2 | 11.99951905 | hydroxymethylbilane synthase |
NM_001005487 | panTro2 | 11.60989607 | olfactory receptor, family 13, subfamily G, member 1 |
NM_016120 | panTro2 | 11.38441703 | ring finger protein, LIM domain interacting |
NM_001031809 | rheMac2 | 11.13806742 | membrane-spanning 4-domains, subfamily A, member 3 |
NM_001900 | panTro2 | 10.76315242 | cystatin D |
NM_173644 | panTro2 | 9.78075552 | chromosome 20 open reading frame 197 |
NM_001012456 | panTro2 | 9.665625341 | Sec61 gamma subunit |
NM_002099 | panTro2 | 9.037433303 | glycophorin A (MNS blood group) |
NM_001006655 | panTro2 | 8.863350751 | family with sequence similarity 149, member A |
NM_207352 | panTro2 | 8.863350751 | cytochrome P450, family 4, subfamily V, polypeptide 2 |
NM_198492 | panTro2 | 8.460999047 | C-type lectin domain family 4, member G |
NM_001004741 | panTro2 | 8.406537948 | olfactory receptor, family 5, subfamily M, member 10 |
NM_001122834 | rheMac2 | 8.315430071 | hedgehog acyltransferase |
NM_006072 | panTro2 | 8.237683925 | chemokine (C-C motif) ligand 26 |
NM_001146333 | panTro2 | 7.835034285 | sulfatase modifying factor 2 |
NM_001007551 | panTro2 | 7.788487115 | cancer/testis antigen family 45, member A5 |
NM_001007272 | panTro2 | 7.712338837 | dual specificity phosphatase 13 |
NM_173487 | rheMac2 | 7.619466474 | chromosome 4 open reading frame 33 |
NM_001912 | panTro2 | 7.541731476 | cathepsin L1 |
NM_152548 | panTro2 | 7.511077919 | family with sequence similarity 81, member B |
NM_058192 | panTro2 | 7.478620447 | RNA pseudouridylate synthase domain containing 1 |
NM_030636 | rheMac2 | 7.47601075 | endonuclease/exonuclease/phosphatase family domain containing 1 |
NM_006117 | panTro2 | 7.390280883 | peroxisomal D3,D2-enoyl-CoA isomerase |
NM_000504 | panTro2 | 7.254875852 | vitamin K-dependent coagulation factor X of the blood coagulation cascade |
NM_002958 | rheMac2 | 7.216840091 | RYK receptor-like tyrosine kinase |
NM_016201 | rheMac2 | 7.203190433 | angiomotin like 2 |

下图别离为长染色体上基因演化速率的漫衍,以及对局部的放大
可以看出,常染色体上基因演化速度平均值较高,但大多会合在0~0.5之间,X染色体上基因的会合水平比常染色体高。
下图横轴为human差异性别间表达差别,纵轴为human和chimpanzee间表达差别。若存在某一值显著大于另一值,则图中点应该漫衍在y=x直线的上方或下方。实际做图后发明,并没有明明的漫衍
纪律,t检讨也证明两组值没有显著差别。
[attach]1885[/attach]
基因成果聚类与通过表达谱举办的物种亲缘性阐明
凡是生物个别基因都有很多差异成果,可是有很多的基因大概有临近的成果,也即它们大概编码了同一个卵白质的差异构成部门,因而其相应的表达量会泛起同步现象。同时,通过基因之间的序列相似间隔,可能我们这里操作同一个基因在三个物种中的表达量间隔作为权衡指标,可以或许成立起各个物种之间的亲缘干系,也就是进化树的构建。 因此这里可以选取基因表达数据的前100个基因作为示范,对其举办条理聚类,以表格中第4:39列为坐标,计较间隔矩阵,并操作hclust(d,”average”)举办聚类,而且凭据在三个物种中大概为高表达和低表达两种程度,即共8中组合举办分别,可获得如下功效图:
如上图所示,一大部门的基因被归为了个中的一类,大概的表明为这一些基因的成果在三个物种傍边较为基本与重要,因此根基上并没有太大差别。而其他少量的其他基因则大概为三个物种在进化中演化出了差异的机制因此呈现了表达的差别。 接着,我们试图鉴定人类是与大猩猩照旧猴子的亲缘性更为临近,而权衡的指标为沟通的基因在三个物种傍边的表达量间的差别,也等于间隔。通过截取源文件1:100行,4:39列,而且转化为矩阵,并举办转置后,构成36个样本100维的样本点。计较间隔矩阵,并做相应的条理聚类阐明,获得的功效如下图:
从图上的功效来看,猴子的基因的表达量应该是与人类和大猩猩的表达量有显著的差异,因此其亲缘干系应该更为远一些,而大猩猩与人类样本的数据则稠浊在一起,说明两个物种间的相应基因的表达量差别并不明明,因此亲缘干系大概更近。总结文中对现有的基因序列比对的数据,通过计较ka/ks作出与表达量的相关干系图,通过对数据中的某些数据总结,找出具有某些特性的基因背后的生物学意义,另外还做了X染色体和常染色体表达的差别等阐明,同时也试图按照表达量这一个指标别离对基因和物种举办聚类阐明。本文回收的阐明要领只是很多要领中的个中一部门,差异的要领大概会有差异的功效。而且对沟通的数据,也可以从其他的角度探讨其他大概埋没的信息。把握更多的阐明手段,将可以或许使得我们更好的处理惩罚大批量数据,总结数据,而且从中挖掘出我们但愿获得的信息。
[1] Yang, Ziheng.Computational molecular evolution. Oxford: Oxford University Press, 2006.
[2]统计建模与R软件 薛毅 2006
[3] http://blog.sciencenet.cn/home.php?mod=space&uid=223428&do=blog&id=511891
原始数据以及剧本下载地点
http://dl.dbank.com/c063st07vi
