https://hbctraining.github.io/DGE_workshop/lessons/01_DGE_setup_and_overview.html
（今天狀態(tài)不佳悄窃，加上統(tǒng)計(jì)學(xué)基礎(chǔ)不怎么樣授霸，如有錯誤請指出狸窘，我將在簡書中更新）
學(xué)習(xí)目標(biāo)
- 解釋實(shí)驗(yàn)及其目標(biāo)
- 描述如何在R中設(shè)置RNA-seq項(xiàng)目
- 描述RNA-seq和差異基因表達(dá)分析工作流程
- 解釋為什么負(fù)二項(xiàng)分布用于模擬RNA-seq計(jì)數(shù)數(shù)據(jù)
差異基因表達(dá)（DGE）分析概述
RNA-seq的目標(biāo)通常是進(jìn)行差異表達(dá)檢測袁辈，以確定哪些基因在不同條件間表達(dá)量有差異恢口。這些基因可以從生物學(xué)角度揭示不同條件下受影響的生物過程绩鸣。
RNA-seq工作流程的詳細(xì)步驟如下圖虏冻，可用于確定基因的表達(dá)水平槽棍。通過生成每個基因的read counts墓怀，在命令行（Linux / Unix）上執(zhí)行所有步驟纳账。差異表達(dá)分析和其他下游功能分析通常用R中的專門的程序包完成，這些R包是為完成差異分析所需的復(fù)雜統(tǒng)計(jì)分析而設(shè)計(jì)的捺疼。
在接下來的幾節(jié)課中疏虫，我們將引導(dǎo)你使用各種R包進(jìn)行end-to-end的基因水平RNA-seq差異表達(dá)工作流程。我們將從計(jì)數(shù)矩陣開始啤呼，進(jìn)行質(zhì)控的探索性數(shù)據(jù)分析卧秘，探索樣本之間的關(guān)系，進(jìn)行差異表達(dá)分析官扣，并在執(zhí)行下游功能分析之前直觀地探索結(jié)果翅敌。
1.檢查示例數(shù)據(jù)集
示例數(shù)據(jù)集是RNA-Seq的完整計(jì)數(shù)矩陣，是Kenny PJ et al惕蹄，Cell Rep 2014(http://www.ncbi.nlm.nih.gov/pubmed/25464849)描述的更大研究的一部分蚯涮。
RNA-Seq實(shí)驗(yàn)對象是HEK293F細(xì)胞，所述HEK293F細(xì)胞用MOV10轉(zhuǎn)基因或siRNA轉(zhuǎn)染以降低Mov10表達(dá)卖陵，或非特異性（無關(guān)）siRNA轉(zhuǎn)染遭顶。這導(dǎo)致了3種情況：Mov10 oe（過表達(dá)），Mov10 kd（擊倒）和不相關(guān)的kd泪蔫。重復(fù)次數(shù)如下所示棒旗。
使用這些數(shù)據(jù)，我們將評估與MOV10表達(dá)的擾動相關(guān)的轉(zhuǎn)錄模式撩荣。請注意铣揉，不相關(guān)的siRNA作為對照組。
這些數(shù)據(jù)集的目的是什么餐曹？Mov10做了什么逛拱？
作者正在研究脆性X綜合征中涉及的各種基因之間的相互作用，這是一種FMRP蛋白異常產(chǎn)生的疾病台猴。
FMRP “最常見于大腦朽合，對正常的認(rèn)知發(fā)育和女性生殖功能至關(guān)重要额湘。該基因的突變可導(dǎo)致脆性X綜合征，精神發(fā)育遲滯旁舰，卵巢早衰，自閉癥嗡官，帕金森病箭窜，發(fā)育遲緩和其他認(rèn)知缺陷⊙苄龋“ - 來自維基百科
MOV10是推定的RNA解旋酶磺樱，其在microRNA途徑的背景下也與FMRP相關(guān)。
該論文正在測試的假說是FMRP和MOV10結(jié)合并調(diào)節(jié)RNA子集的翻譯婆咸。
問題：
- 可以通過MOV10的丟失或獲得來識別哪些表達(dá)模式竹捉？
- 兩種情況之間是否有共同的基因？
2.組織工作目錄
在深入了解分析的細(xì)節(jié)之前尚骄，先打開RStudio并為此分析設(shè)置一個新項(xiàng)目块差。
- 轉(zhuǎn)到
File菜單并選擇New Project。 - 在
New Project窗口中倔丈，選擇New Directory憨闰。然后，選擇Empty Project需五。為新目錄命名DEanalysis鹉动，然后將“創(chuàng)建項(xiàng)目作為子目錄：”桌面（或自定義的位置）。 - 新的project應(yīng)該在RStudio中自動打開了宏邮。
要檢查是否在正確的工作目錄中泽示，請使用getwd()。在控制臺中應(yīng)該會返回路徑Desktop/DEanalysis蜜氨。在工作目錄下使用New folder按鈕創(chuàng)建三個新目錄：data械筛，meta和results。請記住飒炎，好的分析的關(guān)鍵是從一開始就保持井井有條变姨！
轉(zhuǎn)到File菜單并選擇New File，然后選擇R Script厌丑。此時在Rstudio左上角打開了一個腳本編輯器定欧。這是接下來輸入和保存此分析所需的所有命令的地方。在腳本編輯器中輸入標(biāo)題行：
## Gene-level differential expression analysis using DESeq2
現(xiàn)在將文件另存為de_script.R怒竿。完成后砍鸠，工作目錄應(yīng)如下圖：
最后，需要獲取我們將用于分析的文件耕驰。右鍵單擊下面的鏈接爷辱，然后選擇“將鏈接另存為...”選項(xiàng)進(jìn)行下載：
- 將counts矩陣(https://raw.githubusercontent.com/hbc/NGS_Data_Analysis_Course/master/sessionIII/data/Mov10_full_counts.txt)文件保存在
data目錄中。 - 將metadata數(shù)據(jù)表(https://raw.githubusercontent.com/hbc/NGS_Data_Analysis_Course/master/sessionIII/data/Mov10_full_meta.txt)文件保存在
meta目錄中。
4. 加載包
對于這個分析用到的R包饭弓，有的是從CRAN安裝的双饥，有的是從Bioconductor安裝的。要使用這些包（以及它們中包含的函數(shù)）弟断，我們需要加載包咏花。將以下內(nèi)容添加到腳本中，可增加自主注釋
## Setup
### Bioconductor and CRAN libraries used
library(tidyverse)
library(RColorBrewer)
library(DESeq2)
library(pheatmap)
library(DEGreport)
5. 載入數(shù)據(jù)
要將數(shù)據(jù)加載到我們當(dāng)前的環(huán)境中阀趴，我們將使用該read.table功能昏翰。我們需要提供每個文件的路徑，并指定參數(shù)讓R知道我們有一個header（header = T）刘急，第一列是我們的行名（row.names =1）棚菊。默認(rèn)情況下，該函數(shù)需要使用制表符分隔的文件叔汁，這就是我們所擁有的统求。
## Load in data
data <- read.table("data/Mov10_full_counts.txt", header=T, row.names=1)
meta <- read.table("meta/Mov10_full_meta.txt", header=T, row.names=1)
使用class()來檢查我們的數(shù)據(jù)：
### Check classes of the data we just brought in
class(meta)
class(data)
6. 查看數(shù)據(jù)
在繼續(xù)執(zhí)行任何類型的分析之前，請確保數(shù)據(jù)集包含預(yù)期的樣本/信息据块。
View(meta)
View(data)
使用Salmon的豐度估計(jì)值作為DESeq2的輸入
課程中使用的counts是使用RNA-seq分析的標(biāo)準(zhǔn)方法生成的球订，其中使用剪接感知對準(zhǔn)器將樣本與基因組比對，然后計(jì)數(shù)瑰钮。如果使用輕量級算法（如Salmon冒滩，Sailfish或Kallisto）估算豐度，還可以使用DESeq2執(zhí)行基因水平差異表達(dá)分析浪谴。這些轉(zhuǎn)錄本豐度估計(jì)值（通常稱為“偽計(jì)數(shù)”）可以轉(zhuǎn)換為與DESeq2一起使用开睡，但設(shè)置稍微復(fù)雜一些。有關(guān)使用DESeq2的Salmon偽載體的更多信息苟耻，詳見https://hbctraining.github.io/DGE_workshop_salmon/lessons/01_DGE_setup_and_overview.html篇恒。
那么counts數(shù)據(jù)實(shí)際代表什么呢？用于差異表達(dá)分析的counts數(shù)據(jù)表示源自特定基因的序列reads的數(shù)量凶杖。counts越多胁艰，與該基因相關(guān)的reads數(shù)越多，就假設(shè)樣本中該基因的表達(dá)水平越高智蝠。
差異表達(dá)分析可以尋找兩個或更多組（在meta數(shù)據(jù)中定義）之間表達(dá)量發(fā)生變化的基因
- 處理與對照
- 表達(dá)與某些變量或臨床結(jié)果的相關(guān)性
為什么不能根據(jù)兩組之間的差異（基于倍數(shù)變化值）對基因進(jìn)行排序腾么，來識別差異表達(dá)的基因？
通常情況下杈湾，得到的數(shù)據(jù)會比預(yù)期的要多得多解虱。樣品之間表達(dá)水平不同的基因不僅是感興趣的實(shí)驗(yàn)變量的結(jié)果，而且與一些外界因素有關(guān)漆撞。差異表達(dá)分析的目的是確定這些因素的相對作用殴泰，并將“interesting”與“ uninteresting”的因素分開于宙。
“ uninteresting”的因素也會導(dǎo)致數(shù)據(jù)的差異，因此即使樣本組之間的平均表達(dá)水平看起來相差很大悍汛，但實(shí)際上可能差異并不顯著捞魁。這是針對下圖中“GeneA”的“untreated”組和“treated”之間表達(dá)進(jìn)行說明的。'treated'組的geneA的平均表達(dá)水平是'untreated'組的兩倍离咐，但重復(fù)之間的差異表明谱俭，這可能不是顯著差異。在確定基因是否差異表達(dá)時健霹，需要考慮數(shù)據(jù)的變化（以及它可能來自何處）。
差異表達(dá)分析的目標(biāo)是針對每個基因確定組間表達(dá)（計(jì)數(shù)）的差異是否顯著瓶蚂，給定組內(nèi)（重復(fù)）觀察到的變異量糖埋。為了測試顯著性，我們需要一個適當(dāng)?shù)慕y(tǒng)計(jì)模型窃这，準(zhǔn)確地執(zhí)行標(biāo)準(zhǔn)化（以解釋測序深度的差異等）和方差建模（以考慮少數(shù)重復(fù)和大動態(tài)表達(dá)范圍）瞳别。
RNA-seq計(jì)數(shù)分布
為了確定適當(dāng)?shù)慕y(tǒng)計(jì)模型，需要有關(guān)計(jì)數(shù)分布的信息杭攻。繪制樣本'Mov10_oe_1'的計(jì)數(shù)分布圖：
ggplot(data) +
geom_histogram(aes(x = Mov10_oe_1), stat = "bin", bins = 200) +
xlab("Raw expression counts") +
ylab("Number of genes")
如果我們放大到接近零祟敛，我們可以看到大量基數(shù)為零：
ggplot(data) +
geom_histogram(aes(x = Mov10_oe_1), stat = "bin", bins = 200) +
xlim(-5, 500) +
xlab("Raw expression counts") +
ylab("Number of genes")
上圖說明了RNA-seq計(jì)數(shù)數(shù)據(jù)的一些普遍特征：大部分基因計(jì)數(shù)很少，以及由于表達(dá)量低到圖上不顯示而具有很長的右尾兆解。與芯片數(shù)據(jù)不同馆铁，芯片數(shù)據(jù)由于探針有最大限度，counts的動態(tài)范圍有最大值锅睛，而RNA-seq數(shù)據(jù)沒有埠巨。由于這些技術(shù)的不同，RNA-Seq和芯片用于擬合數(shù)據(jù)的統(tǒng)計(jì)模型也是不同的现拒。
注意：芯片數(shù)據(jù)的對數(shù)強(qiáng)度接近正態(tài)分布辣垒。而由于RNA-seq計(jì)數(shù)數(shù)據(jù)的不同性質(zhì)，例如整數(shù)計(jì)數(shù)而不是連續(xù)測量和非正態(tài)分布數(shù)據(jù)印蔬，正態(tài)分布不能準(zhǔn)確地模擬RNA-seq計(jì)數(shù)[ 1 ]勋桶。
7.counts數(shù)據(jù)建模
計(jì)數(shù)數(shù)據(jù)通常使用二項(xiàng)分布建模，這可以讓你有可能在投擲硬幣多次時獲得多個heads（不會翻譯）侥猬。但是例驹，并非所有計(jì)數(shù)數(shù)據(jù)都適合二項(xiàng)分布。二項(xiàng)式基于離散事件退唠，有確定數(shù)量的樣本時使用眠饮。
樣本（即購買彩票的人）數(shù)量非常大，但事件（獲勝的概率）的概率非常小铜邮，泊松分布用于模擬這些類型的計(jì)數(shù)數(shù)據(jù)仪召。泊松與二項(xiàng)式類似寨蹋，但基于連續(xù)事件。詳見Rafael Irizarry在EdX課程(https://youtu.be/fxtB8c3u6l8)
使用RNA-Seq數(shù)據(jù)扔茅，有非常多的RNA已旧，并且抽出特定轉(zhuǎn)錄本的概率非常小。因此召娜，使用泊松分布適合這種情況运褪。但是，此分布的唯一屬性是均值==方差玖瘸。實(shí)際上秸讹，對于RNA-Seq數(shù)據(jù)，重復(fù)中總是存在一些生物學(xué)變異（在樣本類中）雅倒。具有較大平均表達(dá)水平的基因傾向于在重復(fù)中具有較大的觀察到的變異璃诀。
如果mRNA的比例在每個樣本類別的生物學(xué)重復(fù)之間保持完全恒定，我們可以預(yù)期泊松分布（其中均值==方差）蔑匣。Rafael Irizarry在EdX課程(https://youtu.be/HK7WKsL3c2w)中對這個概念進(jìn)行了很好的描述劣欢。但這在實(shí)踐中不會發(fā)生，因此泊松分布被認(rèn)為僅適用于單個生物樣本裁良。
考慮到重復(fù)之間的這種類型的可變性凿将，最適合的模型是負(fù)二項(xiàng)（NB）模型。從本質(zhì)上講价脾，NB模型是平均值<方差的數(shù)據(jù)的良好擬合牧抵，就像RNA-Seq計(jì)數(shù)數(shù)據(jù)的情況一樣。
注意：
- 生物學(xué)重復(fù)代表相同樣品類別的多個樣品（即來自不同小鼠的RNA）
- 技術(shù)重復(fù)代表相同的樣品（即來自相同小鼠的RNA）侨把，但復(fù)制了技術(shù)步驟
- 通常生物學(xué)方差遠(yuǎn)大于技術(shù)方差灭忠，因此我們不需要考慮技術(shù)方差來識別表達(dá)中的生物學(xué)差異
- 不要在技術(shù)重復(fù)上花錢 - 生物復(fù)制更有用
注意： 如果使用細(xì)胞系并且不確定是否準(zhǔn)備了生物學(xué)或技術(shù)重復(fù)，請查看此鏈接座硕。這是一個有用的資源弛作，可幫助確定如何最好地設(shè)置體外實(shí)驗(yàn)。
如何知道我的數(shù)據(jù)是否應(yīng)使用泊松分布或負(fù)二項(xiàng)分布建模华匾？
如果它是計(jì)數(shù)數(shù)據(jù)映琳，它應(yīng)該適合負(fù)二項(xiàng)式蜘拉，如前所述萨西。但是旭旭，繪制平均值與數(shù)據(jù)方差的關(guān)系可能會有所幫助。記住泊松模型持寄，均值=方差源梭，但負(fù)二項(xiàng)分布的NB均值<方差。
運(yùn)行以下代碼以繪制'Mov10過表達(dá)'重復(fù)的平均值與方差：
mean_counts <- apply(data[, 3:5], 1, mean)
variance_counts <- apply(data[, 3:5], 1, var)
df <- data.frame(mean_counts, variance_counts)
ggplot(df) +
geom_point(aes(x=mean_counts, y=variance_counts)) +
geom_line(aes(x=mean_counts, y=mean_counts, color="red")) +
scale_y_log10() +
scale_x_log10()
注意荠卷，在上圖中，重復(fù)的變異傾向于大于平均值（紅線）烛愧，特別是對于具有大平均表達(dá)水平的基因油宜。這是一個很好的跡象，表明我們的數(shù)據(jù)不符合泊松分布慎冤，我們需要使用負(fù)二項(xiàng)模型來解釋這種方差的增加（即泊松會低估導(dǎo)致假陽性DE基因增加的變異性）沧卢。
通過生物學(xué)重復(fù)改善平均值估計(jì)（即減少方差）
隨著生物學(xué)重復(fù)的數(shù)量（隨著重復(fù)次數(shù)的增加，分布將接近泊松分布）违寿，方差或散射傾向于減少熟空，因?yàn)槠骄档臉?biāo)準(zhǔn)偏差小于個體觀察的標(biāo)準(zhǔn)偏差搞莺。額外重復(fù)的價值在于，當(dāng)添加更多數(shù)據(jù)（重復(fù)）時才沧，可以獲得更精確的組均值估計(jì)温圆，并最終對區(qū)分樣本類別（即更多DE基因）之間差異的能力更有信心。
下圖說明了測序深度與重復(fù)數(shù)量之間的關(guān)系岁歉，鑒定了差異表達(dá)基因的數(shù)量 (https://academic.oup.com/bioinformatics/article/30/3/301/228651/RNA-seq-differential-expression-studies-more) 锅移。注意，與增加測序深度相比非剃，重復(fù)數(shù)量的增加傾向于返回更多的DE基因。因此券坞，通常增加重復(fù)比增加測序深度更好，但需要注意的是深浮，檢測低表達(dá)的差異基因和進(jìn)行同種型水平的差異表達(dá)需要更高的測序深度眠冈。一般來說，建議的最小測序深度是每個樣品20-30萬個reads蜗顽，但如果有大量重復(fù)布卡，我們已經(jīng)看到了具有1000萬個reads的優(yōu)秀RNA-seq實(shí)驗(yàn)雇盖。
8.差異表達(dá)分析工作流程
為了在進(jìn)行差異表達(dá)分析時適當(dāng)?shù)貙τ?jì)數(shù)進(jìn)行建模崔挖，已經(jīng)開發(fā)了許多用于RNA-seq數(shù)據(jù)的差異表達(dá)分析的軟件包。即使在不斷開發(fā)新方法的過程中薛匪，通常也會建議使用一些工具作為最佳實(shí)踐，例如DESeq2和EdgeR逸尖。這兩種工具都使用負(fù)二項(xiàng)模型瘸右，采用類似的方法，通常產(chǎn)生類似的結(jié)果苞俘。它們非常嚴(yán)格龄章，并且在靈敏度和特異性之間取得了良好的平衡（減少誤報和假陰性）。
Limma-Voom是另一組常用于DE分析的工具瓦堵，但這種方法對于小樣本量可能不太敏感菇用。當(dāng)每組的生物學(xué)重復(fù)數(shù)量變大（> 20）時，推薦使用該方法惋鸥。
許多描述這些方法之間比較的研究表明悍缠，盡管存在一些一致性耐量，但工具之間也存在很大差異廊蜒。此外，沒有一種方法在所有條件下都能達(dá)到最佳效果（Soneson和Dleorenzi山叮，2013）。
本課程將使用DESeq2進(jìn)行DE分析脑又，使用DESeq2的分析步驟將在下面的流程圖中以綠色顯示问麸。DESeq2首先將計(jì)數(shù)數(shù)據(jù)標(biāo)準(zhǔn)化以解釋樣品之間文庫大小和RNA組成的差異钞翔。然后用標(biāo)準(zhǔn)化計(jì)數(shù)在基因和樣本水平上制作QC的一些圖。最后一步是使用DESeq2包中的相應(yīng)函數(shù)來執(zhí)行差異表達(dá)分析嗅战。
接下來的課程中會深入介紹這些步驟驮捍，但有關(guān)DESeq2的更多詳細(xì)信息和有用的建議可以在DESeq2的vignette中(http://bioconductor.org/packages/devel/bioc/vignettes/DESeq2/inst/doc/DESeq2.html)找到脚曾。如果完成此工作流程出現(xiàn)問題，可參考RStudio中的vignette：
vignette("DESeq2")
通過這個命令很方便地提供了觸手可及的豐富信息珊泳！課程中如有需要可通過上述命令調(diào)用它拷沸。
