菜??的第一次正式上手R，相當(dāng)于翻譯了 edgeR 這個(gè)包的 UserGuide，跑了第一個(gè) case
眼看寫了這么多字，還是發(fā)一下比較好，萬(wàn)一有人看呢

edgeR簡(jiǎn)介

edgeR 可以應(yīng)用于任何可產(chǎn)生基因組特征數(shù)據(jù)（read counts）的技術(shù)，能夠?yàn)?RNA-seq 實(shí)驗(yàn)中評(píng)估差異表達(dá)、ChIP-seq 實(shí)驗(yàn)中差異標(biāo)記提供統(tǒng)計(jì)程序。該R包具備適用于多組實(shí)驗(yàn)的精準(zhǔn)統(tǒng)計(jì)法，同時(shí)還具備廣義線性模型（glms）的統(tǒng)計(jì)學(xué)方法——適用于不同復(fù)雜程度的多因素實(shí)驗(yàn)。有時(shí)人們將前者稱為 classic edgeR，將后者稱為 glm edgeR。然而上述兩種方法是互補(bǔ)的，并且時(shí)常在數(shù)據(jù)分析中被結(jié)合使用。大多數(shù) glm 函數(shù)可以通過函數(shù)名稱中的 "glm" 識(shí)別，這類函數(shù)可利用似然比檢驗(yàn)或擬似然F檢驗(yàn)檢測(cè)差異表達(dá)。edgeR 的功能的一個(gè)重要特點(diǎn)是，不論 classic 和 glm 兩種方法，都屬于經(jīng)驗(yàn)貝葉斯方法，從而能夠在實(shí)驗(yàn)只具有最小水平的生物學(xué)重復(fù)時(shí)，依然能夠判斷出基因特異的生物學(xué)差異。edgeR 可應(yīng)用于不同水平的差異表達(dá)，如基因、外顯子、轉(zhuǎn)錄本、標(biāo)簽等，分析外顯子水平時(shí)可輕易地檢測(cè)出可變剪切或異構(gòu)體特異性差異表達(dá)。

案例分析

RNA-seq —— 口腔腫瘤 vs 對(duì)應(yīng)的正常組織 (RNA-Seq of oral carcinomas vs matched normal)

分析的目的是檢測(cè)腫瘤和正常組織對(duì)比下差異表達(dá)的基因，這個(gè)案例可以體現(xiàn) edgeR 中 GLM 法的工作能力。

1.edgeR 安裝

if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

BiocManager::install("edgeR")

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

下載文章 Tumor Transcriptome Sequencing Reveals Allelic Expression Imbalances Associated with Copy Number Alterations 中的 Table S1。

• 去除多余列項(xiàng)、重命名列名

  > rawdata <- read.delim("TableS1.txt", check.names=FALSE, stringsAsFactors=FALSE)
  

  > rawdata <- rawdata[,c(1:9)]
  

  > names(rawdata)[1:9]<-c("RefSeqID","Symbol","NbrOfExons","8N","8T","33N","33T","51N","51T")
  

  這里不小心把表格的一二行的行數(shù)弄沒了，后面需要注意數(shù)字

  結(jié)果到DEGList這一步時(shí)出現(xiàn)報(bào)錯(cuò)“The count matrix is a data.frame instead of a matrix and the first 6 columns are non-numeric.” 用了as.matrix()后表格里所有內(nèi)容又都變成了factor...??

  **所以，最后還是直接在excel表格里解決了行名??
  更新：其實(shí)這個(gè)很好解決，然這是第一次正式，沒錯(cuò)“正式”，用R....

• 加載edgeR, 用DEGList構(gòu)建列表

  > library(edgeR) Loading 
  > y <- DGEList(counts=rawdata[,4:9], genes=rawdata[,1:3]) 
  

3.注釋

• 篩選轉(zhuǎn)錄本

  這篇文章發(fā)表已經(jīng)是幾年前了，數(shù)據(jù)中一些 RefSeq ID 可能和現(xiàn)在一般使用的 RefSeq ID 有出入，所以只需要保留org.HS.eg.db中有的、具有 NCBI 注釋的那部分轉(zhuǎn)錄本。

  > library(org.Hs.eg.db) 
  > idfound <- y$genes$RefSeqID %in% mappedRkeys(org.Hs.egREFSEQ) 
  > y <- y[idfound,] 
  > dim(y) 
  

  可以看出有15548個(gè)基因

• 將 Entrez Gene ID 加入注釋

  > egREFSEQ <- toTable(org.Hs.egREFSEQ) 
  > head(egREFSEQ) 
  

  > m <- match(y$genes$RefSeqID, egREFSEQ$accession) 
  > y$genes$EntrezGene <- egREFSEQ$gene_id[m] 
  

• 利用 Entrez Gene ID 更新 gene symbol

  > egSYMBOL <- toTable(org.Hs.egSYMBOL) 
  > head(egSYMBOL) 
  

  > m <- match(y$genes$EntrezGene, egSYMBOL$gene_id) 
  > y$genes$Symbol <- egSYMBOL$symbol[m] 
  

5.篩選及歸一化 （Normalization）

• 篩選出 count 數(shù)最多的轉(zhuǎn)錄本

  每個(gè) gene symbol 保留一個(gè)轉(zhuǎn)錄本

  > o <- order(rowSums(y$counts), decreasing=TRUE) > y <- y[o,] 
  > d <- duplicated(y$genes$Symbol) 
  > y <- y[!d,] 
  > nrow(y) 
  

• 重新計(jì)算文庫(kù)大小

  > y$samples$lib.size <- colSums(y$counts) 
  

• 將 Use Entrez Gene ID 設(shè)為行名

  > rownames(y$counts) <- rownames(y$genes) <- y$genes$EntrezGene 
  > y$genes$EntrezGene <- NULL 
  

• TMM歸一化 (trimmed mean of M-values normalization)

  > y <- calcNormFactors(y) 
  > y$samples 
  

6.數(shù)據(jù)挖掘 (Data Exploration)

• 利用plotMDS繪制MDS (Multidimensional scaling) 圖

  首先應(yīng)該分析樣本的離群值和其他關(guān)系，函數(shù)plotMDS能夠生成樣本之間距離與生物差異系數(shù) (biological coefficient of variation, BCV) 相對(duì)應(yīng)的圖表。

  > plotMDS(y) 
  

  圖中橫坐標(biāo)將腫瘤樣本和正常樣本區(qū)分開來(lái)（左T右N），縱坐標(biāo)則大致對(duì)應(yīng)了患者編碼，這些應(yīng)證了樣本的配對(duì)特性?？梢钥闯瞿[瘤樣本的分布比正常樣本更不均勻。

7.設(shè)計(jì)矩陣

• 檢測(cè)差異表達(dá)

  在擬合負(fù)二項(xiàng) GLM 之前，需要根據(jù)實(shí)驗(yàn)設(shè)計(jì)優(yōu)化設(shè)計(jì)矩陣。接下來(lái)需要在同一患者的腫瘤和正常組織樣本中檢測(cè)差異表達(dá)，即調(diào)整患者之間的差異。統(tǒng)計(jì)學(xué)意義上，這是把患者作為區(qū)組因子的一種可加線性模型。

  > Patient <- factor(c(8,8,33,33,51,51)) 
  > Tissue <- factor(c("N","T","N","T","N","T")) 
  > data.frame(Sample=colnames(y),Patient,Tissue)
  

  > design <- model.matrix(~Patient+Tissue) 
  > rownames(design) <- colnames(y)
  > design
  

  這種可加模型適用于成對(duì)設(shè)計(jì)或具批次效應(yīng) (batch effect) 的實(shí)驗(yàn)。

8.估計(jì)數(shù)據(jù)分布

• 估計(jì)數(shù)據(jù)組的負(fù)二項(xiàng)分布

  > y <- estimateDisp(y, design, robust=TRUE) 
  > y$common.dispersion 
  

  公共離散值 (common dispersion) 的平方根，即生物差異 (biological variation) 的差異系數(shù) (coefficient of variation)，為0.4。

  > plotBCV(y) 
  

9.差異表達(dá)

• 擬合基因?qū)用?GLM

  > fit <- glmFit(y, design) 
  

• 腫瘤樣本 vs 正常組織，做似然比檢驗(yàn)，列出差異最為顯著的基因

  > lrt <- glmLRT(fit) 
  > topTags(lrt) 
  

  這個(gè)檢驗(yàn)類似配對(duì)t檢驗(yàn)。排列靠前的標(biāo)簽/基因的p值、FDR都很小，差異倍數(shù) (FC) 很大。

  進(jìn)一步地??

• 單個(gè)樣本中這些基因的每百萬(wàn) count 數(shù) (counts-per-million)

  ps.根據(jù)p值排序

  > o <- order(lrt$table$PValue) 
  > cpm(y)[o[1:10],] 
  

  可以看出，這些基因在腫瘤樣本和正常組織中的變化，在三個(gè)患者中是一致。

• 以 FDR（錯(cuò)誤發(fā)現(xiàn)率）為5%，查看差異表達(dá)基因的數(shù)量

  > summary(decideTests(lrt)) 
  

• 繪制 log-FC 與 log-cpm 均值 - 差異圖 (mean-difference plot)
  圖中的藍(lán)色橫線代表2倍差異 (2-fold changes)。

  > plotMD(lrt) 
  > abline(h=c(-1, 1), col="blue") 
  

10.GO分析

• 針對(duì)生物學(xué)過程 (biological process, BP)做GO分析

  > go <- goana(lrt) 
  > topGO(go, ont="BP", sort="Up", n=30, truncate=30) 
  

  腫瘤中表達(dá)上調(diào)的基因可能與細(xì)胞分化、細(xì)胞遷移、組織形態(tài)發(fā)生相關(guān)。

Setup Info