課程鏈接：生物信息學: 導論與方法
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 然后這是今天的參考文獻，謝謝大家!