基因表達分析（上）- 差異表達分析

基因表達

什么是基因表達，如下是來自于維基百科的解釋：

定義
Flow of genetic information

研究方法

定量PCR

這部分我不太懂脸秽，所以就放幾段百度百科和維基百科的定義。

  • 百度百科

定量PCR（即時聚合酶鏈鎖反應，Real-time Polymerase Chain Reaction讹堤，簡稱 Real-time PCR、即時PCR）厨疙，又稱定量即時聚合酶鏈鎖反應（Quantitative real time polymerase chain reaction洲守，簡稱 Q-PCR/qPCR/rt-qPCR、定量即時PCR沾凄、即時定量PCR）梗醇，是一種在DNA擴增反應中，以螢光染劑偵測每次聚合酶鏈鎖反應（PCR）循環(huán)后產(chǎn)物總量的方法技術撒蟀，有廣義概念和狹義概念叙谨。廣義概念的定量PCR技術是指以外參或內(nèi)參為標準，通過對PCR終產(chǎn)物的分析或PCR過程的監(jiān)測保屯，進行PCR起始模板量的定量手负。狹義概念的定量PCR技術（嚴格意義的定量PCR技術）是指用外標法（熒光雜交探針保證特異性）通過監(jiān)測PCR過程（監(jiān)測擴增效率）達到精確定量起始模板數(shù)的目的涤垫，同時以內(nèi)對照有效排除假陰性結果（擴增效率為零）。

  • 維基百科

A real-time polymerase chain reaction (Real-Time PCR), also known as quantitative polymerase chain reaction (qPCR), is a laboratory technique of molecular biology based on the polymerase chain reaction (PCR). It monitors the amplification of a targeted DNA molecule during the PCR, i.e. in real-time, and not at its end, as in conventional PCR. Real-time PCR can be used quantitatively (quantitative real-time PCR), and semi-quantitatively, i.e. above/below a certain amount of DNA molecules (semi quantitative real-time PCR).

優(yōu)點：靈敏性高竟终，準確性高雹姊，通量也還行。一般而言衡楞，RNA-Seq和microassay分析得到的差異表達基因最終也需要通過這種實驗方法進行驗證吱雏。

但是一般適用于驗證實驗，而不是用于探索性實驗瘾境。

microarray<small>基因矩陣</small>

基因芯片的概念在上個世紀80年代就已經(jīng)提出來了歧杏， 被評為1998年度自然科學領域十大進展之一。他的基本原理通過設計專門的短核苷酸作為探針迷守，把這些探針固定在專門的基片表面犬绒，然后用樣本的cDNA進行雜交，根據(jù)雜交信號的強弱來判斷基因表達的程序兑凿。

microarray

但是microarray檢測的基因數(shù)量完全取決于你的探針設計的數(shù)量凯力，而且難以研究mRNA的可變剪切。

RNA-Seq

RNA-Seq是目前基因表達分析最常用的技術礼华。分為以下幾步

  • 分離所有mRNA
  • 逆轉錄mRNA成cDNA
  • 對cDNA測序
  • 比對參考基因組

RNA-Seq實驗設計中的“重復”包括：技術重復和生物學重復
重復是為了檢測組間和組內(nèi)的變異咐鹤，對于假設檢驗至關重要。

  • 技術重復為了估計測量技術（RNA-Seq）的變異圣絮。
  • 生物學重復是為了發(fā)現(xiàn)生物組內(nèi)的變異祈惶。
    簡單的說，兩組的基因表達的變化只有比組內(nèi)變異還大時才能認為時顯著的扮匠。

RNA-Seq的概率分布
相同基因在不同細胞的表達水平服從log-normal(對數(shù)正態(tài))分布捧请，由定量PCR驗證。（注：這與相同細胞不同基因表達的分布不同）但是大多數(shù)基因表達實驗都是用一群細胞棒搜，幾乎沒有相應分布提出疹蛉。

RNA-Seq試驗中，抽樣得到的raw read counts服從泊松分布力麸。并且同一樣本在兩次試驗中的結果不同可款，這稱為shot noise。這種變異在RNA-Seq技術重復間成為Possion noise末盔。

生物學上不同的樣本間的差異服從負二項(negative binomial)分布,有時稱gamma-Poisson分布筑舅。

由于RNA-Seq count數(shù)據(jù)也表現(xiàn)出zero inflation(大量值為0)的特征，所以很難擬合到負二項分布陨舱，所以有文章認為要用Poisson-Tweedie family建模翠拣。

RNA-seq數(shù)據(jù)和microassay在差異表達分析上的區(qū)別：

  • RNA-Seq觀察到的數(shù)據(jù)是抽樣過程中產(chǎn)生的離散（discrete）count形式。也就是說總體是恒定的游盲，表達量越高的基因在抽樣結果中所占的比例越大误墓。表達量低的基因可能即便有也無法被檢測出來蛮粮。當然，重新對相同文庫進行測序谜慌，還是有可能找到更多表達的轉錄本

  • microassay檢測的是熒光信號的連續(xù)度量然想。由于使用固定的核酸序列去雜交，所以不是一種“零和游戲”欣范，只要能雜交变泄，就能被檢測。（但如果沒有設計相應的引物恼琼，就不能檢測到可能的基因）

研究意義

1.在不同背景下比較mRNA水平

  • 同一物種妨蛹，不同組織：研究基因在不同部分的表達情況
  • 同一物種，同一組織：研究基因在不同處理下晴竞，不同條件下的表達變化
  • 同一組織蛙卤，不同物種：研究基因的進化關系
  • 時間序列實驗： 基因在不同時期的表達情況與發(fā)育的關系

2.基因分類： 找到細胞特異，疾病相關噩死，處理相關的基因表達模式颤难，用于診斷疾病和預測等

3.基因網(wǎng)絡和通路： 基因在細胞活動中的功能，基因間的相互作用已维。

公共數(shù)據(jù)庫

當然行嗤，如果你要研究一個基因的功能時，不要先急著去花錢找公司測序衣摩，先去一些基因表達公共數(shù)據(jù)庫找找看：

  1. http://www.ebi.ac.uk/arrayexpress/
  2. https://www.ncbi.nlm.nih.gov/geo/

差異表達(differential expression,DE)基因分析

通過研究基因的差異表達昂验，我們可以發(fā)現(xiàn)

  • 細胞特異性的基因捂敌；
  • 發(fā)育階段特異性的基因艾扮；
  • 疾病狀態(tài)相關的基因；
  • 環(huán)境相關的基因占婉；
  • ...

基本方法就是以生物學意義的方式計算基因表達量泡嘴，然后通過統(tǒng)計學分析表達量尋找具有統(tǒng)計學顯著性差異的基因，從而

  • 選擇合適的基因
  • 衡量結果的可靠性

分析方法

尋找差異表達基因有三種方式：
第一種是計算Fold change（倍數(shù)變化）逆济，十分簡單粗暴的方法酌予，計算方法如下：

  • E = mean(group1) B=mean(group2)
  • FC = (E-B) / min(E,B)

說人話就是，基因A和基因B的平均值之差與兩者中較小的比值奖慌。選擇2-3倍的基因作為結果（為什么是2-3倍抛虫，就是大家約定俗成）。

fold change

但是簡單粗暴的用2到3倍作為閾值简僧，對于低表達的基因建椰，3倍也是噪音，那些高表達的基因岛马，1.1倍都是生物學顯著了棉姐。更重要的沒有考慮到組內(nèi)變異屠列，沒有統(tǒng)計學意義。所以發(fā)文章肯定這個圖只能作為附錄了伞矩。

第二種就是統(tǒng)計檢驗笛洛，寫文章的時候總需要給出一個p值告訴主編這個結果可信的（雖然p值也存在爭論）。
復習一下：p值指的碰巧是拒絕零假設機會乃坤。P值越大假陽性越低苛让，同時真實結果也可能會剔除。

注： 基因表達分析的零假設是： 基因在不同處理下的表達量相同湿诊。

對于基因芯片的數(shù)據(jù)而言蝌诡，由于樣本服從正態(tài)分布，所以可以用t-test（雙處理）或anova分析（多處理以上）枫吧。

T檢驗適用于只有兩個處理的實驗設計浦旱，如植物葉片在相同處理第一天和第二天的基因表達差異。


T-Test 實驗設計

進行T-test檢驗時要注意：是雙尾檢驗（存在差異）還是單尾檢驗（顯著性上調(diào)或下降）九杂，兩個樣本的總體是不是等方差（標準T檢驗還是Welch’s test）

如果存在多于兩個處理（條件）颁湖，就需要用到ANOVA分析了。ANOVA分析能主要是研究結果之間的差異是如何引起的例隆，具體請移步到我寫方差分析教程甥捺。
對于基因表達而言，研究目標是镀层，對于同一個基因而言镰禾，他們之間的差異是處理不同造成，還是因為系統(tǒng)誤差造成唱逢。

單因素
雙因素

當然你可以研究吴侦，不同基因的表達差異是由因為處理不同，還是基因不同坞古，還是系統(tǒng)誤差备韧，還是其中一些的交互作用。

上面都是針對基因芯片的樣本服從正態(tài)分布進行的統(tǒng)計檢驗』痉悖現(xiàn)在的RNA-Seq织堂，它的抽樣過程是離散的，結果是count奶陈，服從泊松分布易阳，樣本間的差異是服從負二向分布，顯然不能按照上述方法分析吃粒。

方差分析（ANOVA)和線性回歸分析(regression)都是同一時期發(fā)展的兩套緊密相連的理論潦俺。方差分析考量的是離散型自變量(因子）對連續(xù)型應變量（響應變量）的模型分析，而線性回歸分析只要求響應變量是連續(xù)的，對于自變量無要求黑竞。如果響應變量不是連續(xù)型分布捕发，就要使用更加一般化的廣義線性模型（generalized linear model),通過一個連接函數(shù)變換響應變量期望，將響應變量的期望與自變量建立線性關系很魂。

因此扎酷，我們可以用廣義線性模型去分析RNA-seq前期分析得到的離散型結果（count）

方差分析一般用于分析有計劃的實驗結果，比如說不同處理下的水稻產(chǎn)量遏匆》òぃ回歸分析一般分析沒有計劃的數(shù)據(jù)，比如說你可以找到大量體檢的數(shù)據(jù)幅聘，只分析其中性別和身高對體重的影響凡纳。所以兩者各有側重，不要拿大炮轟蚊子帝蒿。

統(tǒng)計檢驗相對于fold change具有統(tǒng)計意義荐糜，不需要參考樣本，需要處理隨機取樣葛超。但是需要重復（ANOVA推薦4-10重復）暴氏，但由于資金和材料等原因，不一定能夠滿足绣张。此外答渔，對于1000個基因，就要做1000次ANOVA或t-test侥涵，最后的p值會有一定的假陽性沼撕，因此要做p值矯正（FDR）篩選。

注：推薦在統(tǒng)計檢驗前過濾表達量低芜飘，也就如果一個基因在所有樣本中count均低于某一閥值务豺，請在分析前剔除。這個閥值也是約定俗成燃箭，一般設置為3.

第三種：Fold Change + 統(tǒng)計檢驗冲呢。說一比較尷尬的事情，在統(tǒng)計檢驗中你找到越多的差異表達基因招狸，在p值矯正之后，你反而找不到差異表達基因邻薯。也就是說裙戏，如果在結果中存在大量濫竽充數(shù)的所謂的DE基因，那么在嚴格的p值矯正篩選后厕诡，反而會誤刪真實的DE基因累榜。

因此在p值矯正之前，你先要手動剔除一部分明顯就是假陽性的DE基因。這個步驟就需要用到前面的fold-change分析壹罚。
我們可以通過火山圖來看看如何確定區(qū)間

火山圖

實戰(zhàn)演練

為了方便理解葛作，我們選用同一組數(shù)據(jù)進行實際操作。默認你懂基本的R語言操作猖凛，如安裝R包赂蠢，查看幫助文件等。
正式學習之前辨泳，先感謝一下Bioconductor虱岂。 根據(jù)他們提供流程，我學習到如何進行RNA-Seq分析菠红。

數(shù)據(jù)介紹和下載

用到的數(shù)據(jù)來自于一篇文獻“A pooled shRNA screen for regulators of primary mammary stem and progenitor cells identifies roles for Asap1 and Prox1”第岖。

實驗目的:

The aim of this study is to determine the absolute and relative expression levels of mRNA transcripts across multiple flow cytometrically sorted epithelial cell types including freshly isolated CD24+CD29hi mammary stem cell-enriched basal cells (MaSC/basal), CD24+CD29loCD61+ luminal progenitor-enriched (LP) and the CD24+CD29loCD61- mature luminal-enriched (ML) cell populations. Additionally, a comparison between these primary cell types and cultured MaSC/Basal-derived mammosphere cells (mammosphere) and the CommaD-βGeo (CommaDβ) cell line was performed.

實驗設計：

Total RNA was extracted and purified from sorted luminal or basal populations from the mammary glands of female virgin 8- to 10-week-old FVB/N mice (3 independent samples for population), MaSC/Basal cells cultured for 1 week under mammsophere conditions and CommaDβ cells grown under maintenance conditions (Deugnier et al. 2006) and their transcriptomes analysed by RNA-Seq. |

根據(jù)文章的介紹，數(shù)據(jù)是100 bp的單端read试溯，用Rsubread比對到mouse reference genome(mm10), 然后使用featureCounts統(tǒng)計每個基因的count數(shù)镊尺。然后用TMM進行標準化，轉換成log2 counts per million.最后用limma包對每個樣本每個基因的平均表達值以觀察水平權重的線性模型進行擬合确封，并用T檢驗找到不同群體的差異表達基因熬甫。以FDR + log2-fold-change對基因排序。

我們的任務是下載他們featureCounts的結果试读，不用limma,改用DESeq2分析杠纵。

# install.packages("R.utils")
library("R.utils")
url <- "https://www.ncbi.nlm.nih.gov/geo/download/?acc=GSE63310&format=file"
utils::download.file(url, destfile="GSE63310_RAW.tar", mode="wb") 
utils::untar("GSE63310_RAW.tar", exdir = ".")
files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", "GSM1545538_purep53.txt",
           "GSM1545539_JMS8-2.txt", "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt",
           "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", "GSM1545545_JMS9-P8c.txt")

for(i in paste(files, ".gz", sep=""))
  R.utils::gunzip(i, overwrite=TRUE)

數(shù)據(jù)讀取

先定一個包含所有文件的向量，方便后續(xù)調(diào)用钩骇。

files <- c("GSM1545535_10_6_5_11.txt", "GSM1545536_9_6_5_11.txt", 
   "GSM1545538_purep53.txt", "GSM1545539_JMS8-2.txt", 
   "GSM1545540_JMS8-3.txt", "GSM1545541_JMS8-4.txt", 
   "GSM1545542_JMS8-5.txt", "GSM1545544_JMS9-P7c.txt", 
   "GSM1545545_JMS9-P8c.txt")

查看下文件的數(shù)據(jù)格式,分別記錄了每個基因的EntrezID比藻，基因長度和數(shù)量

read.delim(files[1], nrow=5)
#   EntrezID GeneLength Count
#1    497097       3634     1
#2 100503874       3259     0
#3 100038431       1634     0
#4     19888       9747     0
#5     20671       3130     1

之后我們需要用DESeqDataSetFromMatrix為DESeq2提供數(shù)據(jù),命令形式如下：

library("DEseq2")
DESeqDataSetFromMatrix(countData, colData, design, tidy = FALSE,ignoreRank = FALSE, ...)

其中countData存放counts數(shù)據(jù)，colData存放樣本信息的數(shù)據(jù)倘屹，design就是實驗設計银亲。

首先從各個count matrix文件中讀取count，基因長度部分可以舍棄纽匙，因為DESeq2只需要為標準化的count數(shù)據(jù)务蝠，不需要提供基因長度信息。
邏輯就是分別讀取每一個文件的count列烛缔，然后賦予文件名馏段。

sample1 <- read.delim(files[1],header = T)[,c(1,3)]
count.table <- data.frame(sample1)
for ( f in files[2:length(files)]){
  column <- read.delim(f,header=T,row.names = 1)[,2]
  count.table <- cbind(count.table, column)
}
head(count.table)
rownames(count.table) <- count.table[,1]
count.table <- count.table[-1]
samplenames <- substring(files, 12, nchar(files)-4)
colnames(count.table) <- samplenames

然后要提供colData,其中colData存放列名要和countData的行名相同。

# colData
condition <- as.factor(c("LP", "ML", "Basal", "Basal", "ML", "LP", 
                         "Basal", "ML", "LP"))
lane <- as.factor(rep(c("L004","L006","L008"), c(3,4,2)))
colData <- data.frame(condition = condition, lane=lane,row.names = samplenames)

all(rownames(colData) %in% colnames(count.table))

最后導入數(shù)據(jù):

dds <- DESeqDataSetFromMatrix(countData = count.table, 
                              colData = colData,
                              design = ~ condition )

數(shù)據(jù)預處理

預過濾

雖然DESeq2會自動屏蔽那些低count的基因践瓷，但是剔除那些幾乎不存在基因的部分能夠提高運行速度院喜。

dds <- dds [ rowSums(counts(dds)) > 1, ]

rlog和方差齊性轉換

許多常見的多維數(shù)據(jù)探索性分析的統(tǒng)計分析方法，例如聚類和主成分分析要求晕翠，在那些同方差性的數(shù)據(jù)表現(xiàn)良好喷舀。所謂的同方差性就是雖然平均值不同，但是方差相同。

但是對于RNA-Seq count數(shù)據(jù)而言硫麻，當均值增加時爸邢，方差期望也會提高。也就說直接對count matrix或標準化count（根據(jù)測序深度調(diào)整）做PCA分析拿愧，由于高count在不同樣本間的絕對差值大杠河，也就會對結果有很大影響。簡單粗暴的方法就是對count matrix取log后加1赶掖。這個1也是約定俗成感猛，看經(jīng)驗了。

隨便舉個栗子看下效果：

lambda <- 10^seq(from = -1, to = 2, length = 1000)
cts <- matrix(rpois(1000*100, lambda), ncol = 100)
library("vsn")
meanSdPlot(cts, ranks = FALSE)
平滑前
log.cts.one <- log2(cts + 1)
meanSdPlot(log.cts.one, ranks = FALSE)
平滑后

DESeq2為count數(shù)據(jù)提供了兩類變換方法奢赂，使得不同均值的方差趨于穩(wěn)定：regularized-logarithm transformation or rlog(Love, Huber, and Anders 2014)和variance stabilizing transformation(VST)(Anders and Huber 2010)用于處理含有色散平均趨勢負二項數(shù)據(jù)陪白。

到底用啥：數(shù)據(jù)集小于30 -> rlog，大數(shù)據(jù)集 -> VST膳灶。還有這個處理過程不是用于差異檢驗的咱士，在DESeq分析中會自動選擇最合適的所以你更不需要糾結了，記得用raw count轧钓。

rld <- rlog(dds, blind = FALSE)
head(assay(rld), 3)
vsd <- vst(dds, blind = FALSE)
head(assay(vsd), 3)

處理的結果如下：

library("dplyr")
library("ggplot2")

dds <- estimateSizeFactors(dds)

df <- bind_rows(
  as_data_frame(log2(counts(dds, normalized=TRUE)[, 1:2]+1)) %>%
    mutate(transformation = "log2(x + 1)"),
  as_data_frame(assay(rld)[, 1:2]) %>% mutate(transformation = "rlog"),
  as_data_frame(assay(vsd)[, 1:2]) %>% mutate(transformation = "vst"))

colnames(df)[1:2] <- c("x", "y")  

ggplot(df, aes(x = x, y = y)) + geom_hex(bins = 80) +
  coord_fixed() + facet_grid( . ~ transformation)
不同轉換的分布

結果就是轉換后更加集中了序厉。

樣本距離

RNA-Seq分析第一步通常是評估樣本間的總體相似度。

  • 那些樣本最接近
  • 那些樣本差異較大
  • 這與實驗設計預期符合么

這里使用R內(nèi)置的 dist 計算 rlog變換數(shù)據(jù)的Euclidean distance毕箍，然后用pheatmap根據(jù)結果畫熱圖

sampleDists <- dist(t(assay(rld)))
library("pheatmap")
library("RColorBrewer")
sampleDistMatrix <- as.matrix( sampleDists )

colors <- colorRampPalette( rev(brewer.pal(9, "Blues")) )(255)
pheatmap(sampleDistMatrix,
         clustering_distance_rows = sampleDists,
         clustering_distance_cols = sampleDists,
         col = colors)
heatmap

同樣的可以用Poisson Distance (Witten 2011)計算距離弛房，計算方式如下：

library("PoiClaClu")
poisd <- PoissonDistance(t(counts(dds)))

PCA分析

還有一種可視化樣本-樣本距離的方法就是主成分分析。PCA分析我打算找點資料好好理解之后再寫而柑，這個說下有這個方法文捶。
DESeq2提供了專門的方法用于作圖，還很好看呢媒咳！

plotPCA(rld,intgroup=c("condition"))
PCA

能夠明顯的發(fā)現(xiàn)不同處理的距離離得很遠粹排。

差異表達基因分析（DEA）

在我們定義好DESeqDataSe 就可以方便的調(diào)用DESeq分析.

dds <- DESeq(dds)

在幾行輸出后信息后，分析就完成了涩澡，更多具體參數(shù)可以用?DESeq查看手冊

構建結果表格

如果直接調(diào)用results顽耳，只會提取出design matrix最后兩個變量（這里是ML vs Basal）的 log2 fold changes and p values等結果。

> results(dds)
log2 fold change (MAP): condition ML vs Basal 
Wald test p-value: condition ML vs Basal 
DataFrame with 27179 rows and 6 columns
             baseMean log2FoldChange     lfcSE        stat       pvalue         padj
            <numeric>      <numeric> <numeric>   <numeric>    <numeric>    <numeric>
497097    118.7522256     -7.0662875 0.5802473 -12.1780608 4.068401e-34 8.556128e-33
100503874   1.3308906     -2.6371037 1.3370584  -1.9723174 4.857338e-02 8.946685e-02
100038431   0.0843771     -0.1541107 1.2443948  -0.1238439 9.014388e-01           NA
19888       1.4833878      2.9503142 1.3351543   2.2097179 2.712475e-02 5.385925e-02
20671      26.4317752     -1.5256766 0.7754853  -1.9673829 4.913908e-02 9.034136e-02
...               ...            ...       ...         ...          ...          ...
100861837   368.84237     0.33334299 0.3226665   1.0330883    0.3015626    0.4139846
100861924    22.79272    -1.25169702 0.8307249  -1.5067528    0.1318740    0.2107990
170942      852.61538    -0.07612745 0.2399318  -0.3172879    0.7510252    0.8235139
100861691     0.00000             NA        NA          NA           NA           NA
100504472     0.00000             NA        NA          NA           NA           NA

實際上results有如下眾多參數(shù)

results(object, contrast, name, lfcThreshold = 0,
  altHypothesis = c("greaterAbs", "lessAbs", "greater", "less"),
  listValues = c(1, -1), cooksCutoff, independentFiltering = TRUE,
  alpha = 0.1, filter, theta, pAdjustMethod = "BH", filterFun,
  format = c("DataFrame", "GRanges", "GRangesList"), test, addMLE = FALSE,
  tidy = FALSE, parallel = FALSE, BPPARAM = bpparam(), ...)

比如可以指定比較對象Basal和PL,可以用mcols查看結果存儲的元數(shù)據(jù)妙同，了解列名的含義射富。

res <- results(dds,contrast = c("condition","Basal","LP"))
mcols(res, use.names = TRUE)

DataFrame with 6 rows and 2 columns
                       type                                   description
                <character>                                   <character>
baseMean       intermediate     mean of normalized counts for all samples
log2FoldChange      results log2 fold change (MAP): condition Basal vs LP
lfcSE               results         standard error: condition Basal vs LP
stat                results         Wald statistic: condition Basal vs LP
pvalue              results      Wald test p-value: condition Basal vs LP
padj                results                          BH adjusted p-valuess

其中l(wèi)fcSE是Log2 fold change的標準誤。其他項都比較好理解渐溶。
最后還可以看一些總結性的內(nèi)容

summary(res)

out of 22026 with nonzero total read count
adjusted p-value < 0.1
LFC > 0 (up)     : 5528, 25% 
LFC < 0 (down)   : 5274, 24% 
outliers [1]     : 43, 0.2% 
low counts [2]   : 2953, 13% 
(mean count < 1)
[1] see 'cooksCutoff' argument of ?results
[2] see 'independentFiltering' argument of ?results

在limma分析結果分別是4127辉浦，4298，稍微多了那么2000個基因茎辐。其他結果也基本上都了2000個。看起來對于22000多個基因而言拖陆，差異不算太大弛槐。

當然我們還可以通過

  • 降低假陽性率（FDR， false discovery rate)閾值
  • 提高log2 fold change的閾值
> res0.5 <- results(dds, contrast = c("condition","Basal","LP"),alpha=0.05)
> table(res0.5$padj < 0.05)

FALSE  TRUE 
 8860  9750 
> resLFC1 <- results(dds, lcontrast = c("condition","Basal","LP"),fcThreshold=1)
> table(resLFC1$padj < 0.1)

FALSE  TRUE 
 8916 10958

于是乎結果就和limma差不多了依啰。

多重試驗

在高通量數(shù)據(jù)分析中乎串，我們通常不是用p值來拒絕原假設，更多是用來進行多重試驗矯正速警。

為什么要做多重試驗矯正叹誉？

如果對一個基因，我有99%的把握認為是差異表達的闷旧，也就是說1%的可能是錯的长豁。那么假設有10000個基因，按照數(shù)學期望忙灼，有100個是假的匠襟。因此為了保證多重試驗結果的可靠性要對結果的p-value做矯正。

DESeq2用的Benjamini-Hochberg (BH) adjustment (Benjamini and Hochberg 1995)该园，計算矯正后的P-value,用于回答如下問題：

如果我們選擇了所有小于或等于矯正p-value閾值的所有顯著性基因酸舍，假陽性比例（ false discovery rate, FDR）是多少？

FDR在高通量試驗中是比較有用的統(tǒng)計值里初，由于我們通常關注一類自己感興趣的基因啃勉，所以我們需要設置一個假陽性上限。

我們從結果中以1%作為顯著性双妨，分別找出顯著性上調(diào)和下調(diào)的基因淮阐。

regSig <- subset(res, padj < 0.01)
# Up
head(resSig[ order(resSig$log2FoldChange), ])
# Down
head(resSig[ order(resSig$log2FoldChange, decreasing = TRUE), ])

可視化

人類對圖像比較敏感，對文字比較遲鈍斥难。所以我們需要一些比較好看的圖來放到文章中解釋說明枝嘶。

最簡單的就是Counts plot，看看特定基因的count數(shù)量哑诊。

topGene <- rownames(res)[which.min(res$padj)]
plotCounts(dds, gene = topGene, intgroup=c("condition"))
count plot

可以用MA-plot（也叫mean-difference plot,Bland-Altman plot）了解模型（如所有基因在不同處理比較的結果）的估計系數(shù)的分布

res <- lfcShrink(dds, contrast=c("contion","ML","PL"), res=res)
plotMA(res, ylim = c(-5, 5))

更加喜聞樂見的是基因聚類所提供的熱圖展示群扶。我們可以找前20個樣本件差異比較大，然后看他們在不同樣本間的表達情況镀裤。

library("genefilter")
topVarGenes <- head(order(rowVars(assay(rld)), decreasing = TRUE), 20)
mat  <- assay(rld)[ topVarGenes, ]
mat  <- mat - rowMeans(mat)
anno <- as.data.frame(colData(rld)[, c("condition","lane")])
pheatmap(mat,annotation_col = anno)
heatmap

大致可以發(fā)現(xiàn)同一組的基因顏色是相同的竞阐，也就是說表達量相近。

結語

找到差異表達的基因只是第一步暑劝，后續(xù)還需要對這些基因進行進一步的分析骆莹，如

  • 基因富集分析
  • 主成分分析
  • 聚類分析

而這些內(nèi)容就是我接下來的學習重點，也是下次更新文章的主題担猛。

最后編輯于 ：
?著作權歸作者所有,轉載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末幕垦，一起剝皮案震驚了整個濱河市丢氢，隨后出現(xiàn)的幾起案子，更是在濱河造成了極大的恐慌先改，老刑警劉巖疚察，帶你破解...
    沈念sama閱讀 206,013評論 6贊 481
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件，死亡現(xiàn)場離奇詭異仇奶，居然都是意外死亡貌嫡，警方通過查閱死者的電腦和手機，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 88,205評論 2贊 382
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進店門该溯，熙熙樓的掌柜王于貴愁眉苦臉地迎上來岛抄，“玉大人，你說我怎么就攤上這事狈茉》蛲郑” “怎么了？”我有些...
    開封第一講書人閱讀 152,370評論 0贊 342
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵论皆，是天一觀的道長益楼。 經(jīng)常有香客問我，道長点晴，這世上最難降的妖魔是什么感凤？ 我笑而不...
    開封第一講書人閱讀 55,168評論 1贊 278
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任，我火速辦了婚禮粒督，結果婚禮上陪竿，老公的妹妹穿的比我還像新娘。我一直安慰自己屠橄，他們只是感情好族跛，可當我...
    茶點故事閱讀 64,153評論 5贊 371
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布。 她就那樣靜靜地躺著锐墙，像睡著了一般礁哄。 火紅的嫁衣襯著肌膚如雪。 梳的紋絲不亂的頭發(fā)上溪北，一...
    開封第一講書人閱讀 48,954評論 1贊 283
  • 城市分裂傳說
    那天桐绒，我揣著相機與錄音，去河邊找鬼之拨。 笑死茉继，一個胖子當著我的面吹牛，可吹牛的內(nèi)容都是我干的蚀乔。 我是一名探鬼主播烁竭，決...
    沈念sama閱讀 38,271評論 3贊 399
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼，長吁一口氣：“原來是場噩夢啊……” “哼吉挣！你這毒婦竟也來了派撕？” 一聲冷哼從身側響起婉弹，我...
    開封第一講書人閱讀 36,916評論 0贊 259
  • 萬榮殺人案實錄
    序言：老撾萬榮一對情侶失蹤，失蹤者是張志新（化名）和其女友劉穎腥刹，沒想到半個月后马胧，有當?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體汉买，經(jīng)...
    沈念sama閱讀 43,382評論 1贊 300
  • ?護林員之死
    正文 獨居荒郊野嶺守林人離奇死亡衔峰，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 35,877評論 2贊 323
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年，在試婚紗的時候發(fā)現(xiàn)自己被綠了蛙粘。 大學時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片垫卤。...
    茶點故事閱讀 37,989評論 1贊 333
  • 活死人
    序言：一個原本活蹦亂跳的男人離奇死亡，死狀恐怖出牧，靈堂內(nèi)的尸體忽然破棺而出穴肘，到底是詐尸還是另有隱情，我是刑警寧澤舔痕，帶...
    沈念sama閱讀 33,624評論 4贊 322
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布评抚，位于F島的核電站，受9級特大地震影響伯复，放射性物質(zhì)發(fā)生泄漏慨代。R本人自食惡果不足惜，卻給世界環(huán)境...
    茶點故事閱讀 39,209評論 3贊 307
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一啸如、第九天 我趴在偏房一處隱蔽的房頂上張望侍匙。 院中可真熱鬧，春花似錦叮雳、人聲如沸想暗。這莊子的主人今日做“春日...
    開封第一講書人閱讀 30,199評論 0贊 19
  • 一樁弒父案帘不，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽说莫。三九已至，卻和暖如春寞焙，著一層夾襖步出監(jiān)牢的瞬間储狭，已是汗流浹背。 一陣腳步聲響...
    開封第一講書人閱讀 31,418評論 1贊 260
  • 情欲美人皮
    我被黑心中介騙來泰國打工棺弊， 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留晶密，地道東北人。 一個月前我還...
    沈念sama閱讀 45,401評論 2贊 352
  • 代替公主和親
    正文 我出身青樓模她，卻偏偏與公主長得像稻艰，于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子侈净，可洞房花燭夜當晚...
    茶點故事閱讀 42,700評論 2贊 345

推薦閱讀更多精彩內(nèi)容