寫在前面：新手老菜鳥，個(gè)人心得伯复，不足之處歡迎交流指正鼠次；
寫作不易更哄，轉(zhuǎn)載前請(qǐng)聯(lián)系我征得同意芋齿，謝謝。
1.基本原理
1.1 概述
- 全稱：DESeq2 package for differential analysis of count data;
- 利用負(fù)二項(xiàng)分布廣義線性模型( negative binomial generalized linear models）成翩，同時(shí)觅捆，還利用了離散型估計(jì)、logFoldChange;
- 負(fù)二項(xiàng)分布是一個(gè)離散分布麻敌，符合測(cè)序reads分布栅炒；
1.2 構(gòu)建dds
- 要求輸入原始 reads count 數(shù)；不接受已經(jīng)做過(guò)處理的FPKM/TPM等术羔，因?yàn)檐浖凶约旱臉?biāo)準(zhǔn)化計(jì)算方法赢赊；
- 構(gòu)建dds。需要設(shè)置design公式级历，即告訴軟件你的數(shù)據(jù)是怎樣來(lái)的释移，基本試驗(yàn)設(shè)計(jì)如何，軟件會(huì)根據(jù)幾個(gè)變量綜合計(jì)算寥殖；
一般：design =~ variable1 + variable2 + ...玩讳；
只有一個(gè)變量時(shí)：design=~ condition；
很多醫(yī)學(xué)分析會(huì)加入年齡嚼贡、性別等：design=~sex+disease+condition熏纯；
可以對(duì)應(yīng)幾個(gè)變量，但如果沒(méi)有額外參數(shù)编曼，log2FC和p值是默認(rèn)對(duì)design公式中的最后一個(gè)變量或者最后一個(gè)因子與參考因子進(jìn)行比較豆巨；
1.3 函數(shù)與計(jì)算
1.3.1 標(biāo)準(zhǔn)化：DESeq函數(shù)
- 不同樣品的測(cè)序量有差異剩辟，最簡(jiǎn)單的標(biāo)準(zhǔn)化方式是計(jì)算counts per million
(CPM) = 原始reads count ÷ 總reads數(shù) x 1,000,000掐场；- 這種計(jì)算方式，易受到極高表達(dá)且在不同樣品中存在差異表達(dá)的基因的影響：這些基因的打開或關(guān)閉會(huì)影響到細(xì)胞中總的分子數(shù)目贩猎，可能導(dǎo)致這些基因標(biāo)準(zhǔn)化之后就不存在表達(dá)差異了熊户，而原本沒(méi)有差異的基因標(biāo)準(zhǔn)化之后卻有差異了；
- RPKM吭服、FPKM嚷堡、TPM 是 CPM按照基因或轉(zhuǎn)錄本長(zhǎng)度歸一化后的表達(dá)，都會(huì)受到這一影響艇棕；
- DESeq2的方法：
- 量化因子 (size factor,SF)蝌戒，首先計(jì)算每個(gè)基因在所有樣品中表達(dá)的幾何平均值；每個(gè)細(xì)胞的SF是所有基因與其在所有樣品中的表達(dá)值的幾何平均值的比值的中位數(shù)沼琉；由于幾何平均值的使用北苟，只有在所有樣品中表達(dá)都不為0的基因才能用來(lái)計(jì)算。這一方法又被稱為RLE(relative log expression)打瘪。
- 不但考慮了測(cè)序深度的問(wèn)題友鼻，還考慮了表達(dá)量超高或者極顯著差異表達(dá)的基因?qū)е耤ount的分布出現(xiàn)偏倚傻昙。
- DESeq函數(shù)分析：
- 三步：estimation of size factors（estimateSizeFactors)， estimation of dispersion（estimateDispersons)彩扔， Negative Binomial GLM fitting and Wald statistics（nbinomWaldTest）妆档；
- 可以分步運(yùn)行，也可一步到位虫碉，最后返回 results可用的DESeqDataSet對(duì)象贾惦。
1.3.2 歸一化：rlog/vst
是我自己取的名字，可能不準(zhǔn)確敦捧，我用在對(duì)dds進(jìn)行vst然后做PCA分析纤虽。
全稱：快速估算離散趨勢(shì)并應(yīng)用方差穩(wěn)定轉(zhuǎn)換陷谱；
若 samples<30 用rlog函數(shù)饲做，>30用vst；
類似的函數(shù)：gmodels - fast.prcomp蜘渣，輸入數(shù)據(jù)為TPM济蝉；或者TMM杰刽；
1.3.3 數(shù)據(jù)收縮：lfcShrink：
shrink the log2 foldchange，不會(huì)改變顯著差異的基因總數(shù)王滤，作者很推薦這個(gè)新功能贺嫂。
- 為何采用lfcShrink？log2FC estimates do not account for the large dispersion we observe with low read counts. 因此雁乡，兩種數(shù)據(jù)特別需要：低表達(dá)量占比高的第喳；數(shù)據(jù)特別分散的。
- 說(shuō)的就是我的數(shù)據(jù)磅馍浴曲饱！但是我只用來(lái)做MA plot并沒(méi)用來(lái)差異分析，因?yàn)椋?/li>
- lfcShrink 不改變p值q值珠月，但改變了fc扩淀，使 foldchange范圍變小，所以選擇DEG時(shí)會(huì)有不同結(jié)果啤挎，一般會(huì)偏少驻谆！所以，根據(jù)數(shù)據(jù)情況庆聘，本次分析DEG還是不做shrink胜臊。
1.3.4 p-value和q-value
- 作者給出的建議：
Need to filter on adjusted p-values, not p-values, to obtain FDR control. 10% FDR is common because RNA-seq experiments are often exploratory and having 90% true positives in the gene set is ok.
即：用padj為標(biāo)準(zhǔn)做結(jié)果篩選。- 事實(shí)上伙判，在軟件計(jì)算過(guò)程中象对，多次以
alpha表示padj，并默認(rèn)alpha=0.1澳腹；
1.3.5 MA plot
- MA plot也叫 mean-difference plot或者Bland-Altman plot织盼，用來(lái)估計(jì)模型中系數(shù)的分布;
- X軸, the "A"(average)杨何；Y軸，the "M"(minus) – subtraction of log values is equivalent to the log of the ratio;
- M表示log fold change沥邻，衡量基因表達(dá)量變化危虱，上調(diào)還是下調(diào)；A表示每個(gè)基因的count的均值唐全；
- 根據(jù)summary(res)可知埃跷，low count的比率很高，所以大部分基因表達(dá)量不高邮利，也就是集中在0的附近（log2(1)=0弥雹，也就是變化1倍），提供了模型預(yù)測(cè)系數(shù)的分布總覽延届。
本次做了lfcshrink+batch之后剪勿，MA圖趨于正常，比較集中在0附近方庭，一些差異基因也可以明顯看出厕吉。
1.4 DESeq(dds)結(jié)果矩陣每一列的含義：
- baseMean： is a just the average of the normalized count values, dividing by size factors, taken over all samples in the DESeqDataSet；是對(duì)照組的樣本標(biāo)準(zhǔn)化counts的均值械念；
- log2FoldChange： the effect size estimate. It tells us how much the gene’s expression seems to have changed due to treatment with dexamethasone in comparison to untreated samples头朱；也不是簡(jiǎn)單的用標(biāo)準(zhǔn)化的counts進(jìn)行計(jì)算，因?yàn)橛?jì)算的時(shí)候需要考慮零值以及其他效應(yīng)龄减；結(jié)果是log2fc(trt/untrt)所以要注意對(duì)照和處理的指定项钮；
- lfcSE： the standard error estimate for the log2 fold change estimate，(the effect size estimate has an uncertainty associated with it,)希停；
- p value: statistical test , the result of this test is reported as a p value. Remember that a p value indicates the probability that a fold change as strong as the observed one, or even stronger, would be seen under the situation described by the null hypothesis烁巫；
- p value有時(shí)候是NA：Sometimes a subset of the p values in res will be NA (not available); This is DESeq's way of reporting that all counts for this gene were zero, and hence no test was applied. In addition, p values can be assigned NA if the gene was excluded from analysis because it contained an extreme count outlier. For more information, see the outlier detection section of the DESeq2 vignette.
- padj: adjusted p-value;
1.5 實(shí)際遇到的其它問(wèn)題
1.5.1 pre-analysis 預(yù)分析
就是開始熟悉你的數(shù)據(jù)，知道ta的特點(diǎn)脖苏，ta是什么樣子什么脾氣秉性程拭，什么方式更能發(fā)現(xiàn)真實(shí)的ta，什么方式能揚(yáng)長(zhǎng)避短棍潘！選擇合適的分析方法；
先做PCA崖媚！
方法包括且不僅包括：
gmodels - fast.prcomp
DESeq2 - vst - plotPCA
1.5.2 批次效應(yīng)
1. 本項(xiàng)目的批次效應(yīng)：
- 根據(jù)預(yù)先進(jìn)行的PCA和correlation結(jié)果亦歉，本項(xiàng)目樣品individual間的差異性大于不同處理的差異；
- 若不去除畅哑，則根據(jù)不同A處理尋找DEG時(shí)必然會(huì)因?yàn)閕ndividual的影響而覆蓋掉一些肴楷，導(dǎo)致結(jié)果偏少；（test事實(shí)證明確是如此）荠呐；（本項(xiàng)目也做了removebatchEffect但是僅僅為了展示和驗(yàn)證）赛蔫；
- 因此砂客，將4個(gè)體當(dāng)做4個(gè)批次，寫入
design公式：
design =~ individual + condition
- 一般批次效應(yīng)：
- 可以用
limma removeBatchEffect或者Combat等去除呵恢；- 但是在做差異分析時(shí)鞠值，
ballgown, DESeq2等軟件建議不要提前去批次，而是將批次作為covariate進(jìn)行分析渗钉。
1.5.3 多組比較
若有好幾組處理彤恶，想兩兩比較，是分開準(zhǔn)備數(shù)據(jù)鳄橘、DESeq声离、results？還是全部數(shù)據(jù)一起瘫怜？（注意：是幾組之間兩兩比較术徊，不是一比多，一比多請(qǐng)移步多重t檢驗(yàn)或者WGCNA鲸湃，作者M(jìn)ike Love也說(shuō)不能做；」亍）
官方建議：
- 如有多個(gè)組需要比較，建議不要將其兩兩分開而是一起分析唤锉，通過(guò)在results時(shí)指定contrast對(duì)象世囊，獲得兩兩的比較結(jié)果，這樣可以綜合考慮所以樣品中的表達(dá)計(jì)算量化因子做DESeq窿祥；
- 但是株憾，如果你通過(guò)PCA/EDA分析發(fā)現(xiàn)某一組或某幾組的within-group variability比其他組的大，那么還是還是兩兩分組分開比較吧晒衩！
舉例嗤瞎，下面的情況就不適合一起做：（這可不就是我的數(shù)據(jù)嘛！所以听系，分開比較分別做DESeq更sensitive贝奇。）
PCA
1.5.4 低表達(dá)量基因過(guò)濾
- 實(shí)際分析完后，我發(fā)現(xiàn)一些假陽(yáng)性DEG靠胜，即由于一個(gè)樣本中出現(xiàn)極高reads而其它樣本都是0 reads導(dǎo)致的DEG掉瞳，在bioconductor上發(fā)帖得到了作者回復(fù)，于是加入附加條件：篩除在少于3個(gè)樣品中表達(dá)量少于10 reads的基因浪漠，再做DESeq標(biāo)準(zhǔn)化和DEG篩選陕习。**
- 舉例：有一個(gè)基因僅僅因?yàn)樵谔幚?中有一個(gè)表達(dá)而其他都是0，就被認(rèn)為是處理2/1的上調(diào)和處理3/2的下調(diào)址愿，很不靠譜该镣；
- 解決：（這個(gè)標(biāo)準(zhǔn)可以根據(jù)數(shù)據(jù)特點(diǎn)，自己制定）
keep <- rowSums(counts(dds) >= 10) >= 3
dds <- dds[keep, ]
1.5.5 outliers離群值處理
這部分提前寫了响谓，outliers是在results之后损合，summary(res)可以看到差異比較的一個(gè)基本結(jié)果省艳，有一項(xiàng)outliers，若占比很高數(shù)量成百上千嫁审，則要引起警惕跋炕。
解決方法：
- 軟件作者提到：
The automatic outlier filtering/replacement is most useful in situations which the number of outliers is limited. When there are thousands of reported outliers, it might make more sense to turn off the outlier filtering/replacement (DESeq with minReplicatesForReplace=Inf and results with cooksCutoff=FALSE) and perform manual inspection.
就是說(shuō)離群值太多的話還是關(guān)閉這一篩選功能，方法也提到了：
minReplicatesForReplace = Inf
results with cooksCutoff = FALSE
- 實(shí)踐證明兩個(gè)條件擇一即可土居，outliers消失枣购，很多假陰性降低。
- 分析前可以做一個(gè)歐氏距離計(jì)算擦耀，看是不是有樣本特別高或者低棉圈，可以剔除，方法見(jiàn)腳本眷蜓。
2. 實(shí)操
經(jīng)過(guò)多次分析和調(diào)整分瘾，最后用的代碼是：
（1）安裝DESeq2包
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("DESeq2")
（2）導(dǎo)入數(shù)據(jù)兩兩比較
setwd('xxx')
colData <- read.table('group.txt', header=TRUE, row.names=1)
readscount <- read.table('readscount.txt', header=TRUE, row.names=1)
condition <- factor(c(rep("treat1",10),rep("treat2",10)))
individual <- factor(c(rep("idv1",3),"idv2",rep("idv3",3),rep("idv4",3),
rep("idv1",3),rep("idv2",3),rep("idv4",4)))
colData
head(readscount)
condition
individual
library(DESeq2)
dds <- DESeqDataSetFromMatrix(readscount, colData, design =~ individual + condition)
keep <- rowSums(counts(dds) >= 10) >= 3 #過(guò)濾低表達(dá)基因，至少在3個(gè)樣品中都有10個(gè)reads
dds <- dds[keep, ]
（3）PCA（還有全部樣品的PCA在另一個(gè)腳本）
vsdata <- vst(dds, blind=FALSE) #歸一化
assay(vsdata) <- limma::removeBatchEffect(assay(vsdata), vsdata$individual) #去批次效應(yīng)
plotPCA(vsdata, intgroup = "condition") #自帶函數(shù)
（4）差異分析
dds_norm <- DESeq(dds, minReplicatesForReplace = Inf) #標(biāo)準(zhǔn)化; 不剔除outliers; 與cookscutoff結(jié)果相同
dds_norm$condition #保證是levels是按照后一個(gè)比前一個(gè)即trt/untrt吁系，否則需在results時(shí)指定
res <- results(dds_norm, contrast = c("condition","treat2","treat1"), cooksCutoff = FALSE) #alpha=0.05可指定padj; cookCutoff是不篩選outliers因?yàn)樘嗔?summary(res)
#resOrdered <- res[order(res$pvalue), ] #排序
sum(res$padj<0.05, na.rm = TRUE)
res_data <- merge(as.data.frame(res),
as.data.frame(counts(dds_norm,normalize=TRUE)),
by="row.names",sort=FALSE)
up_DEG <- subset(res_data, padj < 0.05 & log2FoldChange > 1)
down_DEG <- subset(res_data, padj < 0.05 & log2FoldChange < -1)
write.csv(res_data, "all.csv") #全部基因不篩選德召，做火山圖的背景
write.csv(up_DEG, "up.csv")
write.csv(down_DEG, "down.csv")
（5）判斷歐氏距離，若有異常樣品則不用cooksCutoff汽纤；當(dāng)有上千個(gè)異常值時(shí)也不用：（完全可以不做）
par(mar=c(8,5,2,2))
boxplot(log10(assays(dds_norm)[["cooks"]]), range=0, las=2)
（6）lfcshrink & MA plot
library(apeglm)
resultsNames(dds_norm) #看一下要shrink的維度;shrink數(shù)據(jù)更加緊湊,少了一項(xiàng)stat上岗，但并未改變padj，但改變了foldchange
res_shrink <- lfcShrink(dds_norm, coef="condition_treat2_vs_treat1", type="apeglm") #最推薦apeglm算法;根據(jù)resultsNames(dds)的第5個(gè)維度蕴坪，coef=5肴掷，也可直接""指定;apeglm不allow contrast，所以要指定coef
pdf("MAplot.pdf", width = 6, height = 6)
plotMA(res_shrink, ylim=c(-10,10), alpha=0.1, main="MA plot: ")
dev.off()
（7）火山圖
library(ggplot2)
voldata <-read.csv(file = "all.csv",header = TRUE, row.names =1)
pdf("volcano.pdf", width = 6.13, height = 5.18)
ggplot(data=voldata, aes(x=log2FoldChange,y= -1*log10(padj))) +
geom_point(aes(color=significant)) +
scale_color_manual(values=c("#546de5", "#d2dae2","#ff4757")) +
labs(title="Volcano Plot: ", x=expression(log[2](FC), y=expression(-log[10](padj)))) +
geom_hline(yintercept=1.3,linetype=4) + #反對(duì)數(shù),代表0.05的線
geom_vline(xintercept=c(-1,1),linetype=4) +
theme_bw() + theme(panel.grid = element_blank()) #主次網(wǎng)格線均為空白
dev.off()
3. 心得
- 必須充分了解自己的數(shù)據(jù)背传，分析方法萬(wàn)萬(wàn)千呆瞻，選最合適的才是最好的！
- 任何軟件教程都不如官網(wǎng)
README和文章径玖，遇到問(wèn)題時(shí)要想“一定不是我一個(gè)人”痴脾，善于用Bioconductor論壇。 - 不要鉆牛角尖梳星，完整比完美重要赞赖，任何分析都有不周到，學(xué)會(huì)帶著遺憾繼續(xù)走丰泊，不必十全十美薯定，因?yàn)橛肋h(yuǎn)也不可能做到。
