problem_set_2
7.36/7.91/BE.490
家庭作业 2
上交时间：3月11号下午 1:00
注意：请查阅课程网站获取如何提交程序设计方面问题的电子版。另外，在本次作业结尾有一些关于编程方面的有用的提示。
1．古基因组学 I ——BLAST检索核酸序列
你的实验室已经找到一种从恐龙化石中提取出DNA序列的新方法，现在为了研究恐龙的进化史和鸟类的起源问题，你们从翼龙的化石中提取出了DNA序列。
该序列保存在dino1.fa文件中。
a. 在NR数据库中使用BLASTN进行检索，在使用缺省参数设置时所得到的检索结果
的第一条序列是什么？它看上去正确吗？为什么？
b. 将错配的罚分从缺省值（-3）改成-1，再检索一次，结果的第一条序列是什么？结
果中的E-value值是多少？跟a的结果相比检索质量上是否有不同？降低错配罚分
的值是否得到了你期望的结果？
c. 如何使用BLASTX检索？使用BLASTX检索结果的第一条序列是什么？结果中的
E-value值和bit得分是多少？它看上去是正确的吗？该生物属于哪一类？它能解释
翼龙的进化问题吗？
d. 试解释为什么BLASTN和BLASTX的结果差异很大。
2．古基因组学 II——BLAST用于基因内含子的影响
你得到了一条新的翼龙的序列，保存在dino2.fa中，你尝试使用BLASTN和BLASTX在NR数据库及蛋白质数据库中检索，但未得到相关结果（请尝试一次确认一下）。由于该序列来自翼龙基因组的富含基因区，因此你怀疑也许是序列中包含了一段未知基因使得BLAST无法正确匹配，为了验证该可能性，你决定对翼龙序列的模体结合位点进行研究，看是否能够提取出dino2.fa序列中的基因外显子。在进行了多次检索后，你找到了56个翼龙基因组DNA片段的BLASTX结果及对应的已知蛋白质序列，同时推断这些结果中很可能包含外显子和结合位点，由于BLASTX的结果经常无法准确的对应外显子的边界，（例如，因为在靠近结合区氨基酸的改变造成序列的删节、在结合位点/内含子的转录错误导致序列延长）你建立了一个数据库，大约包含BLASTX检索结果的从前往后数的前25个碱基，命名为“3primesplicesites.txt”（为的是找到3’的结合位点模体），另一数据集包括BLASTX检索结果的从后往前数的25个碱基，命名为“5primesplicesites.txt”（为的是找到5’的结合位点模体）。
a. 使用以上2个数据集运行MEME（https://www.360docs.net/doc/1c19151110.html,/meme/website/meme.html），
找到56个序列中的共同模体作为恐龙3’和5’的结合位点，你找到了哪些模体？使
用如下所示的单字符编码列出相同的模体：
注意：当使用MEME进行“多级相似序列”输出时将会出现问题，第四条核苷序列不会在模体的指定位置输出。因此，检查“Simplified pos.-specific probability matrix”或者“Information content”图表以确定模体的任意位置都包含所有的4个核苷。
b. 编写python程序读入两条相似序列，找到在“dino2.fa”序列中推断的内含子，并输
出预测结合的mRNA。程序应从命令行接收2个参数，即文件名，其一为包含2条
相似序列的文件，每条序列占一行（分别是3’和5’结合位点的相似序列），每条序
列都由上面提到的字符组成，即A，T，C和G。另一文件包含FASTA格式的DNA序
列，你不能假设这条序列只占一行（像上次作业中多次假设的那样），登陆
http://ngfnblast.gbf.de/docs/fasta.html查阅FASTA格式的相关信息，你的程序应能搜
索3’和5’结合位点的相似序列并预测结合的mRNA，我们假设与现在很多的有机体
不同，恐龙中内含子的第一个核苷是5’结合位点模体的5’ 核苷，内含子的最后一
个核苷是3’结合位点模体的3’ 核苷，将程序命名为splicing.py并在线提交，提示：
从翼龙基因组序列的BLAST检索结果可以看出，翼龙内含子的长度一般有300-600
个碱基。本应在程序中加入查错代码，但这次就没必要了。你可以假设程序运行时
都能接收到2个文件名并包含上面要求的格式，登陆课程网站查看程序运行示例
（splicing.run）。
c. 使用你预测的拼接外显子运行BLASTN程序，你发现哪些似是而非的结果了？
d. 进行BLASTX搜索，这次又发现哪些似是而非的结果？假如是的话，解释为什么
在先前的研究中没有出现这个问题（在移除内含子之前）。
3．模体识别
从课程网站下载pombe.fa文件，它包含从分裂酵母Schizosaccharomyces pombe中得来的75个内含子的3’端fasta格式序列，每条序列-35到-6的碱基段与3’端的结合点相关，即该生命体中通常被认为包含分支信息的内含子区域。（a部分问题不参与考评）
a. 检查文件中的序列，通过用眼睛观察能否发现一些共同的模体？
b. 运行https://www.360docs.net/doc/1c19151110.html,/gibbs/gibbs.html上的Gibbs Motif Sampler分析序
列，模体长度设置为7并将指定不同模体数量的位置留空，采用默认设置运行程序
（详情参见https://www.360docs.net/doc/1c19151110.html,/gibbs/bernoulli.html的操作手册），它发现
了哪些共同的模体？在服务器上重复运行4次，每次都得到相同的结果吗？每次的
碱基出现频率都一样吗？解释这一现象。
c. 在https://www.360docs.net/doc/1c19151110.html,/meme/website/intro.html的MEME上分析这些序列，它发现的
共同模体与Gibbs Motif Sampler运行结果有什么不同？解释结果为什么不同或者是
相同。注意MEME给出的模体信息是已发现的模体，这些信息是什么？
d. 编写python程序，输入参数为一个文件名，该文件包含任意长度和任意数量比对好
的特定DNA模体，计算权重矩阵（特定位置的碱基频率）即先农熵以及每个位置
的信息量。示例输入文件在网站上，文件名为motifs.txt，假设每个核苷的频率分布
相同（每个为25%），计算整条模体的总信息量。
e. 扩展你的程序使它能够利用计算出来的权重矩阵并计算以后输入的文件中的模体
序列的bit score（公式在讲座中提到），输入这些序列及其得分，计算平均得分，计
算出来的平均得分跟前面的信息量相比有什么差别？
f. 基于上面计算的信息量，在一个长度为500000个碱基对的随机DNA序列中，你期
望模体出现次数是多少？随机示例模体保存在random.fa文件中，使用门限值9.8
搜索模体，你发现此模体出现了多少次？这两个数据有何差异？讨论你预测的正确
性。（我们建议你使用编写的python来进行，但是这里就不再检验你的程序了）注意：在这里我们假设每行包含一条DNA序列，其间没有空格，每条序列长度相同（因为它们都是由Gibbs sampler计算的相同模体的比对结果），行数没有限制，你可以假设输入文件中确实包含了一个DNA序列列表（在用户输入方面无错误），下面是程序的运行示例，将你的程序命名为motif.py并在线提交。
4．概率与统计
支原体Mesoplasma floracea的基因组比较特殊，它的G+C（25%）含量非常低并且环状基因组非常小（734kb）。有一种称为“顺时针”解读和一种称为“逆时针”解读的方法。下面问题都是关于顺时针解读的。
a. 假设你从基因组序列中的任意位置开始测序，在顺时针测序中，你期望在遇到G+C
之前出现多少碱基对？同样的问题，逆时针测序呢？
b. 在基因组中连续出现G+C的平均间隔是多少（即插入A+T碱基对的平均数量是多
少）？解释为什么该数字与你在a部分中得到的2个数据相同或者不同。
c. 在基因组注释过程中，基因搜索程序将会用来搜索基因组中686个已预测的编码区
域。在这686个假定的区域中，有432个由物种碱基对的编码区域的相似性得以证
实，其余的区域还未被证实。根据编码某一氨基酸的过程，你估计在编码区中G+C
的含量高于非编码区，在432个已证实的编码区中，337个编码区G+C的含量也高
于25%。在203个证实的非编码区中（该区域实验证明未进行表达转录），仅有17
个G+C的含量高于25%，假如一个特定的预测编码区域的G+C的含量高于25%，
那它确实是一个编码区的概率是多少？在联合概率表中写出你的结果，假设未证实
的编码区/非编码区内容相同，分别验证编码区/非编码区。
d. 在注释Mesoplasma floracea基因组的下一步工作中，你希望找到基因组中所有转录
终止子，根据基因组的大小及编码区的数量，粗略计算得到基因组中大概3%为转
录终止子。基于前面的检测，TerminatorScan程序常用于预测DNA序列中的转录终
止子（假设它出现的概率为0.79）。如果该终止子不出现，程序回有0.15的概率预
测错误。使用贝叶斯规则判断，假设程序用于判断基因组中任意片段的转录终止子，
那正确的预测概率是多少？给出计算过程，该程序在注释Mesoplasma基因组时有
多大用处呢？
e. 另一研究小组发表论文称已估计了Mesoplasma floracea基因组核苷之间的突变概
率，他们的分析证实了转换比颠换发生得更频繁，即A+T比G+C更常见和机体缺
乏基本的DNA修复机制。基于此转换概率，计算任意给定的基因组中的核苷是A、
T、G、C的稳定概率，给DNA突变设定一个一次马尔柯夫模型。Mesoplasma floracea
中碱基的稳定含量是多少？（注意：Pxy表示核苷从X突变到Y的概率）
5． Viterbi算法
根据Karyn T. O'Neil和William F. DeGrado的论文“A Thermodynamic Scale for the Helix-Forming Tendencies of the Commonly Occurring Amino Acids, Science, New Series, V ol. 250, No. 4981. (Nov. 2, 1990), pp. 646-651”中提到的α螺旋倾向性的实验数据，你的实验室提出了一个一次隐马尔柯夫模型来预测α螺旋。该模型包含2个隐藏状态——螺旋（H）和非螺旋（N）。
a. 你认为哪个参数来自实验的α螺旋倾向性的值？
b. 在每个状态下每个氨基酸的HMM发散概率列于表1，两个状态之间的转移概率列于
表2，假设第一个位置的预设概率为P H=0.3和P N=0.7，则序列FSRSNHLSPC和片段
HHHHHHHNNN的联合概率是多少？给出计算过程和结果。
表1 螺旋状态和非螺旋状态的发散概率
表2 HMM的转移概率
c. 你正研究一个你估计含有螺旋片段的蛋白质，你希望使用实验室提出的HMM模型
来找到所有螺旋的可能区域，此序列为CGQALFAASP，使用网格图表对序列进行优化分析，将10个圆排成2行，一个接着一个，将第一行标记为（H），最后一行标记为（N），将每一列标记为1-10，将子序列1-k的优化分析概率填入表中，在存放螺旋的行中遇到螺旋状态的圆k时停止，而子序列1-k的优化概率停止于非螺旋行中第k个螺旋状态，将预设概率和发散概率填入第一列的2个圆中，而第二列的螺旋状态的值需要考虑2个概率：1）在H->H转换后，优化位置1截止于H处的概率；2）在N->H的转换后，优化位置1截止于N处的概率，将这两者之间的最大值填入第二列的H圆中，继续下去直至全部的网格被填满，该序列截止于H处的最优联合概率是多少？该序列截止于N处的最优联合概率是多少？总的最优概率是多少？（提示：反向回溯整个过程）
d. 在做完上面的练习后，你发现实验室里已经有人运用python程序编制的Viterbi算
法，你使用该程序计算1000个残基的序列花了0.1秒，基于该信息你现在应该能够
估计该算法的运行时间了，将此算法用在HMM模型中计算SWISSPROT数据库的
所有序列（53,205,430个残基）需要多少时间？假如你再加入2个状态——β折叠
和转角需要计算多少时间（假设这两个状态与前两个状态相同）？
编程相关问题提示：