北邮生物信息基础期末作业.docx
《北邮生物信息基础期末作业.docx》由会员分享，可在线阅读，更多相关《北邮生物信息基础期末作业.docx（13页珍藏版）》请在冰豆网上搜索。
北邮生物信息基础期末作业
•、疾病自动诊断问题
1、设计思路2
2、基于K-L变换的主成分分析法2
3、近邻法判别分析3
二、病毒变异与否的判断4
1、设计思路4
2、基本步骤4
3、星型比对算法4
4、两两比对算法5
5、寻找变异位点5
三、基因编码区域识別5
1、设计思路5
2、问题建模5
3、参数估计与计算过程6
4、实验结果7
四、隐马尔科夫模型的求解7
1、设计思路8
2、问题建模8
3、参数估计与计算过程9
4、计算结果9
五、课堂内容回顾10
《生物信息基础》期末作业
一、疾病自动诊断问题
1、设计思路
根据题目，我们获得80份腹泻类型为细繭型的常规化验结果数据，以及20份腹泻类型为病毒型的化验结果数据，基丁•以上100份数据中的20项数值型结果，我们将其作为分类器训练集，设计出介理的算法将某常规化验结果分类，判断其腹泻类型。
该计算机自动诊断方案包括以下部分：
-4分类器设计.训练
信息获取与预处理——►特征选择与提取
f分类决策（识别）——►结果输岀
其中.信息获取与预处理，即为通过常规化验，采集化验所得的？
0项数据结果X=（Xj,x2,X3,X4,,X17,x18,x19,x20）T»并将其作为观测向屋。
特=征选择与提取采用通过K・L变换实现主成分分析法，分类器设计及决策采用近邻法，基于以上流程，将输出的结果作为判别结果。
2、基于K-L变换的主成分分析法
主成分分析的基本思想是进行特征降维的变换，虽然不能完全地表示原有的对彖,能最总会有损失，但是希望找到一种能最最为集中的变换方法使得能最损失最小，即根据己知向量，将其做变换得到一组互不相关新的特征，对新特征向量进行排序，舍掉瑕小的特征，从而使得损失的能量最小。
1）K・L变换的基本算法
原始输入为x=（x】>X?
X3,X4,...,X】7>X]8,X】9,XgoF，变换后特征为y,变换矩阵（线性变换）为A,则y="x,变换后的相关矩阵则&.="心4。
若不同的y特征互不相关，则对R>.进行分析，以R*的特征向量作为A的列，则满足：
Ry=AtRxA=[aL,a2,...,an]TRx[aL,a2,...,an]=⑺丄心，•••宀卩⑺丄“血“•••人aj=A,其中A为对角矩阵，对角线元素为九,耳…几，以此来达到变换后特征不相关的目的。
2）主成分分析的基本算法
通过K・L变换时特征降维，原冇N维，只保留m维，去^ym+1...Zv,若希望和原来的表示方法差别最小，即E[||x'-x||2]最小，”表示[力在原空间中对应的表示方法。
通过计算不难得出如下公式：
N
e[||x'7『]=》人
匸M+丄
将人从人到小排序，则保留m维后的结果即为损失的能呈垠小。
=KMMGorWT5AMOnttCC**4UMCAno«S
3）简要分析
通过主成分分析可以实现将所得特征降维的目的，同时还可以消除特征中误差的干扰，使得特征提取的结果更为准确，有利于之后分类器的设计及决策。
3、近邻法判别分析
1）基本算法
近邻法的基本思想为：
将于测试样本最近邻样本的类别作为决策的方法。
对于一个c类别的问题，每类有M个样本，二…,c,则第®类的判别函数为：
^（x）=miiifc||x-x^||,k=
其中才为第3类的第k个样本。
决策函数为:
g）M=min^t（x）i.i=1,...,C，则决策xGa）jo
由题目得，该诊断方案需要判断的类别右两类，即细菌型腹泻和病毒型腹泻。
己知细菌型腹泻有80个样本作为训练集，病毒型腹泻有20个样本作为训练集，以此为分类基础设计出一个基于近邻法的分类器。
2）错误率分析
己知待分类样本X,其敲近邻样本为,，已经有N个已知类别的训练样本，而，是其中之一，则特征空间中x点的错分概率：
P（e\X）=jP（e\X,X'）p{X^\X）dJC
如果样本X的两类别后验概率分别为P（5|x）ljP（32|x），那么对X值，在NT8条件下，发生错误决策的概率为：
辄IX）“-0⑷IX）-0（兮I£
而在这*件卞的平均错误率为：
P=粧=lim|p（X）dX叮烈％（€|品£（尉茫=卩-£冲网|卫"的戏
设贝叶斯错误率为严•比较两者错误率，得出如下关系:
PS02■才）
P
3）改进：
K・近邻法
瑕近邻法町以扩展成找测试样本的k个最近样本作为决策依据。
基本规则为:
在所有N个样本中找到与测试样本的k个最近邻者，其中各类别所占个数为心,1=1,…,C,则决策规划为：
kj（X）=mni{k（（X）,i=l.,C,
则做出决策：
xea）J9
二、病毒变异与否的判断
1、设计思路
判断病毒是否发生变异涉及序列比对的问题。
已知15条病毒型腹泻患者的病廷核酸序列，将待测序列与样本集进行多序列比对，判断病毒是否发生了变异。
当待判断序列与样本序列特征差异超过一定阈值时，判断该病毒发生了变异，同时定义另一阈值,当碱基对连续错误超过该阈值时，判断起始位点处发生了变异。
内为基因突变等现象的存在，特征差异及碱基对的阈值设定需要专家进行专业评估后给出。
方案中多序列比对采用星型比对的方法，通过星型比对获得15条对齐的序列，从而总结出该样本集的序列特征，之后将待判断序列与该特征进行对比，将对比的结果与设定的阈值进行比较，从而判断是否变异，如果该病毒发生了变异，巫新遍历对齐序列，记录碱基对连续出现变异的数目，也与设定的阈值进行比较，从而找出变异位点。
2、基本步骤
3.星型比对算法
目前所用的人部分序列多重比对算法需要交多重序列比对转换为两两序列比对，再将两两序列比对组合起来，最终形成完整的参序列比对，这种方法称为渐进法，星型比对是其中的一种。
星形比对的基本思路是：
在给定的若干序列中，选择一个核心序列，通过该序列与其它序列的两两比对，形成所冇序列的多重比对，从而使得该多重比对在核心序列和任何一个其它序列方向的投影是最优的两两比对。
星型比对的基本过程为：
1）选定一个核心序列；
2）计算与核心序列的两两比对：
3）逐对聚合比对结果，获得多重比对。
聚集过程从某一个两两比对开始，比如Sc和S],然后逐步加上其他的两两比对。
在这个过程中，逐步增加Sc中的空位字符，以适应具他的比对，但不删除Sc中已经存在的空位字符。
选择核心序列的方法为：
尝试将每一个序列分别作为核心序列，进行星形多重序列比对，取比对结果最好的一个。
4、两两比对算法
在上面提到的星形比对中，故基本的核心是两两比对，我们使用的两两比对的方法为：
全局最优序列比对的动态规划求解算法。
给出求解过程：
1）初始化DP辅助矩阵；
2）根据状态转移方程递归计算DP辅助矩阵A：
3）确定最优路径，即对应于最优比对。
矩阵更新策略：
M（iJ—l）+p（—,期）
Aftj）=--1J-1）+p^s[i],t[j]）•
、做一1J）+P（s[i]・-），
5、寻找变异位点
待测序列在经过比对之后，己实现对齐。
在此基础之上，我们只需再重新进行一次遍历，记录碱基序列中连续出现变异的碱基个数，并与专家给出的阈值进行比较，从而找出变异的起始位点。
三、基因编码区域识别
1>设计思路
根据已经给出的编码区域与非编码区域的碱基序列样本，结合题目中提及的该病毒RNA序列中基因编码区域的碱基排列顺序具有特定规律，考虑利用马尔科夫链分别来计算所给出的序列片段出现在编码区和非编码区的概率，通过比较概率的人小來判断其足否为编码区域。
利用马尔科夫模型计算概率时，需要统计出各碱基的起始概率及概率转移矩阵，因为题冃中所给的序列长度较短，考虑一阶马氏链的情况，各概率矩阵在参数估计中给出。
问题建模
利用号尔科大链模型求解时，需要考虑一个只有多个状态的系统S,令…4为一系列在各个时刻系统状态的变量，即状态链。
对于一阶马氏链，系统在时间步t处于哪个状态，仅与时间t-1时系统所处的状态何关，P（O1…0n）=P（0jP（02|01）P（03|02）・“P（0nl4』，即系统未來的状态仅依赖于当前状态，马尔科夫链完全决定于初始分布和状态转换概率。
计算儿到勺的转移概率，需要统计令到坷的转移次数，然后除以所有以人为起始状态转变为4k的转移次数：
P⑷坷
k
其中从一坷表示从状态儿转移到状态％的次数。
第一个出现的碱基和最后一个出现的碱基需要单独统计：
P⑷=P⑷[]）=T”/》尢―如
k
P（%）=P（[]\Aj）=G一[]/》N一d
k
即通过上式从训练样本序列中统计碱基儿出现在第一位的频率和碱某勺出现在最后一位的频率。
3.参数估计与计算过程
首先，根据给出的基因编码区域和基因非编码区域的片段，利用频率逼近概率，町以算出两个区域片段的各个碱基对的起始概率分布及各自的一步转移概率矩阵。
利用如下公式计算给出序列在编码区及非编码区出现的概率：
n
P（$|（G）=P（sJ[》P6+i|sJ]P（sJ
4=1
利用贝叶斯公式：
P（S=M4jSeC）P（SeC）
工］p（Sf4isg）p（sy）
如果P（CdS）>P©|S）,则序列S属于类别5否则，属于类别由于计算过程中，贝叶斯公式中编码区和非编码区的分母相同，而P（s（ect）未知，假定p（s（eCi）=p（s（eC2>05,则结果比较只需比较先验概率大小即可。
在編码区中，四个碱基为起始状态的分布概率:
符号
A
U
G
C
概率
3/5
1/5
1/5
0
在非编码区中，四个碱基为起始状态的分布概率：
符号
A
u
G
C
概率
0
0
1/5
4/5
在编码区中，四个碱基为结乡
E状态的分布概率：
符号
A
u
GC
1概率
2/5
1/5
1/5
1/5
在非编码区中，四个碱基为起始状态的分布概率:
符号
A
U
G
C
概率
0
4/5
0
1/5
在编码区中，各个碱基的转移概率矩阵:
A
u
G
C
A
9/23
7/23
7/23
0
u
4/9
0
5/9
0
G
4/19
2/19
3/19
10/19
c
5/9
0
4/9
0
在非编码区中，各个碱基的转移概率矩阵:
A
u
G
c
A
0
0
0
1
U
0
0
1/16
15/16
G
0
1
0
0
c
10/29
15/29
3/29
1/29
4、实验结果
所给待测序列片段为：
[SxJGAAUGUAGCGCAA
[S2]CUCUCUCACACGU
计算S1,其在编码区的概率为：
p⑸G）=G）.（訐（分（分®（訪0*&）*盼（9-備
*（I）*（^）*（D=4246442*10-7
氏在非编码区的概率为：
P（SJC2）=（訴（0）=0
计算S2,其在编码区的概率为：
P^G）=0
其在非编码区的概率为：
p（sm）=G）.禽倚g）.S）.（3.（存匍畑g）.
（1）.
（^）*
（1）*（^）=8.97637171*10-4
根据上述计算结果，得p（s1|c1）>p（51|c2）,PCsjq）o.5的情况下,使得p（G|sj>pgisj,Pfcj^）序列片段S]为编码区域片段，序列片段S2为非编码区域片段。
4.隐马尔科夫模型的求
1、设计思路
根据已知题目，我们需要通过町观测的字符C、G來判断生成其的状态Si、S2、S3,题冃中已经给出需要应用隐马尔科夫的相关知识，计算生成概率及判断最町能的状态序列分别对应其中的评估问题和识别问题，利用所给概率带入公式计算即町。
因为题目中所给的状态转移矩阵中包含一列空状态，为使转移概率矩阵中每行概率和为一，可采用两种方法：
1）不考虑空状态，对每行概率进行归一化。
以第t步时状态为S.则将第t+1步转移至£的概率由03归一为0.3/（0.3+0.2+0.3）,然后将表中将用到的概率均进行归一，采用归一后的概率进行计算。
2）不进行归一，考虑将空状态引入状态集，定义为So,同时需对明字符生成概率进行修正，So生成C、G的概率均为0,即空状态只是表示初始和结尾，不会生成任何字符集。
S】、S2、S3均会生成明字符，即不生成字符的概率为0。
权衡两种方法，归一化之后计算繁琐，且相对不严谨，决定采用第二种方式进行计算。
一阶离散隐马尔科夫模型（HMN）是一个关于时间的序列的随机生成模型，其包含以下要素：
1）有限的隐状态合集：
S={S]…・,Sn};
2）离散明字符集：
V={"…,VM}；
3）初始状态概率向量：
"=5丄，…,ttn};
4）状态转移概率矩阵：
A={atjbj=1…・,N}：
5）明字符生成概率矩阵：
B={bjkj=1，…,N,k=1M}o
HMM记作入=（S,V,ir,A,B）或A=（11,AB）。
HMM将以卜两个序列联系起来：
1）由离散隐状态组成的状态序列：
Q={如，…,他},每个仇6S是一个状态，由初始状态概率及状态转移概率（1V,A）所决定：
2）由明字符组成的观察序列：
0={。
】…0”，每个otev是一个状态，由状态序列及各状态的明字符生成概率（Q,B）所决定。
分析题目，利用采用的第二种方法获得以下信息：
1）有限的隐状态集合S={S0,SpS2,S3}：
2）离散的明字符集合V={",«，◎}，其中”代表生成C,冬代表生成G,%与%相似,为“[]”空状态；
3）初始状态概率向量tt={no，nlfn2>n3},其中心=0,心=0.2,兀2=
0.39心=0.5：
4）状态转移概率矩阵A各元素值如卜表:
So
S\
S2
S3
S。
0
0.2
03
0.5
0.2
03
0.2
0.3
$2
0.3
01
04
0.2
S3
0.4
01
0.2
0.3
5）明字符生成矩阵B各元素值如下表:
[]
牛.成C
生成G
So
1
0
0
Si
0
0.4
06
S2
0
0.3
07
S3
0
06
04
6）O={yyVLfV2fVyV0}9即为CtCtGtCt[]。
根据上述己知信息•计算该观测序列在所给定的HMM卜的生成概率及生成该观测
序列的最佳路径（最M能的状态序列）。
3、参数估计与计算过程
1）评估问题：
计算o出现的概率
定义前向变量为“在时间步t,得到tZ前所右的明符号序列，且时间步t的状态是甲这-事件的概率，贝Ija（ti）=P（ox，…，otfqt=SJA）o算法过程包括以下三步：
a）初始化：
a（l.i）=7r（i）b（i,oJ：
b）递归：
a（t+l.i）=E；=i«（ti）a（tj）]^G'»°t+i）?
c）终结：
P（O|A）=S；=i«（T,i）o
2）识别问题：
判断0是由哪个HMM产生
定义阻i）为t时间步沿状态序列五…qtRqt=St产生出。
丄…・,的最大概率,即8（t,i）=maxg”./.]P=S[,o「入）。
釆用Viterbi算法，具体步骤如卜•:
a）初始化：
&（[）=7r（i）b（i,oJ；
b）递归：
6t（J）=max0a[&7（i）a@]bj（oJ:
c）
终结：
P*=maxg勿[8r（i）]：
4.计算结果
1）评估问题：
计算o出现的概率将每步计算所得的a（ti）绘制如卜表:
t=l
t=2
t=3
t=4
t=5
1=0
0
0
0
0
0004990032
1=1
008
00252
0011304
000292608
0
1=2
009
00336
0.024024
000447408
0
1=3
03
00792
0.015216
0.00765648
0
根据上述表格，得岀观测序列为CCGC的概率为:
3
0.00490032+0+0+0=0.004990032
P（0|入）=》（t（5,i）=
i=Q
2）识别问题：
判断O是由哪个HMM产生将每步计算所得的&ti）绘制如下表：
t=l
t=2
t=3
t=4
t=5（*10-5）
1=0
0
0
0
0
46656
1=1
0.08
00012
0.00324
0.0003888
0
1=2
009
00018
000756
0.0009072
0
1=3
03
0054
0.00648
0.0011664
0
根据表格回溯，将各口计算中使用的6（t-1」）绘制如卜表:
5（ti）
t=l
t=2
t=3
t=4
t=5
1=0
初始化
0
0
0
8（4,3）
1=1
初始化
5（1,3）
8（2,3）
5（3,1）
0
1=2
初始化
8（1,3）
8（2,3）
8（3,2）
0
1=3
初始化
8（1,3）
8（2,3）
8（3,3）
0
注：
恒为0时不计算Mti）
综合以上两表.=5时概率最人的为5（5,3），即取步长为5时对应的状态为小空状
L仆丫成胡F?
；■,卜［制心获得'bI；为4：
M（爪汕丿状态为S3,沁‘灿能h水心乂对为6（3,3）,5（2,3）,8（1,3）o所以最可能的状态路径为:
S3->S3->S3->^3->［］»其中［］代表末尾空状态。
五、课堂内容回顾
大三上学期的专业选修课中，《生物信息基础》与《模式识别》两门课程有很大的共通之处，本门课程中的对基因序列、待测样本等进行分析、比对时使用的方法，大多数都在模式识别踝上也有所提及。
在我看來，两门课相辅相成，融会贯通。
但是，本门课程内容结介了生物学的相关知识，更显生动有趣。
高中在学习生物的过程中，就对基因的表达很感兴趣.通过本门课程的学习获益匪浅。
特别是在比对两个基因序列的时候，通过矩阵的形式计算打分函数，从而获得放优比対，将复杂的问题简单化，按照方法按部就班就町以获得想耍的结果，这一内容是之前没冇接触过的。
此外，隐马尔科夫模型的相关过程中，虽然状态我们无法观测，但是可以通过对其生成的明字符集进行观测，通过计算概率获得我们想要知道的信息，并且利用了前向、后向、VitnI等算法简化时间复杂度，使得我对马尔科夫模型的了解刃为深入，接触了巫深戻次的应用范阳。
这两部分内容小我在本学期滦程过程中较为感兴趣的部分。
本节课程的教学过程中，大多数时间都是在介绍抽象的算法，个人感觉大最时间介绍算法对于我们的理解帮助不人，因为我们在生物信息学方面的知识还相对薄弱，対于一些复杂的算法难以掌握，建议老师在今后的教学过程中，可以结合更多的例子來介绍算法，便于学生的理解，也有利于调动大家的学习兴趣。
例如在模型时，结介了赌场骰子的问题，就使我对模型印象深刻，在之后计算作业及试卷中的题冃时,起到了爭半功倍的效果。
个人建议在之后的教学过程中，老师町以尽量把声音放人一些,让后排的同学能听得更加清楚。
总而言之，一学期的学习开启了我对生物信息学认知的大门，老师在课堂上除了给我们讲授知识，还会穿插一些人生桝理和萌萌哒冷笑话，课堂效果较为活跃，在今后的生活中，我将继续学习这方面的内容，感谢李老师这一学期的循循善诱，让我获益良多！
升级会员 