北邮生物信息基础大作业.docx

北京邮电大学信息与通信工程学院生物信息基础（2014）课程大作业一、 疾病自动诊断问题1、 题目分析根据题目要求，我们需要设计一套计算机自动筛选方案，目的是通过ebov-10得到的十项指标，将疑似患者中的埃博拉病毒可能携带者筛选出来。目前，我们已有的训练集是经过专家筛选后的100例疑似患者的十项指标。其中，20位为埃博拉病毒可能携带者，80位已被排除患病可能性。综上所述，我认为该问题为一个监督下的模式分类问题，两个分类指标为“埃博拉病毒可能携带者”、“非埃博拉病毒可能携带者”，观测向量为经过ebov-10得到的疑似患者的十项指标：x=(x1,x2,x3,x4,x5,x6,x7,x8,x9,x10)t2、 建模流程信息获取与预处理部分，在之前的ebov-10检查中已经较为详细的给出，这一步不再设计；特征选择与提取部分，我计划使用主成分分析方法，通过对十项标准进行线性组合，可以得到更能够体现类间信息的新的一组观测向量；分类器设计采用fisher线性判别分析+最小错误率贝叶斯决策的方法。3、 主成分分析主成分分析的思想是从一组特征中计算出一组重要性按从大到小顺序排列的的新特征，它们是原有特征的线性组合，并且之间是互不相关的。设为x的协方差矩阵，求解出矩阵的各个特征值与特征向量，则特征值最大的特征向量1，为数据集的最佳投影方向。由此方向投影，可获得最大的投影数据的方差。按照这个思路依次找到次大的，第三的特征值对应的特征向量，它们就是次优的，第三优的投影方向。我们可以只提取重要性占前k%个主成分：(i=1ki)/(i=1ni)k我们把原始数据集按这些方向投影，得到的就是降维后的观测向量。选择较少的主成分来表示数据，不但可以用作特征的降维，还可以用来消除特征中的噪声。4、 fisher线性判别分析fisher线性判别分析的基本思想是：将所有的样本投影到一个方向上，然后在这个一维空间确定一个阈值。选择最优的投影方向应该使得各个样本点的类内方差最小，类间方差最大。我们定义类内离散度矩阵sw，类间离散度矩阵sb，投影向量，准则的目标函数为：maxj()=sbsw=wtsbwwtsww这是一个约束条件下的极值问题，我们可以利用拉格朗日乘子法求解。拉格朗日函数：lw,=wtsbw-wtsww-c上式在极值点处，应该满足对w的偏导数等于零。可以解得：w*=sw-1(m1-m2)式中，w*为fisher线性判别准则下的最优投影方向。m1和m2为两类的类均值向量。由于fisher线性判别分析不对样本的分布做任何假设，当样本维数较高样本数也较多的时候，投影到一维空间后样本接近正态分布。这时可以在一维空间中使用正态分布拟合样本，再使用上面提到的最小错误率贝叶斯决策，往往会有很好的效果。5、 最小错误率贝叶斯决策最小错误率贝叶斯决策的决策规则为：如果p(1|x) p(2|x)，则x1，否则x2。其中，后验概率可以使用贝叶斯公式求得：pix=p(x|i)p(i)p(x)=p(x|i)p(i)j=12p(x|j)p(j)但是，将一个疑似患者判为病毒携带者和排除其患病可能性，其代价（损失）是不一样的。对此，我认为可以采用最小风险的贝叶斯决策。这种决策方法中，决策表是需要人为确定的，需要认真分析研究问题的内在特点和分类目的，与疾病防控领域的专家共同决策，设计出适当的决策表。具体的决策步骤：（1） 利用贝叶斯公式计算后验概率（2） 利用决策表，计算条件风险rix=j=1c(i|j)p(j|x)（3） 决策：在各种决策中选择风险最小的决策，即=argmini=1kr(i|x)二、 病毒变异与否的判断1、 问题分析针对病毒变异与否的判断问题，我认为应该采样合适数目的埃博拉病毒dna序列（可以是关键部分的基因），然后使用多序列比对的方案，将多条序列对齐，就可以方便的定位出序列中碱基对的差异位置，从而评估病毒是否发生了变异。由于整条序列做多序列比对，可能时间复杂度较大，最后得到的变异位点信息也和我们想要的相差较大。我们可以选择和病毒的致病能力密切相关的几个基因，作为多序列分析的原材料。2、 求解流程3、 多序列比对下面简要介绍多序列比对的实现方法。由于病毒的dna序列较长，我们采取星形比对的方案。星形比对的基本思想是：在给定的若干序列中，选择一个核心序列，通过该序列与其它序列的两两比对，形成所有序列的多重比对，从而使得该多重比对在核心序列和任何一个其它序列方向的投影是最优的两两比对。下面给出星形比对的基本过程：1. 选择核心序列2. 计算与核心序列的两两比对3. 逐对聚合两两比对的结果，获得多重比对选定一个核心序列，把多重比对转化为k个两两比对聚集过程。从某一个两两比对开始，比如sc和s1，然后逐步加上其他的两两比对。在这个过程中，逐步增加sc中的空位字符，以适应其他的比对，但不删除sc中已经存在的空位字符。选择核心序列的方法为：尝试将每一个序列分别作为核心序列，进行星形多重序列比对，取比对结果最好的一个。在上面提到的星形比对中，最基本的核心是两两比对，我们使用的两两比对的方法为：全局最优序列比对的动态规划求解算法。给出求解过程： 初始化dp辅助矩阵 根据状态转移方程递归计算dp辅助矩阵a 确定最优路径，即对应于最优比对矩阵更新策略：ai,j=ai,j-1+p(-,tj)ai-1,j-1+p(si,tj)ai-1,j+p(si,-)4、 寻找变异位点需要的序列经过了比对，已经对齐。我们可以定义一个阈值，当碱基序列中连续出现的变异碱基个数超过了这个阈值后，我们即可认定这个序列是变异序列，与原序列的差异较大。阈值的选择要与疾病防控领域的专家一起协定，这样才能够准确的发现基因出现的变异情况。5、 可靠性检验下面，我们对上述方法作了可靠性的检验：由于我们缺乏埃博拉病毒的必须dna数据，所以我们选择了来自其他生物的7组dna数据对上述多序列比对的方案作了验证。这七组dna分别是：大鳞泥鳅雌激素受体mrna，斑马鱼雌激素受体b mrna，海猪鱼雌激素受体mrna，鳗鲡雌激素受体mrna，人类14号染色体部分片段，须鱲雌激素受体mrna，拟鲤erb基因雌激素受体 mrna。（注：以上七条数据来自nucleotide数据库）fasta格式的各个序列头标记：gi|145308317|gb|ef530592.1| paramisgurnus dabryanus estrogen receptor beta mrna, partial cdsgi|23466358|gb|af349413.3| danio rerio estrogen receptor beta b mrna, complete cdsgi|32186925|gb|ay305027.1| halichoeres tenuispinis estrogen receptor beta mrna, complete cdsgi|2073112|dbj|ab003356.1| anguilla japonica mrna for estrogen receptor, complete cdsgi|89037528|ref|nw_925528.1| homo sapiens chromosome 14 genomic contig, alternate assembly hs_celera 211000035825320, whole genome shotgun sequencegi|30962102|emb|aj314602.1| candidia barbatus mrna for putative estrogen receptorgi|61097789|dbj|ab190290.1| rutilus rutilus erb mrna for estrogen receptor beta, complete cds 我们利用clustalw2网站，对以上7条序列作了在线多序列比对，比对的结果如下（篇幅所限，这里只给出节选）：gi|30962102|emb|aj314602.1| ctgcccgcc-tcaagtggcctacagc-gaaacacactcacacactgcct 429gi|61097789|dbj|ab190290.1| ctgcccacc-tccgctggcctacagc-gaaacattctcacacactgcct 376gi|23466358|gb|af349413.3| ctgcccgcc-tccgctggcctacagc-gaaacacgttcacacagcgcct 433gi|32186925|gb|ay305027.1| ctgcccaca-gcctctgggctacaat-gaatccggcttacacgcaccct 526gi|2073112|dbj|ab003356.1| cttccagca-gcccctggtgtacaga-gagcccgcc-cactccccgt 677gi|145308317|gb|ef530592.1| -gi|89037528|ref|nw_925528.1| tagataacaagtccggtttcctgcagaagaggcctcatcgccagcaccct 433gi|30962102|emb|aj314602.1| gtgtctggtcatgtgagggatgcaaggcttttttca-aacggagca- 651gi|61097789|dbj|ab190290.1| gtgtctggtcatgtgaggggtgcaaggctttcttca-aacggagca- 598gi|23466358|gb|af349413.3| gtgtctggtcatgtgaagggtgcaaggctttcttca-agcgtagca- 655gi|32186925|gb|ay305027.1| gtgtgtggtcctgcgagggctgtaaagcatttttca-agaggagta- 757gi|2073112|dbj|ab003356.1| gggtgtggtcctgcgaaggctgcaaggccttcttca-agaggagca- 893gi|145308317|gb|ef530592.1| gggtgtggtcatgcgaggggtgcaaggctttcttca-aacgaagca- 55gi|89037528|ref|nw_925528.1| aaaaat-tcctg-gaaaac-ccagacttatctcagacagaggagaaac 662gi|30962102|emb|aj314602.1| cctgagagagccagtgaagaagccat-acacg-gaggc-cagcatgatga 1059gi|61097789|dbj|ab190290.1| cctgaaagagccagtgaaaaagccgt-acact-gaggc-cagcatgatga 1006gi|23466358|gb|af349413.3| cctgagagagccggtgaaaaagccat-acact-gaggc-tagcatgatga 1063gi|32186925|gb|ay305027.1| cctcatggaggagcagaagaagcctt-ttacc-gaggc-cagcatgatga 1174gi|2073112|dbj|ab003356.1| cctcatgaaggagctgaagaagccct-tcacc-gagga-cagcatgatga 1268gi|145308317|gb|ef530592.1| cctgaaagagcaggtgcagaagccgt-acact-gaggc-cagcatgatga 457gi|89037528|ref|nw_925528.1| tttcagacagtttct-acctgtatcacccaaggtgcagtttgatgt 1056 通过以上多序列的比对结果，我们可以轻松地发现多条序列中，发生了变异的部分。由于序列已经对齐，发生了变异的部分必然是空位、插入或者替换中的一种。我们设计算法，遍历一遍对齐后的各条序列，即可轻松的找到变异位点。三、 基因编码区域识别1、 问题建模由于已经给出了基因编码区域和基因非编码区域的片段（训练样本），我考虑使用一阶马氏链来判别两个目标序列的区域，其中，一阶马氏链是根据碱基的排列顺序做转移的。选择一阶马氏链碱基排序的原因有二：一是实验的样本太少，我考虑了对于氨基酸密码子做转移概率矩阵，但是计算得到的矩阵很多都是0元素，这对最终的计算结果影响很大，我也尝试了二阶和更高阶次的马氏链，同样由于训练样本太少，转移概率中概率为0的点太多，故不采纳高阶方案；二是在题目中，明确给出了该病毒rna具有特定的排列顺序的条件，所以不对密码子做转移概率计算。2、 参数估计与计算过程首先，根据给出的基因编码区域和基因非编码区域的片段，可以算出两个区域片段的各个碱基对分布概率。这里利用大数定理，使用频率逼近概率。然后，计算相邻两个碱基对出现的频率，并以此作为一阶马氏链的转移概率。最后，由两个目标序列的排列，我们可以利用公式：ps=p(s1)i=1np(si+1|si)来计算序列s出现的概率。3、 实验结果在编码区中，四个碱基的分布概率：符号augc概率0.3846 0.15380.30770.1538在非编码区中，四个碱基的分布概率：符号augc概率0.1538 0.30770.07690.4615在编码区中，各个碱基的转移概率矩阵：augca0.3913 0.30430.30430u0.4444 00.55560g0.2105 0.10530.15790.5263c0.5556 00.44440在非编码区中，各个碱基的转移概率矩阵：augca0 001.0000u0 00.06250.9375g0 1.000000c0.3448 0.51720.10340.0345计算得，s1序列在编码区出现的概率为1.6332e-06，在非编码区出现的概率为0，s2序列在编码区出现的概率为0，在非编码区出现的概率为6.4733e-04。所以，我们认为s1序列属于编码区，s2序列属于非编码区。在判定过程中，我们发现了两个概率为0的现象，一个是s1在非编码区的概率，另一个是s2在编码区的概率。出现这两个0概率的原因为，由于训练集过小，无法保证所有的碱基组合都出现过，这样就有部分转移概率为0。注：本部分计算所用的matlab源程序已在附录中给出。四、 课堂内容回顾在大三上半学期，我有幸选修了生物信息基础这门专业选修课。选修这门课的原因主要有二，一是我希望未来可以在模式识别领域继续深造，希望可以考取我校模式识别实验室的研究生，所以希望通过这门选修课，让我对模式识别的基本理论和基本算法有所了解；二是去年四月，我参加了2014深圳杯大学生数学建模夏令营，竞赛题目中有一道关于基因组测序的题目（b题），引起了我很大的兴趣。通过这学期生物信息基础课程的学习，我对于生物信息处理有了很深的理解。我们的课程从生命的演化与中心法则讲起，包含生物信息数据库、序列分析、基因组学与基