真核生物基因組的基因分析和預測

一、摘要

1. 加深基因預測基本原理的理解（如
2. 密碼子的偏好性、內含子外顯子剪下識別序列等）；
3. 瞭解同源基因預測的意義所在；
4. 熟悉已有的基因預測的使用（如GenScan、GeneWise等）；

二、材料和方法

1、硬體平臺

處理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
安裝記憶體(RAM)：16.0GB

2、系統平臺

Windows 8.1、Ubuntu

3、軟體平臺

1. BLAST
2. GenScan
3. GeneWise
4. Python3.5
5. Biopython

4、資料庫資源

NCBI資料庫：https://www.ncbi.nlm.nih.gov/
Uniprot資料庫：http://www.uniprot.org/
Genecode資料庫：http://www.gencodegenes.org/

5、研究物件

真菌基因組草圖(課堂提供)
小鼠基因組(fasta)+小鼠基因組註釋檔案(gff)
http://www.gencodegenes.org/mouse_releases/current.html

6、方法

應用部分：

建立本地BLAST資料庫
下載安裝BLAST，使用makeblastdb對基因組草圖建立BLAST資料庫。
已知蛋白序列的下載

根據基因組的物種來源（真菌），從UniProt資料庫下載近緣物種（真菌類）已知蛋白序列，以FASTA格式儲存。
同源基因搜尋
使用tblastn程式，把已知蛋白質序列和基因組草圖序列建立的本地BLAST
資料庫進行比對
Genewise建模
程式設計解析blast結果，取出蛋白質序列與對應基因組區間(上下延伸1000bp)，放入Genewise進行基因結構建模

程式設計建模部分：

基因組不同區域序列的獲取
根據小鼠gff註釋檔案，編寫程式解析小鼠基因組fasta，提取CDS 區域序列、基因內含子序列、基因間序列、外顯子-內含子剪下識別位點附近的序列。
三聯體密碼子偏好性的分析

CDS 區域、基因內含子、基因間區域的序列進行三聯體密碼子偏好性的統計分析。
剪下識別序列鹼基規律分析
對剪下識別位點前後20bp 的序列進行位點鹼基統計，建立序列譜。
剪下識別序列的計算鑑別
分別計算陽性經驗分值和陰性經驗分值，基因結構的計算鑑別。

三、結果

應用部分：

1、建立本地BLAST資料庫

下載安裝BLAST，使用makeblastdb對基因組草圖建立BLAST資料庫。

makeblastdb -in scaffold.fasta -input_type fasta -dbtype nucl -parse_seqids -out fungi

圖表 1建立本地資料庫

2、已知蛋白序列的下載

方法一：
進入uniprot下載頁面，點進按“分類學劃分”Taxonomic_divisions，進入ftp下載頁面，其中sprot為Swiss-prot稽核過的蛋白質(Reviewed)，trembl為TrEMBL資料庫，沒有稽核過的蛋白質(Unreviewed)，選擇下載uniprot_sprot_fungi.dat.gz
ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/taxonomic_divisions/uniprot_sprot_fungi.dat.gz


方法二：
進入uniprot，搜尋taxonomy:fungi，然後選擇reviewed，下載fasta格式。

3、同源基因搜尋

使用tblastn程式，把已知蛋白質序列和基因組草圖序列建立的本地BLAST資料庫進行比對，注意引數設定(如e-value設為0.00001)，建議輸出格式為5，這樣方便匯入Biopython進行後續解析。

tblastn -query uniprot-taxonomy-fungi.fasta -out fungimax1.blast -db fungi -outfmt 7 -evalue 1e-5 -num_alignments 1 -num_threads 8
tblastn -query uniprot-fungi.fasta -out fungimax10.blast -db fungi -outfmt 7 -evalue 1e-5 -max_target_seqs 10 -num_threads 8
tblastn -query uniprot-taxonomy-fungi.fasta -out fungi.blast.max1.xml -db fungi -outfmt 5 -evalue 1e-5 -max_target_seqs 1 -num_threads 8

4、Genewise建模

呼叫Biopython包，解析blast結果，取出同源蛋白質序列，以及蛋白質序列比對到基因組上位置的基因組序列（上下游再延伸1000bp），放入genewise進行驗證：

圖表 2 Genewise結果
可以發現，確實能夠預測出基因，說明擷取序列正確。

程式設計部分：

基因組不同區域序列的獲取

根據小鼠gff註釋檔案，編寫程式解析小鼠基因組fasta，提取CDS 區域序列、基因內含子序列、基因間序列、外顯子-內含子剪下識別位點附近的序列。並根據解析結果分別獲取如下序列集：
（1）編碼蛋白的基因CDS 區域序列；
（2）基因內含子序列；
（3）基因間序列；
（4）外顯子-內含子剪下識別位點附近的序列【前後各10bp】，分為Intron_5和Intron_3 兩類。
PS：注意擷取CDS序列的時候，需要按照每個transcript合併同一transcript下的所有CDS，CDS序列合併以後應該都是3的倍數。有的CDS序列合併後也不是3的倍數，追蹤編號到Ensemble資料庫中，發現是3’ incomplete之類的，研究不透徹，可以去掉這些CDS。

三聯體密碼子偏好性的分析

對於上述提取的 CDS 序列進行三聯體密碼子出現頻數進行統計（WordSize=3，step=3），終止密碼子不計在內；所有序列統計完成後，計算各個密碼子出現的相對比例。
對於基因內含子序列和基因間序列，按照六個閱讀框進行統計分析（WordSize=3，step=1，互補鏈），終止密碼子不計在內；所有序列統計完成後，計算各個密碼子出現的相對比例。
統計CDS,intron,基因間序列密碼子的出現比例，對密碼子偏好性進行視覺化，繪製出現比例條形圖。

圖表 3 密碼子偏好性
紅色是該密碼子在CDS區出現的比例，綠色和藍色是密碼子在intron和基因間區域出現比例，可以看出來，對於每個密碼子，綠色和藍色高度基本一致，而紅色相差較大，很顯然CDS區域有著特殊的密碼子偏好性。
接下來估算不同樣本之間的相似性度量(Similarity Measurement)，使用歐氏距離(Euclidean Distance)為指標
歐式距離計算公式：
向量形式：

相互比較的型別	CDS與Intron	CDS與基因間區域	Intron與基因間區域
歐氏距離	0.06447049	0.07381776	0.01560867

可以從表中看出來，“CDS與Intron”，“CDS與基因間區域”的歐式距離明顯要高於“Intron與基因間區域”，這說明CDS與Intron相似度很低，CDS與基因間區域相似度也很低，Intron與基因間區域相似度卻很高，這從統計指標上就看出來，CDS區域有著明顯的密碼子偏好性。

剪下識別序列鹼基規律分析

對於擷取的長度為20bp的內含子和外顯子識別序列進行位點鹼基統計，建立序列譜；即讀取Intron_5和Intron_3的序列資訊，並統計出每個位點a/g/c/t的頻數和比例，建立如下的序列譜。

圖表 4 Intron_5序列譜

圖表 5 Intron_3序列譜
對序列譜進行視覺化，繪製訊號強度條形圖：

圖表 6 5’端訊號強度圖

圖表 7 3’端訊號強度圖
綜合上圖，初步可以看出來intron的5’端和3’端分別有GT和AG的剪下識別訊號。
接下來，根據資訊熵計算公式：，計算每個位點的混亂程度，混亂程度越低，保守性越好。然後，對於每個位點的資訊熵，同時用每組最大的資訊熵H減，獲得保守性程度。
對保守性進行視覺化，繪製保守程度條形圖。

圖表 8 5’端保守程度條形圖

圖表 9 3’端保守程度條形圖
將保守性程度較高的那些位點找出來，5’端：第9-16位，3’端：第8-10位，可以初步做出序列譜：
5’端保守性序列：AG/GTAAGT
3’端保守性序列：CAG/
與分子生物學課本上寫的保守性序列一致：
外顯子…AG↓GUAAGT…
內含子…Py10CAG↓…外顯子
接下來，使用MEME線上伺服器進行motif識別，由於只能提交1000條序列，因此不一定很準，可是對於intron5’端剪下識別位點，依舊精確地找出了motif

圖表 10 MEME intron5’端motif
對於intron3’端剪下識別位點，也能找出CAG的motif

圖表 11 MEME intron3’端motif

剪下識別序列的計算鑑別

根據intron3’端和5’端20鹼基的序列譜，分別建立打分模型。評分計算方式：
取出20bp的序列，依次找出鹼基在intron5’端中的比例，進行累乘，得出打分值。例如某一20bp序列第一位是A，找出intron5’端20bp第一位鹼基是A的概率P1A，第二位是G，找出intron5’端20bp第二位鹼基是G的概率P2G，第三位是T，找出intron5’端20bp第三位鹼基是T的概率P3T，……第二十位是C，找出intron5’端20bp第二十位鹼基是C的概率P20C，然後，依次累乘P1A*P2G*P3T*…*P20C，得出這一序列的打分值。
如果取出的這一序列是intron5’，對應於intron5’建立的模型；或者這一序列是intron3’，對應於intron3’建立的模型，則求出的是陽性經驗分值；如果模型與序列不對應，則為陰性經驗分值。
接下來繪製打分密度圖：

圖表 12 得分密度圖
- Score3為輸入intron3’剪下位點的序列到對應intron3’模型得出的值；
- Score5為輸入intron5’剪下位點的序列到對應intron5’模型得出的值
- Intronscore3為輸入intron序列到intron3’模型得出的值；
- Intronscore5為輸入intron序列到intron5’模型得出的值；
- Intergenescore3為輸入基因間序列到intron3’模型得出的值；
- Intergenescore5為輸入基因間序列到intron5’模型得出的值；

可以明顯從圖中看出，模型自身對應序列打分會更高；如果關鍵特徵全都對應上，是能用此模型來分類，但是從圖中也能看出，單純根據此模型進行打分，3’剪下位點得分大部分與陰性序列相近，用這個模型不易將3’剪下位點識別出來。

人工神經網路(ANN)

由於之前打分模型過於簡單，不易將3’剪下位點識別出來，下面構建人工神經網路（Artificial Neural Network，即ANN），對序列進行識別分類。由於資料量過大會增加運算時長，因此3’剪下位點序列，5’剪下位點序列，內含子序列，基因間序列一共提取5w條，進行分類。為了方便起見，構造的神經網路為單層的，隱藏的神經元（Hidden Unit）數量分別嘗試3,6,12,18個，學習速率（Learning Rate）分別嘗試了0.003, 0.01, 0.03,0.1，迭代次數為200，尋找最優的神經網路。

圖表 13 每輪迭代的誤差圖
上圖為誤差圖，顯示不同神經元和不同學習速率下，誤差隨迭代次數下降的圖，橫座標為迭代的次數，縱座標為誤差大小，紅線為測試集的誤差，黑線為訓練集的誤差。明顯，黑線一直在紅線之下，即作為生成模型的訓練集分類效果會更好，
這裡還需要解釋一下隱藏的神經元（Hidden Unit）數量和學習速率（Learning Rate）對最終結果的影響。隱藏神經元越多，越容易擬合出複雜的模型，但是在資料量不夠的情況下容易“過擬合”；而學習速率越大，模型的收斂速度會在更少的迭代次數後收斂，但是如果學習速率過大，有可能錯過全域性最小值而無法收斂，就比如上圖中最後一行，學習速率都是各組學習速率最大值0.1，測試集的誤差一直在波動而難以完全收斂。

接下來計算訓練集，測試集的誤差，選取最好的神經網路

圖表 14 誤差計算
可以發現最好的是隱藏層12個神經元，學習速率0.1的那個神經網路。
輸出該模型的基本資訊：
SNNS network definition file V1.4-3D
network name : RSNNS_untitled
no. of units : 36
no. of connections : 288
learning function : Std_Backpropagation
update function: Topological_Order

圖表 15 模型詳細資訊
try my best地來解釋一下，其實結果就是建立了這樣的一個模型：

圖表 16 模型簡圖
20個輸入節點（二十維的評價集，每個節點代表原來20bp序列上的一個鹼基），自己建立的12個隱藏節點，以及輸出的4個節點（四維的結果集，分類出intron，intergene，3’剪下位點，5’剪下位點）
所以總共是36個節點（no. of units: 36）
20*12+4*12=288條連線（no. of connections: 288）。

接下來，使用測試集進行測試，得出預測結果統計表

可以看出來，5’剪下位點、3’剪下位點以及內含子序列識別效果很棒，而有大量基因間序列被識別成為了內含子，因為本身“基因間序列”和“內含子”特徵就不是很明顯，這一結果也是可以接受的。