生物信息学: 导论与方法 | 总结笔记
课程链接：生物信息学: 导论与方法
1 导论与历史
1.1 1-1 什么是生物信息学
1.2 1-2 生物信息学历史
1.3 1-3 中国大陆的生物信息学
2 序列比对
2.1 2-1 序列比对中的基本概念
2.2 2-2 利用动态规划进行全局比对
2.3 2-3 从全局比对到局部比对
2.4 2-4 考虑仿射空位罚分的序列比对，以及如何计算Needleman-Wunsch算法的时间复杂度
2.5 2-S1 学生课堂报告
2.6 2-S2 Interview with M. S. Waterman Waterman访谈
2.7 2-S3 关于同源、相似性、相似性矩阵和点阵图的补充材料
3 序列数据库搜索
3.1 3-1 序列数据库
3.2 3-2 BLAST算法初探
3.3 3-S1 学生课堂报告
4 马尔可夫模型
4.1 4-1 从状态到马尔可夫链
4.2 4-2 隐马尔可夫模型
4.3 4-3 用隐马尔可夫模型建立预测模型
4.4 4-S1 学生课堂报告
5 新一代测序NGS：重测序的回帖和变异鉴定
5.1 5-1 新一代测序
5.2 5-2 序列回帖和变异鉴定
5.3 5-3 序列回帖和变异鉴定的分析演示
5.4 5-S1 关于回帖、变异鉴定的补充材料
5.5 5-S2 关于基因型鉴定的补充材料
5.6 5-S3 Ion Torrent PGM测序介绍
5.7 5-S4 3730 Sanger测序介绍
5.8 5-S5 学生课堂报告
6 变异的功能预测
6.1 6-1 问题概述
6.2 6-2 记录变异的数据库
6.3 6-3 基于保守性和规则的预测方法：SIFT和PolyPhen
6.4 6-4 基于机器学习分类器的预测方法：SAPRED
6.5 6-S1 支持向量机简介
6.6 6-S2 学生课堂报告
7 新一代测序NGS：转录组分析RNA-Seq
7.1 7-1 转录组介绍
欢迎回到北京大学生物信息学：导论与方法网上课程。 我是北京大学生物信息中心高歌
下面，让我们开始课程，利用深度测序技术来研究转录组
类似于前两周，在接下来的两节课中，我们将首先介绍如何处理RNA-seq等转录组测序技术产生的RNA数据 而后以非编码RNA为例，演示如何利用这些分析方法得到的结果，去进一步探索生物学问题
首先，让我们简要介绍一下相关的背景 所谓转录组(transcriptome)，是指特定类型细胞中全体转录本(transcript)的集合。 换句话说，是细胞特定时刻基因表达谱的一个快照(snapshot of expression profile) 0 在转录组中，既包括人们早已熟悉的、编码蛋白的信使RNA(mRNA)， 也包括人们近年来新发现的、不编码蛋白的miRNA、long non-coding RNA(lncRNA)等非编码RNA。 这些RNA转录本彼此协同作用，共同来调控细胞生长、发育、凋亡等一系列重要的生理过程
因此，对于转录组的研究通常包括定性和定量两个方面。 前者主要是要鉴定(identify)出所有表达的转录本， 后者则要确定这些转录本各自的表达量。
通过对经典PCR扩增反应中每一个循环产物荧光信号的实时检测， 我们可以实现对起始模板的定量分析， 这也就是所谓的“实时荧光定量PCR”(Real Time quantitative Reverse Transcription PCR, Real-Time qRT-PCR) 。 通过正确设计引物(primer)和探针(probe)，qRT-PCR技术可以在很大范围内定量的检测目标转录本的拷贝数，也即表达水平， 0 因此常被作为转录组分析中的金标准(Gold Standard)。 然而，qRT-PCR技术一次只能测定一个转录本的表达水平； 同时，这个技术需要事先知道待检测转录本的序列，难以用来发现未知的转录本。（比如，我们茶树基因组没有出来，基因也大多未知，因此不能大范围利用此技术）
在深度测序获得大规模应用前， microarray是主要的高通量转录组表达分析技术。 微阵列(microarray),也称基因芯片(gene chip)，是通过将几十万个不等的探针(probe)分子固定在约1cm见方的固体片基上制成的。 利用核苷酸分子在形成双链时碱基互补原理，microarray可以一次性检测出样本中所有与探针互补的核苷酸片段， 从而快速得到样本中基因的表达谱(expression profile)。 因此，microarray从上世纪90年代初次问世后，即在生物、医学、农学等领域快速获得了广泛应用。 与qRT-PCR技术相比，Microarray虽然在通量上有了显著的提高，但仍然需要事先确定待检测转录本的序列。 0
表达序列标签(EST)技术通过对一个随机选择的cDNA克隆进行单次测序来获得cDNA的部分序列。 与microarray不同，EST是基于测序的，并不需要事先知道待检测转录本的序列， 因此可以被用来发现新的转录本。 事实上，早在1991年，当时还在NIH的Craig Venter等人就开始利用EST技术来寻找人类的新基因。 然而，由于当时测序技术通量的限制，一次EST通常只能得到几千个转录本的序列，远远无法进行全转录组水平的profiling
深度测序技术的出现，使得研究人员首次可以在全转录组水平利用测序技术同时进行定量与定性分析， 也就是所谓的RNA-Seq技术。 具体来说，首先对生物样品中的RNA反转录为cDNA， 而后，将这些cDNA打碎为较小片段后，上机进行测序。 0 一方面，RNA-Seq技术使得研究人员可以快速确定转录组，进而鉴定存在的可变剪切体(Alternative splicing isoform)， 而这是传统的microarray等技术很难做到的。 另一方面，通过对基因组特定位点上reads深度的计数， 可以对表达量水平进行估计。 所以，RNA-Seq技术使得研究人员可以同时对转录组进行定性与定量的研究
需要注意的是，RNA-Seq本质上是对转录本序列的随机抽样(Random sampling)， 因此，其检测效力(power)和灵敏度(sensitivity)高度依赖于测序的深度。 如果深度不够，就会难以检出低拷贝的基因。 原则上，只有在饱和曲线(saturation curve)达到平台期(plateau)后，才能认为深度足够。
对于哺乳动物转录组，一个经验规则是通常要做到100~150x的coverage 0 在随机抽样(random sampling)的情况下，map到特定转录本上的read数目正比于其表达量(transcript abundance)， 因此我们可以利用落在某个转录本上reads的总数来估计其表达量。 但是另外一个方面，落在一个转录本上reads的数目，也与其长度和总测序深度成正比 例如，这里有A和B两个基因， 假设它们表达量相同，都转录2个转录本， 由于A的长度是B的两倍，那么map到A的reads数目就是map到B的reads数目的两倍。 如果我们只看到这些reads的数目，我们会认为A的表达量是B的两倍， 但这显然是不对的 再如，这里有两次RNA-seq的实验， 我们假定基因B的表达量在两次实验中都是一样的。 0 但是由于第一次实验的测序深度是第二次的两倍，那么观察到B基因的reads也是在第二次的两倍。 那么如果我们只数reads数目，也会认为基因B在第一次实验中的表达量是第二次的两倍， 而这也显然是不正确的
所以我们在实际分析中，通常会将原始的reads数目(raw reads count)利用线性放缩(scaling)，转换为RPKM值来进行归一化(normalization)处理。 RPKM就是一个常用的归一化的方法 这个公式里的C是贴到这段转录本上reads的总数目， N是这次实验总的reads数目——也就是测序深度， L是这段序列的长度。 在假定不同样本中RNA的总体分布一致的前提下，RPKM可以正确处理由于转录本长度和测序深度引起的artifact， 从而使得来自不同基因、不同sequencing run乃至不同样本之间的表达数据彼此之间可以比较，
但需要注意的是，RPKM并非是唯一的归一化方法。 通过考虑不同的误差因素(bias effectors)、引入不同的生物学假设，可以构造不同的归一化方法。 事实上，已有研究表明，相比于后续提出的TMM、DESeq等方法，RPKM方法在样本间差异基因表达检验等分析中的效果并不是最理想。
另一个需要在RNA-Seq技术中引起注意的地方是链特异性(strand-specific)。 我们知道，DNA的两条链都可以转录，形成不同的转录本。 然而，常用的illumina RNA-Seq kit是不分链的， 也就是说，我们无法知道配对的reads哪个方向和转录本一致，哪个和转录本反向互补。 而对于分链的数据，又有两种不同的情况， 在利用dUTP技术进行标记(labeling)的方法——也就是Illumina的strand-specific kit使用的方法中， 第二个read和转录本方向一致，第一个read和转录本反向互补 0 而在另一种SOLiD等平台常用的secondstrand分链方法中，就刚好反过来了。 因此，在分析之前一定要弄清楚自己的数据有没有分链，是怎样分链的。
更进一步的细节，可以参看本节课程中由北京大学生物信息中心侯玫与田丰两位同学共同准备的computer lab视频教程。 好，以上我们简要介绍了转录组研究的基本背景与常见的实验测量技术。 这是本节的思考题， 欢迎有兴趣的同学认真思考，并积极通过论坛与其他同学及助教进行交流讨论 在下节里，我们将从针对RNA-Seq的reads mapping开始，具体介绍相关方法， 我们下节见！
7.2 7-2 RNA测序数据回贴与组装
欢迎回到北京大学生物信息学：导论与方法网上课程。 我是北京大学生物信息中心高歌
下面，让我们继续课程，利用深度测序技术来研究转录组
在上一节中，我们简要介绍了转录组研究的基本背景与常见的实验测量技术 在本节里，我们将从针对RNA-Seq的reads mapping开始，具体介绍相关方法
在第五周中，我们曾经提到Reads mapping往往被作为深度测序数据分析的第一步， 其质量的好坏以及速度的快慢，都会直接对后续的分析工作产生影响 同样是基于深度测序技术，因此在reads长度、数量、质量等方面，RNA-Seq产生的reads与之前基因组重测序产生的DNA reads具有相似的性质： 0 它们长度短 ， 数量大， 质量参差不齐，错误率较高
然而，由于RNA-Seq测序数据来源于RNA转录本， 特别的，在DNA转录成mRNA的过程中，内含子被切掉，外显子会在剪切位点连接到一起。 对于这些跨过剪切位点的reads，也就是所谓的junction reads，如果不把它们从中间断开，就无法准确贴到基因组上。
这些junction reads是确定剪切位点的直接证据，对于正确重构转录本结构至关重要。 例如，图中的横跨exon 1和exon 3的junction reads，就可以直接支持存在一个由外显子1和3直接相连、中间不包括exon 2的转录本的存在 类似的，下图中的两个junction reads， 分别支持直接连接exon 1、3，以及exon 3、5的转录本的存在。 0
所以，我们的算法在mapping时，需要“意识到”junction site和intron的存在，从而正确的处理这些junction reads。 具体来说，目前处理这个问题主要有两类策略 其一是join exon。 这个策略的第一步是基于已知转录本中所有的exon来构建所有可能junctions。 需要注意的是，这个library中的junction未必是已知的，而是所有可能的组合。 例如，对于4个exons，就会有6种组合 而后，在mapping时，首先直接采取之前DNA reads类似的unspliced方式基因组比对，将非junction reads map到基因组 对于无法直接map的junction reads，则再与第一步中构建的junction library进行比对。 Join exon策略实质上相当于对之前DNA reads mapping算法的一个”补丁” (patch)， 0 通过构建所有可能的junction库，Join exon策略可以发现新的splicing isoform。然而，对于以前未知的exon，这个策略就无能为力了。
为了克服这个问题，可以使用split reads策略。 与之前类似，split reads策略在mapping时，也是首先采取之前DNA reads类似的unspliced方式将非junction reads map到基因组 而对于无法直接map的junction reads，则参照之前BLAST的方法，切分成若干长度为k的种子(seed)， 再利用这些seed重试mapping 换句话说，也就是在更小的粒度上尝试寻找junction site 最后，再将临近的mapped seeds重新组合起来，来得到最终全read的alignment 与之前的Join exon策略相比，split reads策略因为要在更短的seeds mapping而影响了速度， 但这个策略不依赖于先验(prior)的exon注释，可以发现新的exon乃至新的基因
事实上，目前常用的RNA-Seq工具通常会组合两个策略，以平衡检测的灵敏度与速度。 0
例如，由John Hopkins、Berkeley以及Harvard联合开发的TopHat 2工具中， 就首先采用Join exon策略快速检测已知的junction site， 而后再利用split reads策略来发现新的junctions 值得注意的是，TopHat 2针对不同阶段采用了不同的索引(index)，可以进一步提高mapping的速度
完成mapping只是RNA-Seq数据分析的第一步， 之后我们还需要将这些reads组装成转录本， 并对每个转录本估计相应的表达量
在包括junction reads在内的所有reads均正确map之后， 我们就可以将转录本组装问题描述为一个对有向图(directed graph)的遍历问题。 通过对图中的边给予不同的权值(weight)作为约束(constrain)， 0 就可以利用图论(graph theory)中的寻路算法(path finding algorithm)找到一条或多条符合约束的最优路径及其对应的转录本序列。 以下我们以常用的Cufflinks工具为例，来简要说明一下相关想法
Cufflinks是一个针对RNA-Seq数据进行转录本组装和表达分析的工具软件， 下面我们来看看Cufflinks是怎样工作的。 现在假设，我们不知道这里有这样三个转录本结构， 我们只看到这些reads， 下面怎么做呢？ 首先Cufflinks会去找那些不可能出现在同一个转录本中的片段(fragment) 比如说这里黄色和蓝色的片段，它们就不可能出现在同一个转录本里。 因为如果它们来自同一个转录本，黄色的就应该在蓝色这个相应的位置断开，而不是直接跨过去。 0 同理，红色，黄色，蓝色的片段都是两两互不相容的，而同一颜色的片段则是彼此相容的。 通过将每个相容片段作为节点，并与和它最近且相容的片段连接，就可以得到重叠图(Overlap graph) 基于精简原则(parsimony principle)，Cufflinks在图中选择能覆盖所有reads的路径中互不相连且最少的一组路径，作为最优路径， 并据此来得到最终的三个转录本集合
原则上，在转录本组装正确完成、且每个exon水平的表达量均利用我们在上一单元中提到的归一化(normalization)正确处理后， 转录本的表达量可以直接从各个exon的表达量推出。 例如，我们假定从基因组上的三个exon可以转录出两个转录本t1与t2 同时，每个exon表达水平在归一化后可以确定为e1=20, e2=40, e3=60 那么，根据转录本的结构，我们可以直接得到转录本表达水平与exon表达水平的关系。 例如，exon1只在转录本1中出现，那么它的所有的表达就完全由转录本1来贡献； 0 类似的，exon3在转录本1和2中同时出现，因此这两个转录本都对其表达有贡献。 所以，我们认为，外显子1的表达量就应该对应转录本1的表达量， 而外显子3的表达量应该对应为转录本1和转录本2表达量的和。 以此类推，我们就可以推出转录本1和2各自的表达量：20、40
当然，考虑到对reads的分配本质上是由转录本拼接算法决定的，在实际中这个问题会更加复杂， 例如，在Cufflinks中，对于reads的分配还会同时考虑长度分布等因素。 事实上，为了更为准确的估计表达量，常常会采用EM等方法来反复迭代的考虑转录本组装与表达量估计。
好，以上我们从针对RNA-Seq的reads mapping开始，简要介绍了基于RNA-Seq数据的转录本组装和表达量估计方法。 更进一步的操作细节，请参看本节内容中由北京大学生物信息中心侯玫与田丰两位同学准备的Computer Lab视频教程。 这是本节的思考题，欢迎有兴趣的同学认真思考，并积极通过论坛与其他同学及助教进行交流讨论 0 下节见！
7.3 7-3 RNA-seq 数据分析
大家好，这是computer lab，将由我和田丰给大家介绍RNA-Seq分析当中常用的工具套件TopHat和cufflinks
这是这个视频中主要介绍的工具，前面四个都是在Linux或者Mac OS操作系统上运行的命令行工具，cummeRbund是一个R的包
这是RNA-Seq分析的流程， 从原始数据开始，进行reads回帖，到拼接转录本，计算表达量，分析差异表达，最后可视化分析结果，我们将为大家讲解并演示每一步的过程。
下面先介绍TopHat。TopHat是一个把reads回帖到基因组上的工具。下面来看看它是怎样工作的。 首先，用Bowtie把reads回帖到基因组上， 然后通过拼接，我们就会在基因组上看到一些reads堆叠起来的区域，称为consensus， 0 这些consensus可能是一个真的外显子，也有可能是几个外显子拼在一起的，或者一些别的情况。 我们知道，经典的剪切位点一般都有GT和AG这样的序列标志，接下来，在consensus的边界和内部， TopHat会去找到这样的剪切位点，并且得到它们可能的组合。 然后，对于那些没有被Bowtie贴到基因组上的reads，TopHat会对它们建立索引， 去和这些可能的剪切位点比对，这样就把跨越剪切位点的reads准确地贴到了基因组上。
下面给大家讲讲TopHat中比较重要的一些命令行选项。
首先是关于插入片段长度的选项。在RNA-Seq中，会把mRNA打断成小片段，然后对片段长度进行一定筛选后拿去测序，如果选择的片段长度是300bp， 两端各测序75bp的reads，中间的插入片段长度就应该设为150bp。
然后下面这个是设置插入片段长度的标准差，如果选择的片段长度比较集中，这个值可以设置得小一些，反之应该设置得大一些。
-G选项是提供一个已有的注释文件。如果你分析的基因组被注释得比较好了，最好能够提供这个文件， 0 这时TopHat就会先把reads往转录组上贴，没有贴到转录组上的再往基因组上帖，最后把结果合并起来。 我们知道大多数转录组都是比基因组小得多的，而且junction reads可以直接贴到转录本上，所以这样回帖的效率和准确度都会得到提高。
下面讲讲library type， 标准的Illumina平台是不分链的，也就是说，我们无法知道配对的reads哪个方向和转录本一致，哪个和转录本反向互补。 对于分链的数据，也有两种情况，在firststrand这种分链方法中，第二个read和转录本方向一致，第一个read和转录本反向互补， 而在另一种fr-secondstrand分链方法中，就刚好反过来了。 所以大家在分析的时候一定要弄清楚自己的数据有没有分链，是怎样分链的。
下面是TopHat的演示，本视频所用的数据集是GSE32038，大家可以去GEO的网站上下载练习。 这是一个模拟的RNA-Seq数据集，双端测序，有两种处理，每种处理有3个重复，这里C代表处理，R代表重复，下面用C1 R1进行演示。
在运行TopHat之前，需要有一些准备工作。 0 首先要有参考序列fasta文件，也就是通常说的基因组序列。 这是黑腹果蝇的RNA-Seq数据，所以我们事先下载了果蝇的基因组，也就是这个文件。我们来看一下它。 由于TopHat是用Bowtie2回帖reads，我们首先需要建立Bowtie2的索引文件。 输入这个命令就开始了。 这个过程需要一些时间，我们直接来看运行结果， 这些bt2结尾的就是Bowtie2的索引文件。 我们还需要reads的fastq文件，双端测序的数据，两个fastq文件分别以下划线1和2这样的形式结尾， 给大家看一下这个文件是这样的。 这里每4行代表一个read， 第一行是它的名字，第二行是它的序列，第三行是一个+号，第四行是代表它上面每个碱基质量的一串字符。 0 由于这是模拟的数据，所以都是I。 在实际分析过程中，需要对拿到的数据进行质量评估和过滤等一系列预处理工作，这些工作是非常重要的。
然后需要准备这个注释文件，当然它不是必须的。 它可以是GTF或者GFF3格式的文件，对于注释得比较好的基因组，在UCSC可以下载。 我们来看一下这个文件。 我们来看第一行。比如说这是一个exon，它在这条染色体上， 这是它的起始位置，结束位置，它在正链上，然后它属于这个基因， 这个转录本，是第一个外显子。 这样的注释就可以给TopHat提供剪切位点的信息。
准备好这些之后，就可以运行TopHat了， 0 这里-p是线程数，-G是注释文件，-o是输出文件夹， 选项之后就是参考序列的索引，最后是两个reads的fastq文件。
下面我们输入这个命令来运行TopHat。 由于这个过程时间比较长，我们直接来看之前跑好的结果。 这里C1_R1就是TopHat的输出目录，我们来看一下它。 首先我们来看一下logs， 这个tophat.log文件就是TopHat运行的整个过程， 下面我们来看一下这个align_summary文件，这个文件是reads回帖的一些统计信息， 我们可以看到，reads全都贴到了基因组上，有1.6%的reads贴到了多个位置， 由于这是一组模拟的数据，所以会有100%的回帖比例， 0 这在真实数据当中基本上是不可能的，一般来说真实数据90%以上的回帖比例就非常好了，当然百分之六七十也是一个可以接受的范围。
下面我们来看一下这个bam文件。 这个文件详细地记录了reads回帖到基因组上的情况， 由于这是一个二进制文件，我们需要用samtools查看它。 我们来看一下这个bam文件。 这个bam文件的每一行是对read回帖到基因组上的情况的一个描述，比如说我们来看这一行， 这个16是这个read的名字， 这个数字是对回帖情况的一个描述，比如有没有贴上，是贴到了正链还是负链，和它配对的read有没有贴上等等。 这是贴到的参考序列的名字， 这是这个read最左端在参考序列上的位置， 0 这是回帖的质量， 这个75M表示这个read的75个碱基都贴到了基因组上。 我们可以看到还有这样的情况， 这就是一个read的前面66个碱基贴到了基因组上，然后在基因组上空了76个碱基之后又贴上了这个read的最后9个碱基， 这就是read贴到了剪切位点上。 然后这是与这个read配对的另外一个read在参考序列上的位置， 这是片段的长度。 关于TopHat的介绍就到这里，下面是关于Cufflinks的介绍。
Cufflinks是一套列拼接转录本、计算表达量、计算差异表达的的工具。 下面先介绍拼接转录本、计算表达量这部分。 0 前面已经说过，reads来自于转录本， 如果一个基因有多个转录本结构，就像图里面这样，有红、蓝、黄三种结构的转录本， 除了红蓝黄这三种颜色的reads之外，其他reads我们无法知道它到底是来自于哪个转录本的。 在这样的情况下，我们怎样尽可能地还原事情的真相呢？ Cufflinks就是尽可能地拼接出最有可能的转录本结构，并且估计它的表达量。
下面介绍一下Cufflinks中比较重要的一些命令行选项。
大-G是提供一个注释文件，并且告诉Cufflinks不要去拼接新的转录本，只能用注释文件里提供的转录本。
小-g也是提供一个注释文件，但是Cufflinks会在这些已知转录本的指导下，拼接新的转录本。
-u是告诉Cufflinks用更准确的方法去处理贴到多个位点上的reads。 如果没有-u，Cufflinks只会把这些reads简单地平均分配。比如一个read贴到了10个位置，那么每个位置分得十分之一， 0 加上-u之后，Cufflinks会更准确地分配reads， 比如会先进行平均分配，然后按照这10个位置各自的表达量，计算read被分配到每个位置的概率。 实际上Cufflinks会用EM算法进行迭代，计算在观察到当前数据的情况下，最有可能的reads分配。
library type和TopHat里面差不多。
下面是一个Cufflinks的演示。这几个选项前面都已经介绍过了，最后这个bam文件就是刚刚TopHat的运行结果。 下面我们来运行这个命令。 这需要一定时间，下面我们直接来看运行的结果。 这个文件夹就是刚刚Cufflinks的输出目录，我们来看一下它。 这个文件就是Cufflinks拼接的转录本，我们来看一下。
上一小节我们对Tophat和Cufflinks作了介绍。 现在我们将演示Cuffmerge, Cuffdiff 和CummeRbund的使用。 00
当我们使用Cufflinks处理多个数据之后，我们需要将转录本数据整合为一个全面的转录本集合， 01 Cuffmerge是一个将Cufflinks生成的gtf文件融合成一个更加全面的转录本注释结果的工具。 02
如图所示，图中的6个转录本被整合为一个转录本集合。 03 同时，我们还可以利用基因组注释文件，获得更加准确可靠的结果。 04 合并后的转录本集合为计算每个基因和转录本的表达量提供了一个统一的基础。 05
下面介绍一下Cuffmerge中重要的命令行选项： 06 -g参数指向参考注释GTF文件 07 -p参数决定线程数，-s 参数指向基因组DNA序列。需要注意的是： 08 如果是一个文件夹，则每个contig是一个fasta文件； 09 如果是一个fasta文件，则所有的contigs都需要在里面。
下面是Cuffmerge的一个简单演示 这些参数都已经讲过 最后一项是一个列表，内容包含经过Cufflinks拼接的转录本的文件路径 下面我们来执行这个命令 首先我们需要使用cut命令创建一个所有拼接的转录本的文件路径列表 然后运行cuffmerge 运行后的结果储存在merge_asm这个文件夹里面 其文件夹内包含一个logs文件夹以及一个.gtf文件，也就是我们经过整合的转录本文件
当我们利用Cufflinks获得了拼接的转录本时，我们可以计算不同样品中转录本的表达量 计算的简单原理在于测序深度和外显子长度一定时，Read的数量与对应的转录本数量成正比， 0 通过对reads进行计数计算转录本的表达量。同时Cuffdiff可以计算不同条件下转录本表达水平的显著性差异。
下面介绍一下cuffdiff中重要的命令行选项 -u 命令指Cuffdiff对回贴到基因组中多个位置的read进行一个初步估计，然后加权分配到各个基因组位置 而不是简单地平均分配，其功能与Cufflinks中的u命令相同，这里不再重复讲述。 -L 为每个样品标上名称
下面是Cuffdiff的一个简单演示，以上参数都已经讲过， 接下来是Cuffmerge产生的gtf文件，Cuffdiff需要它提供的注释进行初始转录产物和可变剪切等定量分析。 最后是Tophat产生的bam文件，如果一个样品有多个实验重复， 那么我们需要提供bam文件列表，文件名之间以逗号隔开。 下面我们来执行这个命令 0 这是运行之后的结果 cuffdiff输出的文件在diff_out目录之下 其中包括一些按类别统计的表达水平结果，如具有相同的转录起始位点，或具有相同编码区的转录本的表达水平，我们可以利用它们进行下一步分析。
当我们对数据进行分析之后，我们希望能够将结果可视化，CummeRbund作为一个分析结果数据的软件包，在作图方面具有出色的能力。 我们也可以在R环境下输入以下命令进行下载和安装。 常见的绘图方式有：密度分布图；散点图；火山图以及柱形图 这些命令的输入都是cuffdiff的输出文件，经过处理后的结果，我们会在之后进行详细的介绍 这是为不同条件的样本标记名称 现在我们对CummeRbund进行简单演示 首先我们输入R命令，进入R环境 0 载入预先安装好的软件包CummeRbund，然后使用readCufflink命令载入cuffdiff产生的数据 现在我们开始作图 这张密度分布图表示不同表达水平的转录本的密度分布 这张散点图表示两种条件下转录本表达水平的情况 位于直线两侧的点对应的转录本其表达水平具有不同条件偏好 这张火山图代表不同条件下基因表达是否具有显著性差异，横坐标指不同条件下基因表达水平之比的对数值，纵坐标为T检验中p值的对数值 因此高于一定阈值的点对应的基因存在显著性差异表达 最后我们可以使用getGene命令对特定的基因进行代入分析 这是两种条件下，该特定基因的表达水平的差异 这是两种条件下该特定基因不同的转录本的表达差异 0
更多内容请参看网站和使用手册 最后谢谢大家观看
7.4 7-S1 Maynard Olson教授访谈
我是Maynard Olson，应魏老师邀请来做一个视频。
一些给上北京大学生物信息学MOOC的同学的话，会涉及到人类基因组计划 你们知道，基因组学和信息技术是一起成长起来的， 但基因组学是借着个人电脑、强劲的桌面工作站、数据库技术、网络、互联网这些技术的发展而发展起来的。 我估计对于很多MOOC的同学而言是很难认识到30年前的生物学家压根儿不用电脑，除了极个别的特例，比如蛋白结晶学家。 0年，我还是个副教授，在圣路易斯华盛顿大学做遗传学。 那时我通过和一个结晶学家沟通，从而有机会使用电脑。那时候的电脑还是那种能够放满一整个屋子， 但是计算性能还不如现在一个iPhone来得快的单处理器系统。 0 我有段时间就根本不在实验室，一直在折腾楼里的线路和电脑终端的连接。 我弄好了之后，我部门里的每个人经过时都会停下来看看它，然后问我我打算拿它干什么。 那时候遗传学部门还没有别的电脑。 甚至办公室里都没有文字处理机。 几个月后，我买了第二个终端，我的小伙伴们更加吃惊了。
我打算把一千五百万碱基长的酵母基因组里的[限制性]内切[酶]位点图谱给详细地描绘出来。 那时候测序还是个太遥远的问题；先描绘出详细的物理图谱才是重中之重。 二十年后，我描述了我们的方法， 基本上就是人类基因组[计划]里被自动化的那个方法。 2001年，Nature出了一期专刊，描述人类参考基因组序列的粗略图。 0 你可以在那期Nature里找到这个图谱的制作技术的说明。 我们最近有收到一封学生的来信，希望我能解释最后那张图谱是怎么从展示出来的概要数据造出来的。 我相信MOOC的各位同学应该都可以自己把这张图谱造出来了。 然后我们要说到在一个巨大的规模上造这种图谱，它的算法复杂度会有多高。 如果是酵母的话……差不多一万多的数据。 如果是人的话……那个数大概是……差不多快到一百万/差不多10的13次那么多。
计算生物学到处都是问题，就比如用这种数据来造图谱。 把问题泛型化的时候，就像做这种项目经常会遇到的，一定要注意考虑到数据中那些不可避免的错误。 计算机科学家总想把这类问题过分抽象化甚至是理想化。 而在湿实验生物学的现实世界里，数据不够好、弄得最后造出来的组合序列图谱在逻辑上各种说不通这种情况就是家常便饭，而不是什么偶然情况。 0 我们造物理图谱的目的是帮助人们在较长的尺度上拼接基因组序列。 尽管到2001年时，人们都知道全基因组拼接对于捕获基因组的基本长程结构是够用的了，但是为了把人类基因组计划做完， 还是得用到基于分子克隆造的图谱。 由于我们这里有很多人来帮忙做完人类基因组计划？所以我们的主要问题，也就是由于重复序列引起的一些局部序列问题， 最后我们有几百个专家把人类基因组计划做完——从2001年就开始做，做到2003年才把序列上的细枝末节的问题给解决。 当然，人类参考基因组序列还在继续修正，主要是针对那些高度重复区域，比如那些围绕着着丝粒的 端粒即便是用目前的方法也没法测到。 不过，为了把人类基因组计划完成而花费的巨大投资也在后续使用人类参考基因组序列的研究中提供了大量的生物学证据。 目前第37版人类基因组序列每天都在被世界各地做基因组学研究的地方所使用。 对于研究人类遗传变异，使用下一代测序技术测序人类，还有目前生物信息学中德研究重点而言，这是目前现存的质量最好的几个G的bp的序列。 0
接下来是给各位MOOC同学一些鼓励的话，我就说些我做基因组学这四十年来的一些高级的教训吧。就讲一些跟今天有关的话吧。
首先呢，要深入，无论是生物学，还是计算机科学。 生物学是丰富多彩而又复杂多变的。 去和纯生物学家打交道去吧，没有别的捷径可走。多和他们谈话，读他们的文章，去开他们的会，学学他们的思维方式。
另外一方面，计算机科学家和统计学家提供了生物学所缺乏的坚实可靠的理论基础。 我一开始给我们的物理图谱写模拟软件的时候，我只知道用交换排序来排序。 它的计算复杂度是n方。 我没做出什么来，直到我学了些计算机科学的理论，而不仅仅是实用的编程技术。 不用担心将来会有太多数据。 基因组学家只是IT革命里的一个小小角色。 0 IT产业目前在硬件，网络，还有软件的创新上都没什么进展，估计要等到中国的每个人都可以同时用他/她的移动设备现场录高清视频才算真有进展了。 于是我们就可以跟过去30年一样，继续搭着这股技术创新的旋风。 解决你的问题时，尽量使用比较一般常用的解，而不是陷到过分为某生物学应用特别定制的技术里。 重要的是，坚持你做的东西，和基因组序列紧紧绑在一起。
这30年来我已经涉足了一个又一个不同的生物学领域 蛋白组是个非常棒的例子，描述了DNA变为蛋白质的生化过程。 生物学的信息是数码的， 真正的信息被以四进制的形式编码成长长的数字，写到了双螺旋里面。 这是60年前发现的、但到现在为止都是最伟大的、最出乎意料的科学发现之一。 0 一千年以后的生物学家应该也会像今天的我们这样盯着一样的序列看吧 但是不像天文学家，他们应该没办法造出更高分辨率的望远镜。 DNA序列中德最高分辨率应该就是碱基自己了。 序列就是它自己。 我们基本已经知道它的大部分了。 我们的下一代估计能够用难以想象的方法读出它背后所隐藏的真实。
我这一代的基因组学家基本上就是纠结于获得参考基因组序列了。 而由你们大部分MOOC学生所代表的这一代人正在面对着非常困难的问题。开始尝试弄明白我在说什么吧。 我们只是刚刚涉足到了水面而已。 祝你好运，学业有成。 0
7.5 7-S2 学生课堂报告
今天我们主要讲一下这个RNA-Seq的情况。
首先我们要说一下RNA-Seq可以做些什么事情。
RNA-Seq的话，从一个比较简单的来说，它可以研究一下每个基因表达的多少。 因为比较简单，就说一下RPKM就行了。
第二等级，就是我们要研究一下，我们经过测序之后，发现了有哪些转录本，就是不是基因这个等级，而是这个基因下面，包括转录本这个等级， 有哪些转录本。然后这些转录本各各自表达了多少。
然后，在我们知道这两点的情况下，我们要知道它们内部比较有没有这种可变剪切。 然后它们不同样本比较之间，有木有这种表达差异。
其实RNA-Seq还可以研究更深入的些东西，比如说RNA-editing，还有eQTL。 0 因此，我们今天主要讲的内容是在level2这一部分。
当我们拿出以上这些结果的时候，我们可能会拿到一个gene list，然后我们可以根据这个gene list做GO或者pathway的分析。 恩，然后这就RNA-Seq的一个宏图。 我们今天主要讲的还是这个第二部分。 如果需要回答有哪些转录本，还有它们这些转录本表达量如何，我们就需要一种与之前基因组研究不同的一种回帖方法。
然后我们就介绍回帖方法。
最开始的一种思路就是说我们，就是使用比对基因组的这些回帖软件。它们不能把这些读段给分开。 软件的典型，就比如说是bwa或者bowtie这些软件。 它们一个特点就是，它们不能把读段给分开。比所以说它们就把基因组片段给分开，把它们命名成新的文件， 0
比如说这个样子，这是截取的一小部分，它们都是属于这个基因，但是这个基因有3个不同的转录本， 以说它们实际上是长这个样子，他们的位置基因是重合的，只是在最左边的5’ UTR区域，有点差别，然后所以是他们的 序列也有很多地方是相似的。 这就带来一个非常严重的问题，就是比如说，我在测得了一个读段想往上回帖的时候， 如果是回帖到这个部位的话，我们会发现，我们知道它属于这个基因，当我们不知道它属于哪个转录本。 如果是按照bowtie或者BWA的方法的话，它会在这3个中任意选出一个转录本，把它帖上去， 就像是这样，这个读段在这3条转录本中都出现了， 最后我们可能，如果根据BWA的结果来看它的话，可能就是错误的结果。我们不知道它是哪个转录本。 但是如果我们把转录本上述到基因的话，它的基因表达量是对的。 那我们能不能看它转录本表达情况呢。
以说我们用一种新的比对方法，可以把这个部分reads给拆开， 0 通过拆开reads，我们可以找到一些junction，我们可以把转录本边界给定出来。在这种情况下，我们使用的reference就是基因组。
只有两种不同的策略，一种就是基于外显子的比对，还有一种就是基于种子延伸的这一种策略。 这两种策略的大概流程如下，我们重点介绍一下基于外显子的比对策略 因为它和之前说的基因组的回贴方法比较相似，速度也非常快。 而这种基于种子延伸的方法呢，它是一个动态规划的算法，速度就比较慢，占内存也比较多，不太实用 虽然这种基于外显子的方法非常快，但是有个非常严重的问题，就是如果你是直接比的话，会发现有个假基因的问题 可能它自身的基因这个地方有一个junction，但假基因没有，但是有SNP 在我们打分的时候，一个junction的分数肯定要大于一个SNP的分数 在这种情况下，这个reads就更倾向于贴到一个假基因上面，造成我们看到假基因大量表达，真实基因不表达的假象 那么怎么应对这种情况呢？ 0
实际上有一篇研究提到一种方法 就是通过先比对到cDNA上，在比对的时候，比对时先比对到类似前面给大家看的ref-RNA的那种文件上面 然后等于我们可以不用考虑会回贴到假基因上面的这种情况了 回贴到cDNA之后根据cDNA再回贴基因组 这时候它就会得到一个准确的比对结果 用这种方法我们最推荐的一个软件是TopHat 这款软件； 用TopHat做回贴的结果bam文件时这样的，我们会发现它在这个地方有一个659bp的插入， 然后它在基因组上的659bp的插入就是长这个样子，这就是一个junction； 类似的这样的回贴的结果我们就可以重构它们的转录本； 这种中间有Gap的读段，我们可以把它再拿出来； 0 然后在回贴的时候我们就可以考虑它们是不是在回贴在两个Exon区域中间的地方，是一个junction， 然后我们就可以根据这些reads把junction定出来，定出junction之后， 我们就可以把它们整个Exon连接起来，Exon连接起来就是一个大的转录本；
然后，这是一种策略，另一种策略实际上它考虑的就不是根据回贴这种方法找转录本，而是根据一种de novo的组装的办法 当然这种办法主要还是应用在没有转录本的情况下，有转录本的情况下尽量不要做组装。 它的理论基础是一种图的结构，类似质粒一样的环形DNA做测序的话我们得到这么五个短片段的情况下， 我们可以对它们进行拆分。 就是比如说这五个段片段当中我们可以拆出这么多种dimer，比如说AA、AT,在这些短片段当中都可以发现， 他们一共加起来有这么多种。在得到dimer的基础上我们还可以去拆出一个trimer出来就是这种三联的情况。 然后我们如果去找一下Dimer 和 Trimer的关系，就会找出这样一种对应关系， 0 然后我们就可以找出有前后关系的两个Dimer之间的关系 我们根据这种数据结构可以画出一张图出来 就是De Bruijn这种图 然后我们就想办法怎么把这张图完整的推出来 这时候我们可以用到一条比较长的reads，拿它当参考 比如说这个reads是ATGGCGT，我们可以根据这个找出它的起始 然后把它的起始定在这个图的AT这个地方 然后往前推，推到这时候我们发现有个G它有两条路 这时候怎么办呢？这时候我们同样是拿这个reads来定位 定位的话发现ATG后面跟的是一个G，这时候我们就把他走这条路 0 然后根据这个原理我们就可以把整个图推出来 在这张图推出来之后我们发现是一个循环 其实它的基因组序列就是这样子，这就是组装的大致原理
当然，在能不组装的情况下 特别是我们研究的主要是类似人、小鼠这样的常见物种的话 我们还是推荐尽量不要使用组装，使用Cufflinks这种软件比较好 当然，如果大家需要研究可变剪切几种类型的话 就像我们上节课所说的八种可变剪切 建议大家还是看一下右边这篇08年的Science文章，他对这八种类型做了详细介绍 同时如果是需要做这方面分析的话可以看看左边这篇文章 0 它详细介绍了使用的这些方法，里面使用了一种叫MISO的软件并介绍了它的原理
最后我们就要来说一下表达差异的问题。 首先我们要明白我们对表达做了一个什么样的指标，用什么样的定义表达。 我们就使用，比如在Cufflinks中就使用了FPKM这个概念。 它和RPKM的区别就是它是基于双端的一个概念。 然后我们为什么要做FPKM呢？ 原因就是要对它做一个Nomalization，就是比如说如果是3,4这两个转录本的话，它们的reads 分别是这么多,4的reads明显要高于3。 在这种情况下，如果单纯统计reads数目，我们就会发现，4明显要比3高。 它们实际上根据长度做一下normalization的话，会发现3和4之间它们是没有表达差异的，这也和我们图中预期的他们gaps没有差异这个点是符合的 所以说FPKM是一种gaps的体现，也是反映了两种基因他们不同的转录数目。 0 然后我们就说一下根据这个结果，我们可以怎么去定义这个转录本的有哪些转录本的边界。 同样，我们测到一条reads，我们怎么定义它是在哪条转录本上边? 然后就比如说黄色的这一条双端的序列，如果我们认为它是在C上面的话，我们会发现C非常长。 如果我们定义c长度是500的话，根据这附图我们会发现如果它长度时500，则它的概率会非常非常小，它是一个正态分布。 中间如果是150的话，它如果在C上面，长度500，它会是一个非常小的概率事件。 但如果这个最上面的黄色的这条读段是在B上面的话，我们会发现它实际上概率会上升一点。 然后如果这个黄色的是在黄色的A转录本上的话，它在A上投影的长度会是150，这时候它就会非常符合正态分布的最高点。 所以根据他们这一点，我们就可以分别知道，他们在A, B, C这三个转录本上面的概率是多少。 根据这一条reads的结果，我们可以把它反映到图上所有的reads上面。就可以画出一个这种图。 我们就会知道在这个基因上面，这几种转录本他们之间是占有多大的比例，然后就可以得到这么一个结果。 00 这就是一个同一个基因不同转录本的一个表达差异。 01 然后在同一个基因不同样本我们还可以再结合多个样本来说，可以画出一个图a的这种图示。 02 如果是我们不考虑转录本的问题，然后就拿所有基因，所有的样本我们可以画出右边的图b中这种heatmap这种图。 03 这都是我们RNA-Seq研究可以做到的一些事情。 04
然后如果是有兴趣尝试这个流程的话，有一篇文章描述的非常好，就是这一篇Nature Protocol上的。 05 它大概是描写了这样一个流程，就是从Tophat开始，然后到最后分析，是这样一条主线。 06 然后这是今天的参考文献，谢谢大家!