机器学习&数据挖掘笔记_18（PGM练习二：贝叶斯网络在遗传图谱在的应用）

 

　　前言：

　　这是coursera课程：Probabilistic Graphical Models 上的第二个实验，主要是用贝叶斯网络对基因遗传问题进行一些计算。具体实验内容可参考实验指导教材：bayes network for genetic inheritance. 你们能够去上面的连接去下载实验材料和stard code，如实验内容有难以理解的地方，欢迎私底下讨论。下面是随便写的一些笔记。

　　完成该实验须要了解一些遗传方面的简单知识，可参考：Introduction to heredity(基因遗传简单介绍)

 

　　关于实验的一些约定和术语：

　　同位基因(也叫等位基因)中有一个显现基因（dominant allele）和一个隐性基因（recessive allele）。

　　homozygous：纯合体

　　heterozygous：杂合体

　　一个factor表不只能够表示联合几率，也能够表示条件几率，这里约定规则是： 最左边的那个元素表示条件输出，即它是在后面元素出现条件下的输出。

　　关于变量的命名，由于在图模型中，会出现不少变量，这些变量在实际编程中是要取名字的，通常是按照某种顺序规则来取，本次实验的规则是：依照家谱图中每一个人的名字依次取，假设有k我的，则每一个人对应了2个factor，因此前k个变量表示这k我的的基因型标量，后k个变量表示这k我的的表现型变量。若是把每一个人的基因对应的2个染色体分开考虑，则须要3k个变量，前k个表示每一个人的第一条染色体对应处的基因，中间k个表示每一个人的第二条染色体对应处的基因，后k个变量和之前的同样，表示这k我的的表现型变量。

　　程序中的alleleFreqs是个n*1的向量，n为单个等位基因的个数。

　　程序中的alphaList是个m*1的向量，m为2个等位基因组合的个数。

　　例如，某个基因家谱图可以下：

 　　

 

　　matlab知识：

　　matlab在使用元胞数组时，须要先用cell来定义大小。对元胞使用小括号访问时返回的是对应位置的数据类型，而不是其内容，要想获得其内容，须要用大括号访问。

　　结构体的定义要么直接对其赋值，这时候需调用其field，要么用关键字struct来定义。

　　tf = isfield(S, 'fieldname')：

　　tf=1表示结构体成员S中含有变量fieldname.

 

　　实验中的一些函数：

　　phenotypeFactor = phenotypeGivenGenotypeMendelianFactor(isDominant, genotypeVar, phenotypeVar)：

　　练习1须要完成的内容。该函数功能是生成一个factor，该factor每一项是一个条件几率，表示在基因类型变量（genotypeVar）下的表现类型（phenotypeVar）的几率,即P(person’s phenotype | person’s genotype)。其几率值是按照孟德尔模型(mendelian model)来分配的。

　　v = GetValueOfAssignment(F, A)：

　　获得facotr F在assignment A处的值。

　　F = SetValueOfAssignment(F, A, v)：

　　该函数实现对factor F的一个assignment对应的值进行更新，更新后的值为v.注意必定要将这个函数的返回值赋给F才有效。

　　phenotypeFactor = phenotypeGivenGenotypeFactor(alphaList, genotypeVar, phenotypeVar)：

　　练习2须要完成的内容。这个函数也是实现表现型在基因型下的条件几率，但此时几率不必定为0或者1，由于不严格知足孟德尔模型，因此参数列表中多了一个alphaList变量，该变量存储的是各类基因型下出现某种病的几率值。genotypeVar, phenotypeVar分别为基因类型变量和表现型变量，仅仅是个变量名而已。

　　[allelesToGenotypes,genotypesToAlleles]= generateAlleleGenotypeMappers(numAlleles)：

　　numAlleles为等位基因的个数，假设个数为n，则allelesToGenotypes为一个n*n的矩阵，表明的是2不一样等位基因组合后的编码(行号和列号分别表明这2个基因序号)，由于这个编码矩阵是对称的，因此只有m=n*(n-1)/2+n个编码。genotypesToAlleles为一个m*2大小的矩阵，每一行的两个元素表明该基因编码下的2个基因序号。

　　[allelesToGenotypes,genotypesToAlleles]= generateAlleleGenotypeMappers(numAlleles)：

　　用同位基因的个数numAlleles做为参数，来产生2个矩阵。其中allelesToGenotypes(i,j)表示编号为i的同位基因和编号为j的同位基因在基因类型(每一个基因类型由2个基因构成)中的序号。 genotypesToAlleles由2列构成，每一行对应一个基因类型的序号，这2个列元素表示该基因序列对应的2个同位基因的编号。

　　genotypeFactor = genotypeGivenAlleleFreqsFactor(alleleFreqs, genotypeVar)：

　　练习3须要完成的内容。某个基因可能有多个等位基因，而每一个等位基因的几率不一样，几率值保存在alleleFreqs中。这个函数的做用就是由这些等位基因组合产生的每一个基因类型的几率，这至关于一个先验分布，由于没有去考虑由父母亲的基因交叉获得下一代。genotypeVar表示这个基因类型的变量名。

　　genotypeFactor=genotypeGivenParentsGenotypesFactor(numAlleles, genotypeVarChild, genotypeVarParentOne, genotypeVarParentTwo):

　　练习4须要完成的内容。该函数的做用是已知父亲和母亲某个等位基因的分布，求出儿子对应基因的几率分布。其中numAlleles是等位基因的个数，genotypeVarChild, genotypeVarParentOne, genotypeVarParentTwo分别表示儿子，父亲，母亲对应的基因变量。输出的genotypeFactor是一个变量值val为3维矩阵的factor. 这个函数与上面那个函数不一样，这里的输出factor是在考虑了父母亲factor状况下的。

　　factorList = constructGeneticNetwork(pedigree, alleleFreqs, alphaList)：

　　练习5须要完成的内容。在给定家谱(家谱中包含了成员的名字，父母亲等信息)pedigree、等位基因的几率分布alleleFreqs、不一样基因组合致使生并的几率分布alphaList后，就能够求出这个家庭中每一个成员的2个factor：一个是关于本身的基因类型的，这由其父母亲的基因类型决定;另外一个是关于本身的表现类型的，这与他的基因类型有关。在完成该函数的代码中会调用前面的函数genotypeGivenAlleleFreqsFactor(), genotypeGivenParentsGenotypesFactor(), phenotypeGivenGenotypeFactor().因此前面不少工做都是为后面服务的。

　　phenotypeFactor=phenotypeGivenCopiesFactor(alphaList, numAlleles, geneCopyVarOne, geneCopyVarTwo, phenotypeVar)：

　　练习6所需完成的内容。geneCopyVarOne和geneCopyVarTwo一个来自母亲的基因，一个来自父亲的的基因。phenotypeVar为下一表明现型标量。该函数的做用是给定来自父母亲的2个基因，求这个组合下的下一代出现病的几率，固然了，由于不服从孟德尔模型，因此还须要参考几率值alphaList.

　　geneCopyFactor = childCopyGivenFreqsFactor(alleleFreqs, geneCopyVar):

　　该函数的做用是得得到染色体中一个等位基因的几率，该几率来源于通用的先验alleleFreqs. 该factor表格大小为n*1.

　　geneCopyFactor = childCopyGivenParentalsFactor(numAlleles, geneCopyVarChild, geneCopyVarOne, geneCopyVarTwo):

　　该函数的做用也是得得到染色体中一个等位基因的几率，只是这取决于父母亲各自那个等位基因。因此这个factor表格大小为n*3，且每一个条件几率只可能取3个     值：0, 0.5, 1.

　　factorList = constructDecoupledGeneticNetwork(pedigree, alleleFreqs, alphaList):

　　实验7所需完成内容。该函数是构造Decoupled的贝叶斯网络，这里的Decoupled是指家谱中每一个人染色体的2个等位基因分布来考虑。因此每一个人都对应有3个factor，一个是基因1的几率，n*3大小(若是家谱中没有父母亲，则是n*1大小)，一个是基因2的几率，n*3大小(若是家谱中没有父母亲，则是n*1大小)，另外一个是对应的表现型分布2n^2*3大小。

　　phenotypeFactor=constructSigmoidPhenotypeFactor(alleleWeights, geneCopyVarOneList, geneCopyVarTwoList, phenotypeVar):

　　这是练习8的内容。一种病可能由多个基因决定，而每一个基因又可能有多个等位基因。这里的alleleWeights矩阵中alleleWeights{i}(j)表示的意思是第i个基因的第j号等位基因对这种病的影响权值。当把某种病对应的全部基因的等位基因权值都考虑进来，那么得这种病的几率则与这些权值之和有关，本实验假设它们最后的关系是服从sigmoid分布。geneCopyVarOneList和geneCopyVarTwoList都是列向量(由于一种病可能由多个基因决定)，因此在计算phenotypeFactor.var时须要这2个向量转置。

 

　　课件中的一些内容：

　　贝叶斯推理的3种形式：因果推理，证据推理，内部因果推理(注意explean away现象)。

　　active trail(路径中不能出现V型结构，除非V型中间节点或者其子孙节点被观测到)。

　　变量之间的独立性判断、条件独立性判断(能够经过联合几率与某2个factor的乘积的关系来判断)。

　　一个图可以被因式分解，则有可能推导出独立性。若是有D-separation，则必定能够推导出独立性。

　　分布P可以分解图G ，则G是P的一个I-map，反之，若是G是P的一个I-map，则P能分解G.

　　为了简化时序模型，通常会有马尔科夫和时序不变性这2个假设。

　　时序转移模型(Template Transition Model)，能够包含内部时序转移和外部时序转移。

　　Ground Network(unrolled network)指的是不一样时间内的结构依次复制。

　　2TBN是一个动态贝叶斯网络，它有一个初始化以及unrolled network,因此有2个转移几率。

　　HMM中的状态转移几率是稀疏的。

　　Plate model能够嵌套，能够重叠

 

　　参考资料：

       Introduction to heredity(基因遗传简单介绍)

       coursera课程：Probabilistic Graphical Models