參考學(xué)習(xí)《R語言與Bioconductor生物信息學(xué)應(yīng)用》第六章

前言

Y叔的公眾號biobabble發(fā)過一篇【聽說你想學(xué)R宇植？】,七月份發(fā)的但昨天推送了一次，所以我看到了埋心，看到了對《R語言與Bioconductor生物信息學(xué)應(yīng)用》這本書的強力吐槽，而我發(fā)的這篇筆記闲坎，連同上次發(fā)的Bioconductor分析基因芯片數(shù)據(jù)都是來自于對該書內(nèi)容的提煉和學(xué)習(xí)茬斧。所以呢我覺得有必要在這里說幾句腰懂。

坦白說项秉，這本書確實有很多問題，我自己講講幾點吧：首先娄蔼，它有點過時了，公眾號評論就有人說基因芯片分析過時了锚沸，我個人覺得不客觀涕癣，但確實這本書有些過時哗蜈，由此產(chǎn)生了一系列問題坠韩，特別要提到的就是代碼的可重復(fù)性，我個人在運行書中一些代碼時很多時候會不work绽昼，然后我會自己思考怎么讓它work，實在不行就放棄；第二點，它不適合R入門的人菱农，特別是你去看它的第一章柿估，云里霧里，好歹我也摸了半年R了啊秫舌，好吧的妖，我就直接跳過了；第三點足陨，它的流程不完整嫂粟，就是不是很連續(xù)的都能讓你從頭到尾的go through下去，自然心里會感覺不爽快墨缘。如此種種星虹。

接著你心里已經(jīng)有準備聽我下面要將的“but”的跳轉(zhuǎn)了，我為什么會學(xué)它并做筆記镊讼，乃至分享出來宽涌？首先，在我需要學(xué)習(xí)芯片和RNAseq分析的時候身邊剛好有這本書蝶棋，我也不知道是誰的，好吧兼贸，那就拿起來看看献酗，發(fā)現(xiàn)正是我需要的，所以我看它罕偎，第一章看不太懂颜及，沒意思，我也不想看俏站，直接跳過從第二章看起，到現(xiàn)在整本書基本看了大半墨林，看過的代碼都嘗試著去運行過，確實有所收獲旭等，所以我會寫前言搔耕，算是對這本書的客觀評價吧；其次弃榨，我想談?wù)動心男┦斋@，我本人可以算是有編程基礎(chǔ)的娜饵，算不得菜鳥腊凶，但是對于基因芯片的基礎(chǔ)也好，RNAseq乃至基因組分析流程褐缠、背景等等可以說是菜鳥的不能再菜鳥风瘦，這本書給了我對芯片數(shù)據(jù)來源、處理流程的一些基本認知万搔，其實這在一些國內(nèi)資源上是找不到的瞬雹。有一點我心里非常的不服氣，為什么我聽說中國做生信非常厲害的人很多酗捌，找得到的中文資料卻很少？為什么百度其他搜索引擎很強大尚镰，一涉及科研領(lǐng)域就非常之垃圾哪廓？這也是我挺佩服生信技能樹或類似的這樣的團隊以及相關(guān)個人，當(dāng)我們在噴一些書籍很垃圾分俯，而實際它確實有很多問題的時候，我們能不能貢獻自己的力量呢澳迫？幾個月的學(xué)習(xí)里我深知自己才能有限橄登，所知甚少讥此，所以不斷模仿和記錄。我把這些筆記陳列出來并不是它寫的有多好萄喳，多值得模仿他巨，而是它能夠給予我們新的知識，又能夠在我們忘記時方便查詢染突。學(xué)習(xí)必然是一個探尋和思索的過程，技能的掌握它不是一本書可以帶給你的也榄，特別是一本技術(shù)類書籍司志，它給了你一個看似可行的方案，你要實際去操作它囚霸，然后心里給予評價激才，在你不確定時，需要多方面整理實踐不同的解決方案吨述，然后找到自己的出路钞脂，建立自己對該某個問題處理的完整體系。

這篇筆記并不會帶你真正學(xué)好RNA-seq的分析冰啃，至少我看完之后沒有刘莹，但它確實可以補一些知識的模塊点弯。它不適應(yīng)入門R矿咕，也不適應(yīng)完全模仿做具體的分析，而是適合你在掌握R之后碳柱，你在做測序分析之前想了解的一些知識莲镣。當(dāng)你知道它非你所需，你可以完全不看它瑞侮。

整體而言，這本書非常短越妈，整體評價偏差慈缔，但國內(nèi)在這個方面學(xué)習(xí)恐怕沒有比這好的中文書籍吧？（所以我建議多看網(wǎng)絡(luò)資料瓤檐，這也是我在交互學(xué)習(xí)的娱节，比如生信媛公眾號文章目錄）我希望那些厲害的人物（教授級人物們）能夠多拓展一些中文科研的視界（提升國內(nèi)人員對生信的整體認知與分析能力，加速學(xué)習(xí)周期）谴古，我也會持續(xù)記錄這樣一些知識稠歉，與生信技能樹里面的小伙伴一起從不同的研究方向，角度去拓展形形色色的基礎(chǔ)知識與理論带饱。我再次強調(diào)，我專注于筆記的目的除了自身學(xué)習(xí)以外勺疼，是當(dāng)你在面對一些概念或者問題的疑惑時执庐，你能鍵入百度搜索后快速地鏈接到本文，并從中找到可執(zhí)行的方案或者幫助你理解轨淌，而不是完完全全整篇的通讀。而當(dāng)你確實是需要對所有的知識點有學(xué)習(xí)的需求時婚被，你再選擇讀它梳虽，不僅僅是這一篇筆記或博文灾茁。不要浪費自己的時間，也要耐心地投入自己的時間禀挫。

如果想了解測序基本原理和知識拓颓，查看我整理的測序與平臺的幾個鏈接。

使用的學(xué)習(xí)數(shù)據(jù)：NCBI SRA (Sequence Read Archive)數(shù)據(jù)庫砰左，數(shù)據(jù)集編號SRA091277

使用的是菊花轉(zhuǎn)錄組樣品场航，分析過程包括原始數(shù)據(jù)獲取、數(shù)據(jù)清理僻造、質(zhì)量控制孩饼、轉(zhuǎn)錄組拼接、轉(zhuǎn)錄本定量立膛、標準化和表達差異分析等過程汽畴。

高通量測序基礎(chǔ)知識

（這里只記錄書中重要的知識點并加以理解）

基于第二代測序建立起來的基因組測序、RNA-seq和Small RNA-seq等應(yīng)用鲁猩，都由樣本收集、文庫制備和測序三個過程組成搅窿，不同之處在于樣品收集和文庫制備隙券。

rna-seq.jpg

第二代測序儀（Illumina）測的序列，無論來自DNA-seq文庫還是RNA-seq文庫渔欢，從左到右依次分為3個區(qū)域：5'接頭（Adapter）區(qū)、目標序列區(qū)和3'接頭區(qū)。當(dāng)多個樣品在一個泳道（Lane）中同時測序時盹憎，我們可以使用多樣品（Multiplex）技術(shù)陪每，具體而言是給每個樣品分配一個不重復(fù)的條形碼（Barcode），其實質(zhì)是一個6-8位的DNA序列檩禾，測序后可以通過這個序列將不同的樣品分開。單端測序（Single end）指僅從正向測序竹握；雙端測序（Paired end）指先從正向測序辆飘，然后從反向測序。Barcode則根據(jù)另一引物“Sb”獨立測序得到芹关。

理論上紧卒，由于制備的RNA-seq文庫插入長度的峰值常常為200或300bp，所以測序應(yīng)該只得到文庫中目標序列5'端開始的部分（這就是常說的read了轴总，以前總搞不懂~）。但是呢功偿，文庫中會有少量目標序列不到測序長度（比如測100bp實際目標序列只有幾十）往堡，那么測序就可能會測到3'端接頭序列，這就是所謂的接頭污染吨瞎。數(shù)據(jù)預(yù)處理時穆咐，如果發(fā)現(xiàn)接頭序列過多，一般是RNA-seq文庫插入長度沒有控制好对湃；如果出現(xiàn)大量全長的3'接頭熟尉，一般是接頭過量洲脂，導(dǎo)致了大量接頭自連（Self ligation）。

實際應(yīng)用中恐锦，估計測序深度使用更多的是達到質(zhì)量標準的有效數(shù)據(jù)量一铅，而不是原始數(shù)據(jù)量。當(dāng)RNA-seq有效數(shù)據(jù)比例過低時潘飘，無法檢測一些低豐度的轉(zhuǎn)錄本卜录，要考慮重新測序。

測序深度艰毒，也叫乘數(shù)，指每個堿基被測序的平均次數(shù)柑土，是用來衡量測序量的首要參數(shù)。測序覆蓋度扮宠，也叫覆蓋率诫欠，指被測序到的堿基占全基因組大小的比率。假如用Illumina 2000測序儀完成一次人類基因組（3G大薪钨恕）單端測序被廓，即可得到300G數(shù)據(jù)（假設(shè)全部是有效數(shù)據(jù)），估計的測序深度即為100倍（300G/3G）昆婿，常見表示為100X蜓斧。將所有讀段比對到人類基因組，如果發(fā)現(xiàn)只有2.7G的堿基至少有1個讀段覆蓋到挎春，其實際測序深度為111X（300G/2.7G），測序覆蓋度為90%（2.7G/3G）能庆。

不同的測序目的要使用不同的測序策略脚线。如DNA組裝使用較多的是2X100bp或更長的雙端測序；RNA-seq使用較多的是100bp或更長的單端鏈特異性測序渠旁；small RNA-seq多用50bp單端測序斯碌。

從測序得到的讀段組裝成目標基因組或者轉(zhuǎn)錄組的基因策略是比對和拼接，比對是把讀段定位到參考基因組或者轉(zhuǎn)錄組上投慈，然后再拼接成連續(xù)序列；拼接也膠從頭組裝（Denovo assembly）加袋，是在沒有參考基因組或者轉(zhuǎn)錄組前提下抱既，根據(jù)讀段之間的重疊區(qū)，把所有讀段拼接起來蚀之，直接獲得基因組或者轉(zhuǎn)錄組（參考文章）捷泞。

轉(zhuǎn)錄組比對常用的軟件有BWA、Bowtie和Tophat失受；拼接常用的軟件是Trinity咏瑟。

基因組和轉(zhuǎn)錄組組裝的不同點：基因組組裝希望盡量獲得唯一或較少的組裝結(jié)果，即一致性序列（Consensus sequence）兄旬。一致性序列上并不是每個位點都只有一種堿基余寥，它實際上只代表該位點出現(xiàn)頻率最高的堿基，存在兩種以上堿基的位點叫做雜合位點。注意呐馆，過分追求一致性會導(dǎo)致過拼接，即來自不同基因的相似序列會被誤拼接到一起续膳∈瞻啵基因組和轉(zhuǎn)錄組組裝可以用一個非常重要的指標N50來評價，即將所有組裝后的序列按照長度從大到小排列社付，累加值接近所有序列長度總和一半時的那個位置對應(yīng)的序列長度。N50越大燕鸽，組裝的結(jié)果越好啼辣，類似的有N90。

測序的質(zhì)量分數(shù)

Phred分數(shù)

測序中常用錯誤概率（Error probability）來表示每個核苷酸測量的準確性党远，還可以賦予一個數(shù)值來更簡便地表示這個意思富弦，叫做測序質(zhì)量分數(shù)（Quality score）。因為這個分數(shù)最開始通過Phred軟件從測序儀生成的色譜圖中得到的花沉，所以也叫做Phred分數(shù)（）媳握。Phred分數(shù)的取值范圍是0到93，可以表示很寬的誤差范圍娩脾，即從1（完全錯誤）到非常低的錯誤率打毛。Phred分數(shù)是最基本的質(zhì)量分數(shù)，其他的質(zhì)量計分標準都來自Phred分數(shù)碰声。

		Base call accuracy
0	1	0
10	0.1	0.9
20	0.01	0.99
30	0.001	0.999
40	0.0001	0.9999
50	0.00001	0.99999

Sanger分數(shù)（Phred+33）

Phred分數(shù)包括2位數(shù)字熬甫，還需要用空格分隔，不方便閱讀瞻颂，又要占用大量存儲空間郑象，實際上文件中不采用它。為了在文件中方便地表示質(zhì)量盖矫，常常將Phred分數(shù)加上33（從33到126變化，ASCII碼正好覆蓋了可打印區(qū)）吐根，并用其ASCII碼值對應(yīng)的字符表示辐马，這就是Sanger分數(shù)。Sanger分數(shù)常用于FASTQ格式的文件冗疮。

Illumina/Solexa分數(shù)（Phred+64）

分數(shù)之間的轉(zhuǎn)換公式：

Solexa分數(shù)的取值范圍是-5到62檩帐，它在FASTQ文件中需要加上64并轉(zhuǎn)換為相應(yīng)的ASCII碼值（59到126）對應(yīng)的字符來表示質(zhì)量湃密。2006年，Illumina公司收購Solexa公司后繼續(xù)沿用其標準泛源。顯著Illmina采用新的標準，采用了Phred分數(shù)（范圍0-62）加64的質(zhì)量分數(shù)没龙。

Sanger分數(shù)（Phred+33）和Illumina分數(shù)（Phred+64）是當(dāng)前應(yīng)用最為普遍的質(zhì)量分數(shù)系統(tǒng)缎玫。

以Phred=20（即常見的Q20標準）為例赃磨，其Sanger分數(shù)為53，對應(yīng)數(shù)字5邻辉；其Illumina分數(shù)為84，對應(yīng)字母T。

Bioconductor中的ShortRead包提供了SolexaQuality和PhredQuality函數(shù)分別生成Illumina分數(shù)和Sanger分數(shù)雷客。

source("http://Bioconductor.org/biocLite.R")
biocLite("ShortRead")
library(ShortRead)
Q=20
PhredQuality(as.integer(Q))
SolexaQuality(as.integer(Q))

> Q=20
> PhredQuality(as.integer(Q))
  A PhredQuality instance of length 1
    width seq
[1]     1 5
> SolexaQuality(as.integer(Q))
  A SolexaQuality instance of length 1
    width seq
[1]     1 T

高通量測序文件格式

FASTQ格式

FASTQ格式是序列文件中常見的一種桥狡，它一般包括四部分：第一部分是由“@”開始皱卓，后面跟著序列的描述信息（對于高通量數(shù)據(jù)部逮，這里是讀段的名稱）兄朋，這點跟FASTA格式是一樣的（起頭的符號不一樣）；第二部分是DNA序列颅和；第三部分是由“+”號開始峡扩，后面或者是讀段的名稱，或者為空教届；第四部分是DNA序列上每個堿基的質(zhì)量分數(shù)案训，每個質(zhì)量分數(shù)對應(yīng)一個DNA堿基。

Bioconductor中的ShortRead包提供了quality函數(shù)可以自動識別FASTQ文件中的質(zhì)量分數(shù)的種類萤衰。

我隨便寫了一個fastq的demo文件，內(nèi)容

@HWUSI-EAS100R:123:COEPYACXX:6:73:941:1973#0/1
GATTTGGGGTTCAAAGCAGTATCGATCAAAATAGTAAAATCCATTTGTTCAACTCACAGTTT
+ HWUSI-EAS100R:123:***********************
!"************************************************************

進行一些操作：

library(ShortRead)
# 讀入FASTQ文件
reads <- readFastq("./demo.fastq")
# 得到質(zhì)量分數(shù)的類型
score_sys <- data.class(quality(reads))
# 得到質(zhì)量分數(shù)
qual <- quality(quality(reads)) # 這里好像一個quality函數(shù)就夠了倦卖，還是尊重原文吧
# 質(zhì)量分數(shù)轉(zhuǎn)為16進制表示
myqual_mat <- charToRaw(as.character(unlist(qual)))

# 如果是Phred+64分數(shù)表示系統(tǒng)
if(score_sys=="SFastqQuality"){
  # 顯示分數(shù)系統(tǒng)類型
  cat("The quality score system is Phred+64", "\n")
  # 輸出原始分數(shù)值
  strtoi(myqual_mat, 16L)-64
}
# 如果是Phred+33分數(shù)表示系統(tǒng)
if(score_sys=="FastqQuality"){
  # 顯示分數(shù)系統(tǒng)類型
  cat("The quality score system is Phred+33", "\n")
  # 輸出原始分數(shù)值
  strtoi(myqual_mat, 16L)-33
}

The quality score system is Phred+33 
 [1] 0 1 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9 9

NCBI中的FASTQ與SRA格式

NCBI的Sequence Read Archive （SRA）數(shù)據(jù)庫怕膛，接受FASTQ格式的高通量數(shù)據(jù)上傳秦踪，并將分數(shù)標準從開始的Illumina分數(shù)轉(zhuǎn)換成了Sanger分數(shù)。而且柠逞，NCBI在讀段名稱前面增加了數(shù)據(jù)在SRA庫的編號和版本景馁，后面增加了讀段的長度。SRA數(shù)據(jù)庫為了節(jié)省存儲空間绰精，將FASTQ文件壓縮為二進制的SRA格式進行保存。用戶如果下載SRA格式數(shù)據(jù)卿樱，可以使用工具軟件fastq-dump將數(shù)據(jù)從SRA格式轉(zhuǎn)回FASTQ格式硫椰。

QUAL格式文件

Solid測序儀產(chǎn)生分離的序列文件（CSFASTA格式）和質(zhì)量文件（QUAL格式）繁调，兩者必須成對出現(xiàn)。QUAL文件采用Phred分數(shù)最爬，而且行必須與FASTA文件中的行一一對應(yīng)涉馁。CSFASTA文件與FASTA格式看似相同，但實際不同爱致，Solid不是用核苷酸殘基表示序列數(shù)據(jù)烤送，而是采用了顏色空間的表示方法。Solid的序列文件如果和質(zhì)量文件合并帮坚，可以產(chǎn)生CSFASTAQ格式的文件，也可以根據(jù)顏色編碼轉(zhuǎn)為真正的FASTAQ格式的文件互艾。

Solid的結(jié)果文件轉(zhuǎn)為標準FASTQ格式的文件需要注意兩個問題：第一是第1個堿基由于來自測序引物试和，必須刪除；第二就是由于Solid顏色空間的編碼是前后依賴的纫普，一旦錯一個阅悍，會導(dǎo)致后面連續(xù)錯誤，一般都將參考基因組反轉(zhuǎn)為顏色空間編碼再進行比對等分析昨稼，而不主張將Solid的結(jié)果文件直接轉(zhuǎn)換為FASTQ文件节视。

關(guān)于常見的生信數(shù)據(jù)文件格式，參見生信常見數(shù)據(jù)格式假栓。

RNA-seq技術(shù)的特點

RNA-seq對芯片的優(yōu)勢

RNA-seq檢測基因表達主要在7個方面比基因芯片有優(yōu)勢寻行。

RNA-seq存在的問題

• RNA-seq測序之前需要一個比較復(fù)雜的文庫構(gòu)建過程，這個過程的每一步都可能帶來誤差甚至導(dǎo)致實驗失敗匾荆。如cDNA片段化拌蜘、PCR擴增等都會帶來偏倚，最終導(dǎo)致有的片段被反復(fù)測了多次牙丽，有的沒有測到简卧。rRNA去除不干凈等因素也會帶來大量污染（即非目標序列）。還有其他一些實驗問題烤芦。
• RNA-seq檢測靈敏度和最大值是隨測序深度變化的举娩，深度不夠，不能發(fā)現(xiàn)超低表達的轉(zhuǎn)錄本，需要在測序前預(yù)估轉(zhuǎn)錄本大小晓铆。而由于復(fù)雜的RNA編輯等原因，高等生物的轉(zhuǎn)錄組與其編碼基因數(shù)量沒有固定比例關(guān)系绰播，所以預(yù)估容易產(chǎn)生較大誤差骄噪。
• 參考基因組或轉(zhuǎn)錄組不準確、測序誤差蠢箩、錯誤拼接或比對帶來的錯誤會大大影響各種變異或者可變剪切事件的識別链蕊。
• 各種其他的問題。比如整個實驗流程可能引進了各種污染谬泌；原始數(shù)據(jù)預(yù)處理的數(shù)據(jù)模型不完善等等滔韵。

下面看一下如何計算由測序誤差引起的Barcode的錯誤分配，假設(shè)Barcode（Barcode唯一確定樣本）長度為6個堿基掌实，每個Barcode兩兩之間兩個堿基不同陪蜻，所有的Barcode都用滿，同時假設(shè)錯誤的發(fā)生符合二項分布贱鼻，那么只要2個堿基錯誤宴卖，就會發(fā)生一次錯誤分配，在Illumina測序儀每個堿基的平均錯誤率0.5%的前提下邻悬，下面例子可以計算出一個泳道的錯誤錯誤分配概率症昏。

> p = 0.0005
> sum(sapply(2:6, FUN=function(k) choose(6,k)*p^k*(1-p)^(6-k)))
[1] 3.745003e-06

也就是說，在一個泳道父丰，每百萬讀段就會有370個讀段分配錯誤肝谭。如果還考慮DNA簇混合和跳躍PCR引起的Barcode錯誤分配，這個數(shù)值還要高很多蛾扇。這種分配錯誤對一般的轉(zhuǎn)錄組分析沒有影響攘烛，但是對一些高靈敏度的突變檢測項目影響很大。

Illumine 2000測序儀中屁桑，一次運行（Run）可以使用2個流動槽（Flow cell）医寿，每個流動槽包括8個泳道（Lane），一個泳道包含2個面（Surface）蘑斧，每個面還有3個條（Swath）也叫列（Column）靖秩，每一列由16個小區(qū)（Tile）組成，后者又由大量DNA簇（Cluster）組成竖瘾。

Illumine 2000測序儀每次運行（單端測序）理論上可以產(chǎn)生大約30億個DNA簇沟突，每個DNA簇理論上可以產(chǎn)生一條讀段（Read），如果測序長度為100bp捕传，一次運行可以得到3G個讀段惠拭，其原始數(shù)據(jù)量為300G個堿基。

RNA-seq數(shù)據(jù)預(yù)處理

RNA-seq數(shù)據(jù)預(yù)處理與基因芯片預(yù)處理的目的一致，都是要得到基因表達數(shù)據(jù)职辅，這里確切說是轉(zhuǎn)錄本棒呛。但它們的實現(xiàn)細節(jié)和方法則有很大不同，下面逐一介紹域携。

質(zhì)量控制

數(shù)據(jù)質(zhì)量分析報告可以調(diào)用ShortRead包中的qa函數(shù)得到簇秒，另一個常用的工具是FastQC。

下面是一個用qa函數(shù)生成質(zhì)量分析報告的實例秀鞭，示例數(shù)據(jù)來自SRA數(shù)據(jù)庫的RNA-seq數(shù)據(jù)SRR921344趋观。

library(BiocInstaller)
# 安裝包
biocLite("SRAdb")
biocLite("R.utils")
# 導(dǎo)入包
library(ShortRead)
library(SRAdb)
library(R.utils)
# download data
db_dir <- "~/Projects/coGECP-pro/db/"
sra_dbname <- paste0(db_dir, 'SRAmetadb.sqlite')    
sra_con <- dbConnect( dbDriver("SQLite"), sra_dbname )
getFASTQfile("SRR921344", sra_con, destDir = "~/Projects/coGECP-pro/fast-format/", srcType = 'ftp', ascpCMD = NULL )
dbDisconnect( sra_con )

# 解壓后改名
gunzip("../fast-format/SRR921344.fastq.gz", destname="../fast-format/T2-1.fastq")

# 需要分析的數(shù)據(jù)文件名稱
fastqfile="T2-1.fastq"
# 得到質(zhì)量分析的記過
qa <- qa(dirPath = "../fast-format/", pattern=fastqfile, type="fastq")

# 輸出html格式的分析報告
report(qa, dest="qcReport", type="html")

執(zhí)行完畢后會得到名為qcReport的結(jié)果目錄，主要包括質(zhì)量分析的報告文件“index.html”和文件夾“image”锋边，image文件夾下包括所有的統(tǒng)計信息圖示皱坛。

比較重要的圖示信息有前20個高頻出現(xiàn)的讀段統(tǒng)計信息，這對于確定接頭或者其他污染很重要豆巨。FastQC軟件帶有一個文件包括了常用的Illumina接頭序列剩辟，會把高頻出現(xiàn)的序列比對到這些接頭序列，并給予提示搀矫；Biocondutor則需要編程比對到NCBI UniVec數(shù)據(jù)庫（http://www.ncbi.nlm.nih.gov/tools/vecscreen/univec）來確定接頭序列抹沪。

接下來是四種堿基的逐點質(zhì)量圖，該圖橫坐標是測序的循環(huán)數(shù)瓤球，對應(yīng)測序時5'端開始的每個位點融欧，縱坐標是所有該位點測量的堿基的質(zhì)量分數(shù)平均數(shù)。

質(zhì)量控制中最重要的一個圖是逐點質(zhì)量圖卦羡，它與四種堿基逐點質(zhì)量圖的區(qū)別在于不區(qū)分堿基種類噪馏，給出每個位點的測序質(zhì)量分數(shù)的平均數(shù)、中位數(shù)和上下四分位線绿饵。

通過質(zhì)量控制欠肾，可以確定當(dāng)前樣品的數(shù)據(jù)是否應(yīng)該保留進入下一步分析或者丟棄。對于通過質(zhì)量控制的數(shù)據(jù)拟赊，還要進行讀段處理刺桃，清理后的數(shù)據(jù)才是實際分析中使用的數(shù)據(jù)。

讀段清理

讀段清理主要去掉讀段中多個“N”堿基吸祟，讀段兩端的低質(zhì)量區(qū)域（質(zhì)量分數(shù)少于Q20）瑟慈，讀段3'端可能混入的接頭序列，還有可能污染進來的rRNA和病毒序列屋匕。

轉(zhuǎn)錄組組裝

轉(zhuǎn)錄組組裝包括從頭組裝和基于參考序列的組裝葛碧。從頭組裝最常用的軟件是Trinity，其優(yōu)點是：不依賴參考序列过吻；能較好地重建變異的进泼、可變剪切的或者來自染色體重組的轉(zhuǎn)錄本。從頭組裝的缺點是會消耗大量的內(nèi)存資源，測序深度的需求也很高乳绕，對測序錯誤很敏感绞惦，高相似度的轉(zhuǎn)錄本可能會被誤拼到一起。基于參考序列的組裝常使用TopHat和Cufflinks（主要用于轉(zhuǎn)錄本的識別洋措、定量翩隧、標準化與差異分析）兩種軟件的組合來完成。其優(yōu)點是：內(nèi)存需求猩胛啤；污染影響小专缠，因為污染讀段不能比對到參考序列雷酪；靈敏度高，所需測序深度低涝婉，能檢測低豐度的轉(zhuǎn)錄本哥力；組裝的轉(zhuǎn)錄本序列更完整；可以增加參考基因組中的轉(zhuǎn)錄本注釋墩弯》园希基于參考序列的組裝的缺點：嚴重依賴參考序列及其注釋信息等。

轉(zhuǎn)錄組定量和標準化

有參的轉(zhuǎn)錄組定量渔工，需要利用TopHat比對軟件將所有讀段比對到參考基因組上锌钮，然后由Cuffinks軟件完成定量；無參的轉(zhuǎn)錄組定量需要利用Bowtie或BWA比對軟件將所有讀段比對到組裝得到的轉(zhuǎn)錄組上直接計數(shù)引矩。

單端測序讀段的比對比較簡單梁丘，雙端測序的質(zhì)量不一致，往往是反向一端測序質(zhì)量低旺韭，如果按照同樣的標準要求兩端測序的讀段都比對上氛谜，會丟失很多比對結(jié)果。一般采取的方式是兩端讀段分別依據(jù)不同的標準（例如正向允許錯配一個区端，反向允許錯配兩個）做單端比對值漫，然后根據(jù)兩端對齊后中間距離抽取成對的比對結(jié)果。

轉(zhuǎn)錄組定量得到的基因表達矩陣是簡單計數(shù)得來的织盼，因此稱作原始計數(shù)（Raw counts）杨何。有參的轉(zhuǎn)錄組定量，可以用Bioconductor的GenomicRanges/Rsamtools軟件包中的summarizeOverlaps函數(shù)實現(xiàn)悔政。輸入的數(shù)據(jù)為比對后得到的Bam或者Sam文件晚吞，經(jīng)過基因水平或者外顯子水平的計數(shù)，可以直接輸出為某種預(yù)定義對象谋国，便于下游軟件（如edgeR包和DESeq包）繼續(xù)處理槽地。

下面例子使用summarizeOverlaps函數(shù)從Bam文件中獲取原始計數(shù)，并輸出為edgeR包和DESeq包中的數(shù)據(jù)對象。

require(BiocInstaller)
# biocLite("TxDb.Dmelanogaster.UCSC.dm3.ensGene")
# biocLite("DESeq")
require(Rsamtools)
require(DESeq)
require(edgeR)
require(pasillaBamSubset)
library(GenomicAlignments)

#此處數(shù)據(jù)已經(jīng)不能用了捌蚊，所以從其他包里弄了個bam格式文件試試
# bamfile <- system.file("extdata", "ex1.bam", package="Rsamtools")
# bf1 <- BamFileList(bamfile, index=character())
# features <- GRanges(seqnames = c())

# from GenomicRanges::GenomicRangesHOWTOs   
reads <- c(untrt1=untreated1_chr4(),
             untrt3=untreated3_chr4())
# single-end reads
# paired-end reads

library(TxDb.Dmelanogaster.UCSC.dm3.ensGene)
exbygene <- exonsBy(TxDb.Dmelanogaster.UCSC.dm3.ensGene, "gene")
se <- summarizeOverlaps(exbygene, reads, mode="IntersectionNotEmpty")

# back 
deseq <- newCountDataSet(assays(se)$counts, rownames(colData(se)))
edger <- DGEList(assays(se)$counts, group=rownames(colData(se)))

類似基因芯片集畅，RNA-seq定量后的數(shù)據(jù)需要標準化，使得所有的樣品具有可比性缅糟。最常見的一個指標是RPKM（Reads Per Kilo base per Million reads）挺智，即每百萬讀段中來自某一個基因每千堿基長度的讀段數(shù)目。具體計算時使用比對到某個基因的讀段個數(shù)除以比對到基因組或者轉(zhuǎn)錄組的所有讀段個數(shù)（以百萬為單位）窗宦，再除以基因或轉(zhuǎn)錄本的長度（以KB為單位）赦颇。基因表達差異分析的常用軟件edgeR軟件和DESeq都自帶了數(shù)據(jù)標準化功能赴涵，可以直接處理用原始計數(shù)表示的基因表達矩陣媒怯，RPKM表示的基因表達矩陣的使用主要是為了方便用戶直接觀察數(shù)據(jù)。

RNA-seq數(shù)據(jù)分析

RNA-seq數(shù)據(jù)與表達譜芯片數(shù)據(jù)在基因表達差異的顯著性分析流程基本相同髓窜，不同的地方只在確定差異表達基因方面扇苞。

基因表達差異的顯著性分析

表達差異分析只對比不同樣本之間的同一個轉(zhuǎn)錄本，所以不需要考慮轉(zhuǎn)錄本長度寄纵，只考慮總讀段數(shù)鳖敷。一個最簡單思想就是，樣本測序得到的總讀段數(shù)（實際上是可以比對到轉(zhuǎn)錄組的總讀段數(shù)）越多程拭，則每個基因分配到的讀段越多定踱。因此最簡單的標準化因子就是總讀段數(shù)，用總讀段數(shù)作標準化的前提是大部分基因的表達是非顯著變化的恃鞋，這與基因芯片中的基本假設(shè)相同屋吨。但是實際工作中發(fā)現(xiàn)很多情況下總讀段數(shù)主要是一小部分大量表達的基因貢獻的。Bullard等（2010）在比較了幾種標準化方法的基礎(chǔ)上發(fā)現(xiàn)在每個泳道內(nèi)使用非零計數(shù)分布的上四分位數(shù)（Q75%）作為標準化因子是一種更穩(wěn)健的選擇山宾，總體表現(xiàn)是所研究方法中最優(yōu)的至扰。

Bioconductor的edgeR包和DESeq包分別提供了上四分位數(shù)和中位數(shù)來作為標準化因子，就是出于這個思想资锰。

edgeR提供了三種標準化算法敢课，分別是M值加權(quán)截斷均值法（Weighted trimmed mean of M-values, TMM），相對對數(shù)表達值法（Relative log expression, RLE）和上四分位法（Upperquartile）绷杜，其中TMM是默認設(shè)定直秆。這些標準化方法大同小異，其基本思想就是去除表達值較大的少數(shù)基因的影響鞭盟，而保留大部分沒有顯著變化的基因圾结。

由于基因芯片檢測雜交的熒光強度信號是連續(xù)值，往往假設(shè)它符合正態(tài)分布齿诉；而RNA-seq測量的是離散值筝野，最簡單的假設(shè)就是二項分布晌姚。由于RNA-seq讀段數(shù)量非常多，而且一條讀段映射到一個給定基因的概率足夠小歇竟，在實際計算中挥唠，二項分布常用它的極限形式泊松分布來替代。泊松分布的一個性質(zhì)是其方差等于均值焕议，但是當(dāng)有生物學(xué)重復(fù)時宝磨，RNA-seq數(shù)據(jù)會表現(xiàn)出比泊松分布期望的更高的變異性，對相當(dāng)多的基因來說方差可能超過均值盅安，這種現(xiàn)象叫過離散唤锉。對過離散數(shù)據(jù)，基于泊松分布假設(shè)的分析容易低估不同生物學(xué)重復(fù)帶來的取樣誤差而得到過多的假陽性的差異表達基因别瞭。為了允許額外的變異腌紧，一個自然的想法就是給均值加上一個散度參數(shù)，以使方差可以大于均值畜隶。于是作為泊松分布的推廣，又引入了負二項分布（NB）來作為基本假設(shè)号胚，負二項分布是當(dāng)前基因表達的顯著性分析中最常用假設(shè)籽慢。

現(xiàn)在從百度百科摘取二項分布與負二項分布的公式，以加強理解猫胁。

二項分布箱亿，即重復(fù)n次的貝努利試驗，用表示隨機試驗的結(jié)果弃秆。如果事件發(fā)生的概率是p,不發(fā)生的概率是1-p届惋，那么N次獨立重復(fù)試驗中發(fā)生k次的概率是：

期望：

方差：

而負二項分布的公式（概率密度）為：

它表示，已知一個事件在貝努利試驗中每次的出現(xiàn)概率是p菠赚，在一連串貝努利試驗中脑豹，一件事剛好在第r+k次試驗出現(xiàn)第r次的概率。

寫作：X ~ NB(r; p)

基于負二項分布的edgeR衡查、DESeq包是當(dāng)前最主要的分析程序瘩欺。在edgeR中，對于任一樣品i中的任一個基因g拌牲，假設(shè)它的分布符合負二項分布俱饿。

其中是樣品i中的讀段總數(shù)（實際中是可以比對上的讀段總數(shù)）；就是基因g的散度塌忽；是基因g在某個實驗條件j下或者分組j中的相對豐度拍埠。第g個基因在某個實驗條件j下或分組j中，NB分布的均值為土居，方差為枣购，對于表達差異分析嬉探，需要估計的是散度，當(dāng)它趨于0時坷虑，負二項分布退化為泊松分布甲馋，這時方差退化為第一項，一般認為來自技術(shù)重復(fù)迄损，方差第2項來自生物學(xué)重復(fù)定躏。

RNA-seq數(shù)據(jù)分析往往只有很小的樣本量，為每個基因估計一個散度非常困難芹敌。比較好的策略是允許不同的基因有不同的個體散度痊远，而這些個體散度的估計可以用一些合適的統(tǒng)計方法借助基因間的信息來改進。相對于edgeR氏捞，DESeq默認設(shè)置采取了最保守的估計策略碧聪，即選取每個基因的經(jīng)驗散度和擬合得到的散度趨勢線取值中最大的作為最終的散度估計值，因此DESeq往往選出更少的差異表達基因（為什么呢液茎？）逞姿。DESeq由于可以利用同一個樣本基因間的數(shù)據(jù)估計散度，而不一定需要重復(fù)樣本來計算捆等，因此可以直接用于無重復(fù)實驗的表達差異分析滞造。

后面一些實例代碼不好重復(fù)，不再描述栋烤。