真核生物基因组的基因分析和预测
一、摘要
- 加深基因预测基本原理的理解（如
- 密码子的偏好性、内含子外显子剪切识别序列等）；
- 了解同源基因预测的意义所在；
- 熟悉已有的基因预测的使用（如GenScan、GeneWise等）；
二、材料和方法
1、硬件平台
处理器:Intel（R） Core（TM）i7-4710MQ CPU @ 2.50GHz
安装内存(RAM)：16.0GB
2、系统平台
Windows 8.1、Ubuntu
3、软件平台
- BLAST
- GenScan
- GeneWise
- Python3.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 84、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，然后，依次累乘P1AP2GP3T*…*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）
2012+412=288条连线（no. of connections: 288）。接下来，使用测试集进行测试，得出预测结果统计表
可以看出来，5’剪切位点、3’剪切位点以及内含子序列识别效果很棒，而有大量基因间序列被识别成为了内含子，因为本身“基因间序列”和“内含子”特征就不是很明显，这一结果也是可以接受的。
