前言:編程
這是coursera課程:Probabilistic Graphical Models 上的第二個實驗,主要是用貝葉斯網絡對基因遺傳問題進行一些計算。具體實驗內容可參考實驗指導教材:bayes network for genetic inheritance. 你們能夠去上面的連接去下載實驗材料和stard code,如實驗內容有難以理解的地方,歡迎私底下討論。下面是隨便寫的一些筆記。數組
完成該實驗須要瞭解一些遺傳方面的簡單知識,可參考:Introduction to heredity(基因遺傳簡單介紹)網絡
關於實驗的一些約定和術語:app
同位基因(也叫等位基因)中有一個顯現基因(dominant allele)和一個隱性基因(recessive allele)。dom
homozygous:純合體ide
heterozygous:雜合體函數
一個factor表不只能夠表示聯合機率,也能夠表示條件機率,這裏約定規則是: 最左邊的那個元素表示條件輸出,即它是在後面元素出現條件下的輸出。編碼
關於變量的命名,由於在圖模型中,會出現不少變量,這些變量在實際編程中是要取名字的,通常是按照某種順序規則來取,本次實驗的規則是:依照家譜圖中每一個人的名字依次取,假設有k我的,則每一個人對應了2個factor,因此前k個變量表示這k我的的基因型標量,後k個變量表示這k我的的表現型變量。若是把每一個人的基因對應的2個染色體分開考慮,則須要3k個變量,前k個表示每一個人的第一條染色體對應處的基因,中間k個表示每一個人的第二條染色體對應處的基因,後k個變量和之前的同樣,表示這k我的的表現型變量。spa
程序中的alleleFreqs是個n*1的向量,n爲單個等位基因的個數。code
程序中的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能夠嵌套,能夠重疊
參考資料: