|
差異表達(dá)分析之FDR
隨著測(cè)序成本的不斷降低，轉(zhuǎn)錄組測(cè)序分析已逐漸成為一種很常用的分析手段。但對(duì)于轉(zhuǎn)錄組分析當(dāng)中的一些概念，很多人還不是很清楚。今天，小編就來談?wù)勗谵D(zhuǎn)錄組分析中，經(jīng)常會(huì)遇到的一個(gè)概念FDR，那什么是FDR?為什么要用FDR呢?一起來學(xué)習(xí)吧！ 什么是FDR FDR (false discovery rate)，中文一般譯作錯(cuò)誤發(fā)現(xiàn)率。在轉(zhuǎn)錄組分析中，主要用在差異表達(dá)基因的分析中，控制最終分析結(jié)果中，假陽性結(jié)果的比例。 為什么要用FDR 在轉(zhuǎn)錄組分析中，如何確定某個(gè)轉(zhuǎn)錄本在不同的樣品中表達(dá)量是否有差異是分析的核心內(nèi)容之一。一般來說，我們認(rèn)為，不同樣品中，表達(dá)量差異在兩倍以上的轉(zhuǎn)錄本，是具有表達(dá)差異的轉(zhuǎn)錄本。為了判斷兩個(gè)樣品之間的表達(dá)量差異究竟是由于各種誤差導(dǎo)致的還是本質(zhì)差異，我們需要根據(jù)所有基因在這兩個(gè)樣本中的表達(dá)量數(shù)據(jù)進(jìn)行假設(shè)檢驗(yàn)。常用的假設(shè)檢驗(yàn)方法有t-檢驗(yàn)、卡方檢驗(yàn)等。很多剛接觸轉(zhuǎn)錄組分析的人可能會(huì)有這樣一個(gè)疑問，一個(gè)轉(zhuǎn)錄本是不是差異表達(dá)，做完假設(shè)檢驗(yàn)看P-value不就可以了么？為什么會(huì)有FDR這樣一個(gè)新的概念出現(xiàn)？這是因?yàn)檗D(zhuǎn)錄組分析并不是針對(duì)一個(gè)或幾個(gè)轉(zhuǎn)錄本進(jìn)行分析，轉(zhuǎn)錄組分析的是一個(gè)樣品中所轉(zhuǎn)錄表達(dá)的所有轉(zhuǎn)錄本。所以，一個(gè)樣品當(dāng)中有多少轉(zhuǎn)錄本，就需要對(duì)多少轉(zhuǎn)錄本進(jìn)行假設(shè)檢驗(yàn)。這會(huì)導(dǎo)致一個(gè)很嚴(yán)重的問題，在單次假設(shè)檢驗(yàn)中較低的假陽性比例會(huì)累積到一個(gè)非常驚人的程度。舉個(gè)不太嚴(yán)謹(jǐn)?shù)睦印?/p> 假設(shè)現(xiàn)在有這樣一個(gè)項(xiàng)目： ● 包含兩個(gè)樣品，共得到10000條轉(zhuǎn)錄本的表達(dá)量數(shù)據(jù)， ● 其中有100條轉(zhuǎn)錄本的表達(dá)量在兩個(gè)樣品中是有差異的。 ● 針對(duì)單個(gè)基因的差異表達(dá)分析有1%的假陽性。 由于存在1%假陽性的結(jié)果，在我們分析完這10000個(gè)基因后，我們會(huì)得到100個(gè)假陽性導(dǎo)致的錯(cuò)誤結(jié)果，加上100條真實(shí)存在的結(jié)果，共計(jì)200個(gè)結(jié)果。在這個(gè)例子中，一次分析得到的200個(gè)差異表達(dá)基因中，有50%都是假陽性導(dǎo)致的錯(cuò)誤結(jié)果，這顯然是不可接受的。為了解決這個(gè)問題，F(xiàn)DR這個(gè)概念被引入，以控制最終得到的分析結(jié)果中假陽性的比例。 如何計(jì)算FDR FDR的計(jì)算是根據(jù)假設(shè)檢驗(yàn)的P-value進(jìn)行校正而得到的。一般來說，F(xiàn)DR的計(jì)算采用Benjamini-Hochberg方法（簡(jiǎn)稱BH法），計(jì)算方法如下： 將所有P-value升序排列.P-value記為P，P-value的序號(hào)記為i，P-value的總數(shù)記為m FDR(i)=P(i)*m/i 根據(jù)i的取值從大到小，依次執(zhí)行FDR(i)=min{FDR(i),FDR(i+1)}
注：實(shí)際上，BH法的原始算法是找到一個(gè)最大的i，滿足P≤i/m*FDR閾值，此時(shí)，所有小于i的數(shù)據(jù)就都可以認(rèn)為是顯著的。在實(shí)踐中，為了能夠在比較方便的用不同的FDR閾值對(duì)數(shù)據(jù)進(jìn)行分析，采用了步驟3里的方法。這個(gè)方法可以保證，不論FDR閾值選擇多少，都可以直接根據(jù)FDR的數(shù)值來直接找到所有顯著的數(shù)據(jù)。 下面我們以一個(gè)包含10個(gè)數(shù)據(jù)的例子來看一下FDR計(jì)算的過程差異表達(dá)分析之FDR 在這個(gè)例子中，第一列是原始的P-value，第二列是排序后的序號(hào)，第三列是根據(jù)P-value校正得到的初始FDR，第四列是最終用于篩選數(shù)據(jù)的FDR數(shù)值。如果我們?cè)O(shè)定FDR<0.05，那么綠色高亮的兩個(gè)數(shù)據(jù)就是最終分析認(rèn)為顯著的數(shù)據(jù)。 FDR的閾值選擇在轉(zhuǎn)錄組分析中是非常重要的一個(gè)環(huán)節(jié)，常用的閾值包括0.01、0.05、0.1等。實(shí)踐中也可以根據(jù)實(shí)際的需要來靈活選擇。例如，在做真菌或者原核生物的轉(zhuǎn)錄組分析時(shí)，由于這些物種轉(zhuǎn)錄本數(shù)量較少，假陽性累積的程度較低，所以可以適當(dāng)將FDR閾值設(shè)置的較高一些，這樣可以獲得較多的差異表達(dá)結(jié)果，有利于后續(xù)的分析。 淺談多重檢驗(yàn)校正FDR
Posted: 四月 12, 2017 Under: Basic By Kai no Comments 例如，在我們對(duì)鑒定到的差異蛋白做GO功能注釋后，通常會(huì)計(jì)算一個(gè)p值。當(dāng)某個(gè)蛋白的p值小于0.05（5%）時(shí)，我們通常認(rèn)為這個(gè)蛋白在兩個(gè)樣本中的表達(dá)是有差異的。但是仍舊有5%的概率，這個(gè)蛋白并不是差異蛋白。那么我們就錯(cuò)誤地否認(rèn)了原假設(shè)（在兩個(gè)樣本中沒有差異表達(dá)），導(dǎo)致了假陽性的產(chǎn)生（犯錯(cuò)的概率為5%）。 如果檢驗(yàn)一次，犯錯(cuò)的概率是5%；檢測(cè)10000次，犯錯(cuò)的次數(shù)就是500次，即額外多出了500次差異的結(jié)論（即使實(shí)際沒有差異）。 為了控制假陽性的次數(shù)，于是我們需要對(duì)p值進(jìn)行多重檢驗(yàn)校正，提高閾值。 第一種方法Bonferroni，最簡(jiǎn)單嚴(yán)厲的方法。
例如，如果檢驗(yàn)1000次，我們就講閾值設(shè)定為5% / 1000 = 0.00005；即使檢驗(yàn)1000次，犯錯(cuò)誤的概率還是保持在N×1000 = 5%。最終使得預(yù)期犯錯(cuò)誤的次數(shù)不到1次，抹殺了一切假陽性的概率。但是該方法雖然簡(jiǎn)單，但是檢驗(yàn)過于嚴(yán)格，導(dǎo)致最后找不到顯著表達(dá)的蛋白（假陰性）。 第二種方法FDR（False Discovery Rate）
相對(duì)Bonferroni來說，F(xiàn)DR用比較溫和的方法對(duì)p值進(jìn)行了校正。其試圖在假陽性和假陰性間達(dá)到平衡，將假/真陽性比例控制到一定范圍之內(nèi)。例如，如果檢驗(yàn)1000次，我們?cè)O(shè)定的閾值為0.05（5%），那么無論我們得到多少個(gè)差異蛋白，這些差異蛋白出現(xiàn)假陽性的概率保持在5%之內(nèi)，這就叫FDR＜5%。 那么我們?cè)趺磸膒 value 來估算FDR呢，人們?cè)O(shè)計(jì)了幾種不同的估算模型。其中使用最多的是Benjamini and Hochberg方法，簡(jiǎn)稱BH法。雖然這個(gè)估算公式并不夠完美，但是也能解決大部分的問題，主要還是簡(jiǎn)單好用！ FDR的計(jì)算方法
除了可以使用excel的BH計(jì)算方法外，對(duì)于較大的數(shù)據(jù)，我們推薦使用R命令p.adjust。 p.adjust(p, method = p.adjust.methods, n = length§) p.adjust.methods c(“holm”, “hochberg”, “hommel”, “bonferroni”, “BH”, “BY”,“fdr”, “none”)我們還可以從R命令p.adjust的源代碼，了解其運(yùn)行的機(jī)制是什么。 p.adjust
function (p, method = p.adjust.methods, n = length§){
method <- match.arg(method)
if (method == “fdr”)
method <- “BH”
nm <- names§
p <- as.numeric§
……
BH = {
i <- lp:1L
o <- order(p, decreasing = TRUE)
ro <- order(o)
pmin(1, cummin(n/i * p[o]))[ro]
}
……
p0
}
其實(shí)該函數(shù)表達(dá)的意思是這樣的：
我們將一系列p值、校正方法（BH）以及所有p值的個(gè)數(shù)（length§）輸入到p.adjust函數(shù)中。
將一系列的p值按照從大到小排序，然后利用下述公式計(jì)算每個(gè)p值所對(duì)應(yīng)的FDR值。 公式：p * (n/i)， p是這一次檢驗(yàn)的p value，n是檢驗(yàn)的次數(shù)，i是排序后的位置ID（如最大的P值的i值肯定為n，第二大則是n-1，依次至最小為1）。
將計(jì)算出來的FDR值賦予給排序后的p值，如果某一個(gè)p值所對(duì)應(yīng)的FDR值大于前一位p值（排序的前一位）所對(duì)應(yīng)的FDR值，則放棄公式計(jì)算出來的FDR值，選用與它前一位相同的值。因此會(huì)產(chǎn)生連續(xù)相同F(xiàn)DR值的現(xiàn)象；反之則保留計(jì)算的FDR值。
將FDR值按照最初始的p值的順序進(jìn)行重新排序，返回結(jié)果。 p值還是 FDR ?
差異分析
如何篩選顯著性差異基因，p value, FDR 如何選
經(jīng)常有同學(xué)詢問如何篩選差異的基因（蛋白）。已經(jīng)計(jì)算了表達(dá)量和p value值，差異的基因（蛋白）太多了，如何篩選。其中最為關(guān)鍵的是需要對(duì)p value進(jìn)行校正。 基本概念： 零假設(shè)：在隨機(jī)條件下的分布。 p值：在零假設(shè)下，觀測(cè)到某一特定實(shí)驗(yàn)結(jié)果的概率稱為p值。 假陽性：得到了陽性結(jié)果，但這個(gè)陽性結(jié)果是假的。 假陰性：得到了陰性結(jié)果，但這個(gè)陰性結(jié)果是假的。 單次檢驗(yàn)： 針對(duì)單個(gè)基因（蛋白），采用統(tǒng)計(jì)檢驗(yàn)，假設(shè)采用的p值為小于0.05，我們通常認(rèn)為這個(gè)基因在兩個(gè)（組）樣本中的表達(dá)是有顯著差異的，但是仍舊有5%的概率，這個(gè)基因并不是差異基因。 單多次檢驗(yàn)： 當(dāng)兩個(gè)（組）樣本中有10000個(gè)基因采用同樣的檢驗(yàn)方式進(jìn)行統(tǒng)計(jì)檢驗(yàn)時(shí)，這個(gè)時(shí)候就有一個(gè)問題，單次犯錯(cuò)的概率為0.05， 進(jìn)行10000次檢驗(yàn)的話，那么就有0.05*10000=500 個(gè)基因的差異被錯(cuò)誤估計(jì)了。 多重檢驗(yàn)矯正： 為了解決多次檢驗(yàn)帶來的問題，我們需要對(duì)多次檢驗(yàn)進(jìn)行校正。那如何校正呢？在此介紹兩種方法： Bonferroni 校正法
Bonferroni校正法：如果進(jìn)行N次檢驗(yàn)，那么p值的篩選的閾值設(shè)定為p/N。 比如，進(jìn)行10000次檢驗(yàn)的話，如果p值選擇為0.05， 那么校正的p值篩選為0.000005。 p值低于此的基因才是顯著性差異基因。
該方法雖然簡(jiǎn)單，但是過于嚴(yán)格，導(dǎo)致最后找的差異基因很少，甚至找不到差異的基因。 FDR(False Discovery Rate) 校正法
FDR錯(cuò)誤控制法是Benjamini于1995年提出的一種方法，基本原理是通過控制FDR值來決定p值的值域。相對(duì)Bonferroni來說，F(xiàn)DR用比較溫和的方法對(duì)p值進(jìn)行了校正。其試圖在假陽性和假陰性間達(dá)到平衡，將假/真陽性比例控制到一定范圍之內(nèi)。
那么怎么從p值來估算FDR呢，人們?cè)O(shè)計(jì)了幾種不同的估算模型。其中使用最多的是Benjamini and Hochberg方法，簡(jiǎn)稱BH法。該方法分兩步完成，具體如下：
2.1 假設(shè)總共有m個(gè)候選基因，每個(gè)基因?qū)?yīng)的p值從小到大排列分別是p(1),p(2),…,p(m)
2.2 若想控制FDR不能超過q，則只需找到最大的正整數(shù)i，使得 p(i)<= (i*q)/m . 然后，挑選對(duì)應(yīng)p(1),p(2),…,p(i)的基因做為差異表達(dá)基因，這樣就能從統(tǒng)計(jì)學(xué)上保證FDR不超過q。 如何實(shí)現(xiàn)多重檢驗(yàn)： 如果你了解R語言的話，那么采用p.adjust方法就可以了。 簡(jiǎn)單使用DESeq做差異分析
Posted: 五月 06, 2017 Under: Transcriptomics By Kai no Comments DESeq這個(gè)R包主要針對(duì)count data，其數(shù)據(jù)來源可以是RNA-Seq或者其他高通量測(cè)序數(shù)據(jù)。類似地，對(duì)于CHIP-Seq數(shù)據(jù)或者質(zhì)譜肽段數(shù)據(jù)也是使用的。 由于DESeq是一個(gè)R包，因此使用它需要一點(diǎn)點(diǎn)R基礎(chǔ)語法。 首先需要讀入一個(gè)數(shù)據(jù)框，列代表每個(gè)sample，行代表每個(gè)gene database_all <- read.table(file = “readcount”, sep = “\t”, header = T, row.names = 1)
database <- database_all[,1:6]
這里主要對(duì)于兩兩比較的數(shù)據(jù)，因此我取了數(shù)據(jù)的前6列，分別是兩組樣品，每組3個(gè)生物學(xué)重復(fù) 設(shè)定分組信息，也就是樣本分組的名稱 type <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
我這里是樣品1是LC_1，樣品2是LC_2 由于DESeq包要求接下來的count data必須要整數(shù)型，因此我們需要對(duì)數(shù)據(jù)進(jìn)行取整,然后將數(shù)據(jù)database和分組信息type讀入到cds對(duì)象中 database <- round(as.matrix(database))
cds <- newCountDataSet(database,type)
接下里對(duì)于不同類型的數(shù)據(jù)要進(jìn)行不同的處理，可以粗略分為有生物學(xué)重復(fù)數(shù)據(jù)、有部分生物學(xué)重復(fù)數(shù)據(jù)以及無生物學(xué)重復(fù)數(shù)據(jù) 4.1 有生物學(xué)重復(fù) cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
res <- nbinomTest(cds,“LC_1”,“LC_2”)
其通過estimate the dispersion并對(duì)count data進(jìn)行標(biāo)準(zhǔn)化，然后得到每個(gè)gene做T test檢驗(yàn) 4.2 對(duì)于部分樣品有生物學(xué)重復(fù) cds <- estimateSizeFactors(cds)
cds <- estimateDispersions(cds)
res <- nbinomTest(cds,“LC_1”,“LC_2”)
其步驟跟有上述的一樣的，DESeq會(huì)根據(jù)有生物學(xué)重復(fù)的樣品來estimate the dispersion，當(dāng)然要保證unreplicated condition does not have larger variation than the replicated one 4.3 對(duì)于沒有生物學(xué)重復(fù) cds <- estimateDispersions(cds, method=“blind”, sharingMode=“fit-only” )
res <- nbinomTest(cds,“LC_1”,“LC_2”)
注意參數(shù)method=”blind” 和 sharingMode=”fit-only”即可 最后就是查看符合閾值的差異基因有多少個(gè)即可，然后將結(jié)果輸出到csv文件中方便查看 table(res
p
a
d
j
<
0.05
)
r
e
s
<
?
r
e
s
[
o
r
d
e
r
(
r
e
s
padj <0.05) res <- res[order(res
padj<0.05)res<?res[order(respadj),]
sum(res$padj<=0.01,na.rm = T)
write.csv(resdata,file = “LC_1_vs_LC_2_DESeq.csv”)
Summary
DESeq在前幾年的文章中經(jīng)常被使用，但是現(xiàn)在有了其升級(jí)版DESeq2，后者相比前者對(duì)于犯第一類錯(cuò)誤卡的并不是那么嚴(yán)格了，所以在同樣的padj的閾值下，篩選到的差異基因的數(shù)目也會(huì)多一點(diǎn)。 并且上述中，我都只使用nbinomTest來獲得由顯著差異的gene，其實(shí)DESeq包還提供了其他檢驗(yàn)方法，比如nbinomGLMTest函數(shù)，具體可以查看DESeq文檔，還有其他對(duì)應(yīng)對(duì)因素分析的使用方法。 DESeq2和EdgeR都可用于做基因差異表達(dá)分析，主要也是用于RNA-Seq數(shù)據(jù)，同樣也可以處理類似的ChIP-Seq,shRNA以及質(zhì)譜數(shù)據(jù)。 這兩個(gè)都屬于R包，其相同點(diǎn)在于都是對(duì)count data數(shù)據(jù)進(jìn)行處理，都是基于負(fù)二項(xiàng)分布模型。因此會(huì)發(fā)現(xiàn)，用兩者處理同一組數(shù)據(jù)，最后在相同閾值下篩選出的大部分基因都是一樣的，但是有一部分不同應(yīng)該是由于其估計(jì)離散度的不同方法所導(dǎo)致的。 DESeq2的使用方法：
輸入矩陣數(shù)據(jù)，行名為sample，列名為gene；DESeq2不支持無生物學(xué)重復(fù)的數(shù)據(jù)，因此我選擇了2個(gè)樣本，3個(gè)生物學(xué)重復(fù)的數(shù)據(jù)；并對(duì)count data取整（經(jīng)大神指點(diǎn)，這里需要說明下，我的測(cè)試數(shù)據(jù)readcount是RSEM定量的結(jié)果，并不是常見的htseq-count的結(jié)果，所以count值會(huì)有小數(shù)點(diǎn)，而DESeq2包不支持count數(shù)有小數(shù)點(diǎn)，所以這里需要round取整）。 database_all <- read.table(file = “readcount”, sep = “\t”, header = T, row.names = 1)
database <- database_all[,1:6]
#type <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
database <- round(as.matrix(database))
設(shè)置分組信息以及構(gòu)建dds對(duì)象 condition <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
coldata <- data.frame(row.names = colnames(database), condition)
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
使用DESeq函數(shù)進(jìn)行估計(jì)離散度，然后進(jìn)行標(biāo)準(zhǔn)的差異表達(dá)分析，得到res對(duì)象結(jié)果 dds <- DESeq(dds)
res <- results(dds)
最后設(shè)定閾值，篩選差異基因，導(dǎo)出數(shù)據(jù) table(res
p
a
d
j
<
0.05
)
r
e
s
<
?
r
e
s
[
o
r
d
e
r
(
r
e
s
padj <0.05) res <- res[order(res
padj<0.05)res<?res[order(respadj),]
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by=“row.names”,sort=FALSE)
write.csv(resdata,file = “LC_1_vs_LC_2.csv”)
EdgeR的使用方法：
跟DESeq2一樣，EdgeR輸入矩陣數(shù)據(jù)，行名為sample，列名為gene；DESeq2不支持無生物學(xué)重復(fù)的數(shù)據(jù)，因此我選擇了2個(gè)樣本，3個(gè)生物學(xué)重復(fù)的數(shù)據(jù)。 exprSet_all <- read.table(file = “readcount”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
exprSet <- exprSet_all[,1:6]
group_list <- factor(c(rep(“LC_1”,3), rep(“LC_2”,3)))
設(shè)置分組信息，去除低表達(dá)量的gene以及做TMM標(biāo)準(zhǔn)化 exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,]
exprSet <- DGEList(counts = exprSet, group = group_list)
exprSet <- calcNormFactors(exprSet)
使用qCML（quantile-adjusted conditional maximum likelihood）估計(jì)離散度（只針對(duì)單因素實(shí)驗(yàn)設(shè)計(jì)） exprSet <- estimateCommonDisp(exprSet)
exprSet <- estimateTagwiseDisp(exprSet)
尋找差異gene(這里的exactTest函數(shù)還是基于qCML并且只針對(duì)單因素實(shí)驗(yàn)設(shè)計(jì))，然后按照閾值進(jìn)行篩選即可 et <- exactTest(exprSet)
tTag <- topTags(et, n=nrow(exprSet))
tTag <- as.data.frame(tTag)
write.csv(tTag,file = “LC_1_vs_LC_2_edgeR.csv”)
Summary
以上我主要針對(duì)單因素兩兩比較組進(jìn)行差異分析，其實(shí)DESeq2和EdgeR兩個(gè)R包都可以對(duì)多因素進(jìn)行差異分析。 DESeq2修改以上代碼的分組信息design參數(shù)以及在差異分析results函數(shù)中添加所選定的分組因素，其他代碼基本一樣，具體參照DESeq2手冊(cè) EdgeR則需要用Cox-Reid profile-adjusted likelihood (CR)方法來估算離散度，y <- estimateDisp(y, design)或者分別使用三個(gè)函數(shù)（y <- estimateGLMCommonDisp(y, design)，y <- estimateGLMTrendedDisp(y, design)， ）y <- estimateGLMTagwiseDisp(y, design)；然后差異表達(dá)分析也跟單因素分析不同，主要使用generalized linear model (GLM) likelihood ratio test 或者 quasi-likelihood(QL) F-test，具體代碼可以參照EdgeR手冊(cè)。 簡(jiǎn)單使用limma做差異分析
Posted: 五月 12, 2017 Under: Transcriptomics By Kai no Comments 首先需要說明的是，limma是一個(gè)非常全面的用于分析芯片以及RNA-Seq的差異分析，按照其文章所說： limma is an R/Bioconductor software package that provides an integrated solution for analysing data from gene expression experiments. 在這我只是對(duì)其中的一種情況進(jìn)行簡(jiǎn)單的總結(jié)，比如這個(gè)包可以處理RNA-Seq數(shù)據(jù)，我簡(jiǎn)單的以兩個(gè)比較組進(jìn)行分組為例，至于其他分組情況，請(qǐng)看limma說明文檔，有非常詳細(xì)的說明，非常親民。 首先我們還是輸入count矩陣，這里也跟其他差異分析R包一樣，不要輸入已經(jīng)標(biāo)準(zhǔn)化的數(shù)據(jù)。順便也加載下edgeR這個(gè)R包 library(limma)
library(edgeR)
counts <- read.table(file = “conut_all.txt”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
接著按照文檔的說明以及l(fā)imma包的習(xí)慣，我們需要對(duì)count進(jìn)行標(biāo)準(zhǔn)化以及轉(zhuǎn)化為log2的值，這里標(biāo)準(zhǔn)化的方法為TMM，使用edgeR里面的calcNormFactors函數(shù)即可 dge <- DGEList(counts = counts)
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)
這里prior.count值我粗略理解為是為了防止log2()遇到過于小的值 然后跟其他包一樣，設(shè)置分組信息 group_list <- factor(c(rep(“control”,2), rep(“siSUZ12”,2)))
design <- model.matrix(~group_list)
colnames(design) <- levels(group_list)
rownames(design) <- colnames(counts)
接下來就是常規(guī)的差異分析 fit <- lmFit(logCPM, design)
fit <- eBayes(fit, trend=TRUE)
output <- topTable(fit, coef=2,n=Inf)
sum(output$adj.P.Val<0.05)
到這里為止，我們主要是用了limma包里對(duì)RNA-Seq差異分析的limma-trend方法，該方法主要適用于樣本間測(cè)序深度基本保持一致的情況，也就是說兩個(gè)樣本的文庫(reads數(shù)目)大小相差的不懸殊（說明文檔中是默認(rèn)3倍以內(nèi)？） 當(dāng)文庫大小在樣本間變化幅度相當(dāng)大的話，可以使用limma的voom方法，可按照下面的代碼流程： count數(shù)據(jù)的輸入以及數(shù)據(jù)標(biāo)準(zhǔn)化還是跟之前的一樣 counts <- read.table(file = “conut_all.txt”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
dge <- DGEList(counts = counts)
dge <- calcNormFactors(dge)
還是一樣需要分組信息 group_list <- factor(c(rep(“control”,2), rep(“siSUZ12”,2)))
design <- model.matrix(~group_list)
colnames(design) <- levels(group_list)
rownames(design) <- colnames(counts)
接下來進(jìn)行voom轉(zhuǎn)化 v <- voom(dge, design, plot=TRUE)
其實(shí)可以不進(jìn)行TMM標(biāo)準(zhǔn)化，直接用count數(shù)據(jù)進(jìn)行voom轉(zhuǎn)化，如： v <- voom(counts, design, plot=TRUE)
最后就是普通的差異分析過程 fit <- lmFit(v, design)
fit <- eBayes(fit)
output <- topTable(fit, coef=2,n=Inf)
sum(output$adj.P.Val<0.05)
Summary
Limma長(zhǎng)久以來就是一個(gè)非常流行的差異分析R包，其內(nèi)容涉及的非常廣泛，用于RNA-Seq只是其內(nèi)容的一小部分，并且使其處理RNA-Seq數(shù)據(jù)也使用芯片類似線性模型下，并且按照其說法，limma包比其他基于負(fù)二項(xiàng)式分布模型的差異分析R包更加的優(yōu)秀。 其實(shí)差異分析不外乎數(shù)據(jù)的標(biāo)準(zhǔn)化以及統(tǒng)計(jì)模型分析差異兩個(gè)方面的作用，每個(gè)差異分析R包都有其自身的優(yōu)點(diǎn)，個(gè)人理解，取舍在于自己的理解以及想法即可。 其實(shí)，自己對(duì)于limma包的理解還是比較粗淺的 轉(zhuǎn)錄組差異表達(dá)分析小實(shí)戰(zhàn)（一）
Posted: 七月 28, 2017 Under: Transcriptomics By Kai no Comments 讀文獻(xiàn)獲取數(shù)據(jù)
文獻(xiàn)名稱：AKAP95 regulates splicing through scaffolding
RNAs and RNA processing factors 查找數(shù)據(jù)：Data availability
The RIP-seq an RNA-seq data have been deposited in the Gene
Expression Omnibus database, with accession code GSE81916. All other data is
available from the author upon reasonable request. 獲得GSE號(hào)：GSE81916 下載測(cè)序數(shù)據(jù)
https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE81916獲取數(shù)據(jù)信息，并點(diǎn)擊網(wǎng)址下方的ftp，下載測(cè)序數(shù)據(jù) 從https://trace.ncbi.nlm./Traces/study/?acc=PRJNA323422可知我們需要的mRNA測(cè)序編號(hào)為SRR3589956到SRR3589962 通過Apera下載SRR數(shù)據(jù)，這里以SRR3589956為例： ascp -T -i /home/anlan/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp-private.ncbi.nlm.:sra/sra-instant/reads/ByRun/sra/SRR/SRR358/SRR3589956/SRR3589956.sra ./
轉(zhuǎn)化fastq測(cè)序數(shù)據(jù)
通過sratoolkit工具將SRR文件轉(zhuǎn)化為fastq格式的測(cè)序數(shù)據(jù)（寫了個(gè)shell循環(huán)） for i in
(
s
e
q
5662
)
;
d
o
n
o
h
u
p
f
a
s
t
q
?
d
u
m
p
?
?
s
p
l
i
t
?
3
S
R
R
35899
(seq 56 62);do nohup fastq-dump --split-3 SRR35899
(seq5662);donohupfastq?dump??split?3SRR35899{i} &;done
通過fastqc對(duì)每個(gè)fastq文件進(jìn)行質(zhì)檢，用multiqc查看整體質(zhì)檢報(bào)告（對(duì)當(dāng)前目錄下的fastq測(cè)序結(jié)果進(jìn)行質(zhì)檢，生成每個(gè)fq文件的質(zhì)檢報(bào)告總multiqc整合后統(tǒng)計(jì)查看） fastqc *.fastq
multiqc ./
點(diǎn)擊這個(gè)url可以查看我這個(gè)multiqc報(bào)告：http://www./data/multiqc_report.html 如果有接頭或者質(zhì)量值不達(dá)標(biāo)的需要進(jìn)行過濾，這次的數(shù)據(jù)質(zhì)量都不錯(cuò)，因此直接進(jìn)行比對(duì)即可
序列比對(duì)
安裝hisat2軟件，下載人類的hiast2索引文件 hisat2下載并安裝： ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/downloads/hisat2-2.1.0-Linux_x86_64.zip
unzip hisat2-2.1.0-Linux_x86_64.zip
下載hisat2的human索引 ftp://ftp.ccb.jhu.edu/pub/infphilo/hisat2/data/hg19.tar.gz
tar zxvf hg19.tar.gz
用hisat2進(jìn)行比對(duì)，測(cè)序數(shù)據(jù)放在data目錄下，索引文件放在reference/index/hisat2/hg19目錄下，SRR3589956-SRR3589958為人的測(cè)序數(shù)據(jù) for i in KaTeX parse error: Undefined control sequence: \ at position 28: …do hisat2 -p 4 \? ?-x ~/reference/…{i}_1.fastq -2 ./data/SRR35899KaTeX parse error: Undefined control sequence: \ at position 13: {i}_2.fastq \? ?-S SRR35899i.sam >SRR35899${i}.log;done
用samtools將sam文件轉(zhuǎn)化為bam文件，并使用默認(rèn)排序 for i in
(
s
e
q
5658
)
;
d
o
s
a
m
t
o
o
l
s
s
o
r
t
?
@
5
?
o
S
R
R
35899
(seq 56 58);do samtools sort -@ 5 -o SRR35899
(seq5658);dosamtoolssort?@5?oSRR35899{i}.bam SRR35899${i}.sam;done
reads計(jì)數(shù)
用htseq對(duì)比對(duì)產(chǎn)生的bam進(jìn)行count計(jì)數(shù) htseq安裝，使用miniconda，省事！唯一的問題是htseq版本不是最新的，是0.7.2。想要最新版還是要正常安裝，可參考http://www./thread-1847-1-2.html conda install -c bioconda htseq
用htseq將對(duì)比后的結(jié)果進(jìn)行計(jì)數(shù) for i in KaTeX parse error: Undefined control sequence: \ at position 48: …m -r pos -s no \? ?SRR35899{i}.bam ~/reference/genome/hg19/gencode.v26lift37.annotation.gtf
1>SRR35899
i
.
c
o
u
n
t
2
>
S
R
R
35899
{i}.count 2>SRR35899
i.count2>SRR35899{i}_htseq.log;done
將3個(gè)count文件（SRR3589956.count，SRR3589957.count，SRR3589958.count）合并成一個(gè)count矩陣，這是就需要腳本來解決這個(gè)問題，不然其他方法會(huì)稍微麻煩點(diǎn) #!/usr/bin/perl -w
use strict; my $path = shift @ARGV;
opendir DIR, $path or die;
my @dir = readdir DIR; my $header;
my @sample;
my %hash;
foreach my KaTeX parse error: Expected '}', got 'EOF' at end of input: …dir) { if (file =~ /^\w+.*.count/) {
push @sample, $file;
KaTeX parse error: Undefined control sequence: \t at position 12: header .= "\?t?file";
open my $fh,
f
i
l
e
o
r
d
i
e
;
w
h
i
l
e
(
<
file or die; while (<
fileordie;while(<fh>) {
chomp;
next if ($_ =~ /^\W+/);
my @array = split /\t/, $_;
KaTeX parse error: Expected '}', got 'EOF' at end of input: hash{array[0]} -> {$file} = $array[1];
}
close KaTeX parse error: Expected 'EOF', got '}' at position 9: fh; }? } print "header\n";
map{
my $gene =
;
p
r
i
n
t
"
_; print "
;print"gene";
foreach my KaTeX parse error: Undefined control sequence: \t at position 33: … print "\?t?".hash{KaTeX parse error: Expected 'EOF', got '}' at position 5: gene}? -> {file};
}
print “\n”;
}keys %hash;
按照接下來的劇本，應(yīng)該講count_matrix文件導(dǎo)入DESeq進(jìn)行差異表達(dá)分析。但是從這篇文章的Bioinformatic analyses部分可以發(fā)現(xiàn)，作者的control組的2組數(shù)據(jù)是來自2個(gè)不同的批次（一個(gè)是SRR3589956，另外一個(gè)來源GSM1095127 in GSE44976），treat組倒是同一個(gè)批次（SRR3589957和SRR3589958）。但是對(duì)于Mouse cells來說，倒是滿足2個(gè)control和2個(gè)treat都正常來自同個(gè)批次，因此打算重新用SRR3589959-SRR3589962重新做個(gè)一個(gè)count_matrix進(jìn)行后續(xù)差異分析 轉(zhuǎn)錄組差異表達(dá)分析小實(shí)戰(zhàn)（二）
Posted: 八月 14, 2017 Under: Transcriptomics By Kai no Comments 差異基因表達(dá)分析
我按照前面的流程轉(zhuǎn)錄組差異表達(dá)分析小實(shí)戰(zhàn)（一），將小鼠的4個(gè)樣本又重新跑了一遍，從而獲得一個(gè)新的count文件：mouse_all_count.txt，有需要的話，可以下載下來進(jìn)行后續(xù)的差異分析。 一般來說，由于普遍認(rèn)為高通量的read count符合泊松分布，所以一些差異分析的R包都是基于負(fù)二項(xiàng)式分布模型的，比如DESeq、DESeq2和edgeR等，所以這3個(gè)R包從整體上來說是類似的（但各自標(biāo)準(zhǔn)化算法是不一樣的）。 當(dāng)然還有一個(gè)常用的R包則是Limma包，其中的limma-trend和limma-voom都能用來處理RNA-Seq數(shù)據(jù)（但對(duì)應(yīng)適用的情況不一樣） 下面準(zhǔn)備適用DESeq2和edgeR兩個(gè)R包分別對(duì)小鼠的count數(shù)據(jù)進(jìn)行差異表達(dá)分析，然后取兩者的結(jié)果的交集的基因作為差異表達(dá)基因。 DEseq2 library(DESeq2)
##數(shù)據(jù)預(yù)處理
database <- read.table(file = “mouse_all_count.txt”, sep = “\t”, header = T, row.names = 1)
database <- round(as.matrix(database)) ##設(shè)置分組信息并構(gòu)建dds對(duì)象
condition <- factor(c(rep(“control”,2),rep(“Akap95”,2)), levels = c(“control”, “Akap95”))
coldata <- data.frame(row.names = colnames(database), condition)
dds <- DESeqDataSetFromMatrix(countData=database, colData=coldata, design=~condition)
dds <- dds[ rowSums(counts(dds)) > 1, ] ##使用DESeq函數(shù)估計(jì)離散度，然后差異分析獲得res對(duì)象
dds <- DESeq(dds)
res <- results(dds) #最后設(shè)定閾值，篩選差異基因，導(dǎo)出數(shù)據(jù)(全部數(shù)據(jù)。包括標(biāo)準(zhǔn)化后的count數(shù))
res <- res[order(res$padj),]
diff_gene <- subset(res, padj < 0.05 & (log2FoldChange > 1 | log2FoldChange < -1))
diff_gene_DESeq2 <- row.names(diff_gene)
resdata <- merge(as.data.frame(res), as.data.frame(counts(dds, normalized=TRUE)),by=“row.names”,sort=FALSE)
write.csv(resdata,file = “control_vs_Akap95.csv”,row.names = F)
最終獲得572個(gè)差異基因（篩選標(biāo)準(zhǔn)為padj < 0.05, |log2FoldChange| > 1） edgeR library(edgeR)
##跟DESeq2一樣，導(dǎo)入數(shù)據(jù)，預(yù)處理（用了cpm函數(shù)）
exprSet<- read.table(file = “mouse_all_count.txt”, sep = “\t”, header = TRUE, row.names = 1, stringsAsFactors = FALSE)
group_list <- factor(c(rep(“control”,2),rep(“Akap95”,2)))
exprSet <- exprSet[rowSums(cpm(exprSet) > 1) >= 2,] ##設(shè)置分組信息，并做TMM標(biāo)準(zhǔn)化
exprSet <- DGEList(counts = exprSet, group = group_list)
exprSet <- calcNormFactors(exprSet) ##使用qCML（quantile-adjusted conditional maximum likelihood）估計(jì)離散度（只針對(duì)單因素實(shí)驗(yàn)設(shè)計(jì)）
exprSet <- estimateCommonDisp(exprSet)
exprSet <- estimateTagwiseDisp(exprSet) ##尋找差異gene(這里的exactTest函數(shù)還是基于qCML并且只針對(duì)單因素實(shí)驗(yàn)設(shè)計(jì))，然后按照閾值進(jìn)行篩選即可
et <- exactTest(exprSet)
tTag <- topTags(et, n=nrow(exprSet))
diff_gene_edgeR <- subset(tTagKaTeX parse error: Expected 'EOF', got '&' at position 19: …le, FDR < 0.05 &? (logFC > 1 | l…table,file = “control_vs_Akap95_edgeR.csv”)
最終獲得688個(gè)差異基因（篩選標(biāo)準(zhǔn)為FDR < 0.05, |log2FC| > 1） 取DESeq2和edgeR兩者結(jié)果的交集 diff_gene <- diff_gene_DESeq2[diff_gene_DESeq2 %in% diff_gene_edgeR]
最終的差異基因數(shù)目為545個(gè) head(diff_gene)
[1] “ENSMUSG00000003309.14” “ENSMUSG00000046323.8” “ENSMUSG00000001123.15”
[4] “ENSMUSG00000023906.2” “ENSMUSG00000044279.15” “ENSMUSG00000018569.12”
其他兩個(gè)R包（DESeq和limma）就不在這嘗試了，我之前做過對(duì)于這4個(gè)R包的簡(jiǎn)單使用筆記，可以參考下：
簡(jiǎn)單使用DESeq做差異分析
簡(jiǎn)單使用DESeq2/EdgeR做差異分析
簡(jiǎn)單使用limma做差異分析 GO&&KEGG富集分析
以前一直沒有機(jī)會(huì)用Y叔寫的clusterProfiler包，這次正好看說明用一下。 GO富集，加載clusterProfiler包和org.Mm.eg.db包（小鼠嘛），然后將ENSEMBL ID后面的版本號(hào)去掉，不然后面不識(shí)別這個(gè)ID，然后按照clusterProfiler包的教程說明使用函數(shù)即可。 library(clusterProfiler)
library(org.Mm.eg.db) ##去除ID的版本號(hào)
diff_gene_ENSEMBL <- unlist(lapply(diff_gene, function(x){strsplit(x, “\.”)[[1]][1]}))
##GOid mapping + GO富集
ego <- enrichGO(gene = diff_gene_ENSEMBL, OrgDb = org.Mm.eg.db,
keytype = “ENSEMBL”, ont = “BP”, pAdjustMethod = “BH”,
pvalueCutoff = 0.01, qvalueCutoff = 0.05)
##查看富集結(jié)果數(shù)據(jù)
enrich_go <- as.data.frame(ego)
##作圖
barplot(ego, showCategory=10)
dotplot(ego)
enrichMap(ego)
plotGOgraph(ego)
KEGG富集，首先需要將ENSEMBL ID轉(zhuǎn)化為ENTREZ ID，然后使用ENTREZ ID作為kegg id，從而通過enrichKEGG函數(shù)從online KEGG上抓取信息，并做富集 library(clusterProfiler)
library(org.Mm.eg.db) ##ID轉(zhuǎn)化
ids <- bitr(diff_gene_ENSEMBL, fromType = “ENSEMBL”, toType = “ENTREZID”, OrgDb = “org.Mm.eg.db”)
kk <- enrichKEGG(gene = ids[,2], organism = “mmu”, keyType = “kegg”,
pvalueCutoff = 0.05, pAdjustMethod = “BH”, qvalueCutoff = 0.1)
##查看富集結(jié)果數(shù)據(jù)
enrich_kegg <- as.data.frame(kk)
##作圖
dotplot(kk)
到這里為止，轉(zhuǎn)錄組的差異表達(dá)分析算是做完了，簡(jiǎn)單的來說，這個(gè)過程就是將reads mapping 到reference上，然后計(jì)數(shù)獲得count數(shù)，然后做差異分析，最后來個(gè)GO KEGG，over了。。。 對(duì)于mapping和計(jì)數(shù)這兩部還有其實(shí)還有好多軟件，具體可見文獻(xiàn)：Gaining comprehensive biological insight into the transcriptome by performing a broad-spectrum RNA-seq analysis，有時(shí)間可以都嘗試下。 至于GO && KEGG這兩步，對(duì)于人、小鼠等模式物種來說，不考慮方便因素來說，完全可以自己寫腳本來完成，數(shù)據(jù)可以從gene ontology官網(wǎng)下載，然后就是GO id與gene id相互轉(zhuǎn)化了。KEGG 也是一樣，也可以用腳本去KEGG網(wǎng)站上自行抓取，先知道URL規(guī)則，然后爬數(shù)據(jù)即可。
|
