?

請一定看這里：寫下來只是為了記錄一些自己的實踐，當然如果能對你有所幫助那就更好了，歡迎大家和我交流

三者區別

?

三者區別

差異分析流程:

1 初始數據

2 標準化(normalization)：DESeq、TMM等

為什么要標準化？
消除文庫大小不同，測序深度對差異分析結果的影響
怎樣標準化？
找到一個能反映文庫大小的因子，利用這個因子對rawdata進行標準化

3 根據模型檢驗求p value：泊松分布(poisson distribution)、負二項式分布(NB)等

4 多重假設得FDR值：

檢驗方法：Wald、LRT
多重檢驗：BH

5 差異基因篩選：pvalue、padj

---

💥💥💥一、 edgeR的使用💥💥💥

因為目前沒有合適的數據，所以數據來源于這里?參考這篇：劉堯科學網博客

0. 前期工作

用到的gene.txt文件，內容如下

gene.txt文件內容

c1表示為一組，c2表示為另一組，.后為第幾個樣本
讀取數據并設置分組，要保證樣本名稱和分組名稱的順序是一一對應的。

#加載edgeR包 library(edgeR) #讀進來文件 targets <- as.matrix(read.delim("你的路徑/gene.txt", sep = '\t', row.names = 1)) #分組，這里是兩組，每組5個樣本 group <- rep(c('c1','c2'),each = 5)

1. 構建DGEList對象

根據基因表達量矩陣以及樣本分組信息構建DGEList對象

dgelist <- DGEList(counts = targets, group = group)

2. 過濾低表達的基因

DESeq2能夠自動識別這些低表達量的基因的，所以使用DESeq2時無需手動過濾。
edgeR推薦根據CPM（count-per-million）值進行過濾，即原始reads count除以總reads數乘以1,000,000,使用此類計算方式時，如果不同樣品之間存在某些基因的表達值極高或者極低，由于它們對細胞中分子總數的影響較大（也就是公式中的分母較大), 有可能導致標準化之后這些基因不存在表達差異，而原本沒有差異的基因在標準化之后卻顯示出差異

這里參考這篇：當我們在說RNA-seq reads count標準化時，其實在說什么？
為了解決上述問題，BSM(between-sample normalization)類分出control set去評估測序深度而不是用所有數據，主要分三種：
(1)?TMM(trimmed mean of M-values): TMM是M-值的加權截尾均值，即選定一個樣品為參照，其它樣品中基因的表達相對于參照樣品中對應基因表達倍數的log2值定義為M-值。隨后去除M-值中最高和最低的30%，剩下的M值計算加權平均值，權重來自Binomial data的delta方法 (Robinson and Oshlack, 2010)。
(2)?RLE (relative log expression):首先計算每個基因在所有樣品中表達的幾何平均值。然后再計算該值與每個樣品的比值的中位數，也叫被稱為量化因子scale factor (Anders and Huber 2010)。
(3)?UQ (upper quartile): 上四分位數 (upper quartile, UQ)是樣品中所有基因的表達除以處于上四分位數的基因的表達值。同時為了保證表達水平的相對穩定，計算得到的上四分位數值要除以所有樣品中上四分位數值的中位數。
以上三種方法效果大同小異，通常比較流行的是TMM和DESeq normalization

CPM 按照基因或轉錄本長度歸一化后的表達即 RPKM (Reads Counts Per Million)、FPKM (Fragments Per Kilobase Million)和 TPM (Trans Per Million)，推薦使用TPM######

1)直接選某個值過濾

keep <- rowSums(dgelist$counts) >= 50 dgelist <- dgelist[keep, ,keep.lib.sizes = FALSE]

2)利用cpm過濾

keep <- rowSums(cpm(dgelist) > 1 ) >= 2 dgelist <- dgelist[keep, ,keep.lib.sizes = FALSE]

實際的數據分析中，還需多加嘗試，選擇一個合適的過濾條件

3. 標準化

calcNormFactors()函數對數據標準化，以消除由于樣品制備或建庫測序過程中帶來的影響。
這里選的是TMM標準化方法，還有其他的可以?calcNormFactors進行查看

TMM法

?

dgelist_norm <- calcNormFactors(dgelist, method = 'TMM')

dgelist_norm

plotMDS()是limma包中的方法，繪制MDS圖，使用無監督聚類方法展示出了樣品間的相似性（或差異）。可據此查看各樣本是否能夠很好地按照分組聚類，評估試驗效果，判別離群點，追蹤誤差的來源等。

plotMDS(dgelist_norm, col = rep(c('red', 'blue'), each = 3))

可以考慮一下出現較大偏差的原因

4. 估算離散值

負二項分布（negative binomial，NB）模型需要均值和離散值兩個參數。
edgeR中，分組矩陣使用model.matrix()獲得，并可以通過estimateDisp()估算離散值。

design <- model.matrix(~group) #構建分組矩陣 dge <- estimateDisp(dgelist_norm, design, robust = TRUE) #估算離散值plotBCV(dge) #作圖查看

design中用0和1表示是哪一組，比如第二列1表示屬于c2組

?需要注意，標識好0、1類型后，后面的差異分析將以分組1的基因表達量相較于分組0是上調還是下調為準進行統計。因此在本示例中，后續分析獲得的基因表達量上調/下調均為分組c2相較于分組c1而言的。實際的分析中，切記不要搞反了。(有時會出現兩組順序相反的情況，還沒找到方法怎么實現)

estimateDisp()實際上是個組合函數，可以一步得到多個計算結果，例如在上文我們使用分組矩陣design通過estimateDisp()估算了3個值，其實就是estimateGLMTagwiseDisp()、estimateGLMCommonDisp()和estimateGLMTrendedDisp()這3個結果的組合。如果不指定分組矩陣，則將得到estimateCommonDisp()和estimateTagwiseDisp()的結果組合。

一定要記住是誰較誰

5. 差異基因分析

前面都是準備工作，現在可以開始正式分析了！

1) 負二項式廣義對數線性模型（edgeR）

首先擬合負二項式廣義對數線性模型（negative binomial generalized log-linear model），獲取差異基因。這種方法大致可以這樣理解，如果某個基因的表達值偏離這個分布模型，那么該基因即為差異表達基因。

使用edgeR包中的函數glmFit()和glmLRT()實現，其中glmFit()用于將每個基因的read count值擬合到模型中，glmLRT()用于對給定系數進行統計檢驗。

fit <- glmFit(dge, design, robust = TRUE) #擬合模型 lrt <- glmLRT(fit) #統計檢驗 topTags(lrt) #查看前10行 -n可修改查看前幾行

topTags(lrt)

write.csv(topTags(lrt, n = nrow(dgelist$counts)), 'glmLRT.csv', quote = FALSE) #輸出主要結果 dge_de <- decideTestsDGE(lrt, adjust.method = 'fdr', p.value = 0.05) #查看默認方法獲得的差異基因 summary(dge_de) plotMD(lrt, status = dge_de, values = c(1, -1), col = c('blue', 'red')) #作圖觀測 abline(h = c(-1, 1), col = 'gray', lty = 2) #在圖后加輔助線

decideTestsDGE() 的結果

decideTestsDGE()可用于統計差異基因數量。屏幕輸出了其默認值（供參考，大多數情況下我們還是優先根據Fold Change值以及p值等手動去篩選，而不會在意這個程序自己判斷的數值）
-1表示下調基因數量，1表示上調基因數量，0表示無差異基因數量。注意，對于這里的示例數據，基因表達量上調/下調均為“分組c2”相較于“分組c1”而言的。

輸出的 glmLRT.csv

logFC即log2轉化后的 Fold Change值，但是要注意，這里不是簡單的將基因的read count值直接對比，而是分別計算了基因在兩組中的CPM值，然后據此計算的logFC
logCPM是log2轉化后的CPM值
LR，似然比統計
PValue，差異表達的p值
FDR，FDR校正后的p值

最后結果，也可以畫火山圖

?

2) 類似然負二項式廣義對數線性模型（edgeR）

對于類似然負二項式廣義對數線性模型（quasi-likelihood negative binomial generalized log-linear model），使用edgeR包中的函數glmQLFit()和glmQLFTest()實現，同樣地，glmQLFit()用于將每個基因的read count值擬合到模型中，glmQLFTest()用于對給定系數進行統計檢驗，如果某個基因的表達值偏離這個分布模型，那么該基因即為差異表達基因。
相較于上一個模型，作者提到這個方法更嚴格一些。當然實際分析中還得視情況考慮了。

fit <- glmQLFit(dge, design, robust = TRUE) #擬合模型 lrt <- glmQLFTest(fit) #統計檢驗 topTags(lrt) #查看默認前10行

topTags(lrt)

write.csv(topTags(lrt, n = nrow(dgelist$counts)), 'glmQLFTest.csv', quote = FALSE) #輸出主要結果 dge_de <- decideTestsDGE(lrt, adjust.method = 'fdr', p.value = 0.05) #查看默認方法獲得的差異基因 summary(dge_de)

summary(dge_de)

plotMD(lrt, status = dge_de, values = c(1, -1), col = c('blue', 'red')) #作圖觀測 abline(h = c(-1, 1), col = 'gray', lty = 2)

跟第一種方法只有細微差別，大部分都是一樣的

3） 配對檢驗（edgeR）

除了擬合模型的方法外，在edgeR中還可使用exactTest()直接執行兩組負二項分布count之間基因均值差異的精確檢驗。

dge_et <- exactTest(dge) #檢驗 topTags(dge_et) write.csv(topTags(dge_et, n = nrow(dgelist$counts)), 'exactTest.csv', quote = FALSE) #輸出主要結果

topTags(dge_et)

dge_de <- decideTestsDGE(dge_et, adjust.method = 'fdr', p.value = 0.05) #查看默認方法獲得的差異基因 summary(dge_de)

summary(dge_de)

detags <- rownames(dge)[as.logical(dge_de)] plotSmear(dge_et, de.tags = detags, cex = 0.5) #作圖觀測 abline(h = c(-1, 1), col = 'gray', lty = 2)

因limma包的plotMD()函數無法在此處適用，這里使用的作圖函數plotSmear()是edgeR包中的方法
圖中縱軸為log2 Fold Change值；橫軸為log2 CPM值，反映了基因表達量信息；紅色的點表示差異基因（未使用顏色進一步區分上調/下調），黑色的點為無差異基因。

結果是這樣

?

4) voom線性建模（limma）

limma包可以說是處理RNA-seq數據上的“老大”了，功能強大自然無需多說。因此也很容易得知，limma包中同樣提供了多種差異基因分析的方法，其中最常用的就是voom方法（請允許我這么稱呼它）
我們仍可以基于上文前幾步獲得的預處理結果（DGEList對象、標準化數據、估算的離散值等），繼續使用limma包voom方法來完成后續的差異基因分析

將read count數據轉換為log2-counts per million（logCPM），通過估計均值-方差（mean-variance）關系并使用它來計算合適的observation-level weights，然后，數據就可以進行線性建模。好吧具體它怎么工作的咱也看不懂（voom參考文獻來源）……不過搞懂它的分析流程，以及結果怎么解讀，還是可以的

limma_voom <- voom(dgelist_norm, design, plot = TRUE)

limma_voom

fit <- lmFit(limma_voom, design) #擬合 fit <- eBayes(fit) topTable(fit, coef = 2) write.csv(topTable(fit, coef = 2, number = nrow(dgelist$counts)), 'limma_voom.csv', quote = FALSE) #輸出主要結果

topTable(fit, coef = 2)

---

💥💥💥二、DESeq2的簡潔使用💥💥💥

參考這篇

很慢，可以下這個

devtools::install_github('mikelove/DESeq2@ae7c6bd')

如果跟已安裝的包沖突的話，

remova.packages('xxx') BiocManager::install('xxx')

開始：

library(DESeq) x <- as.matrix(read.delim("你的路徑/gene.txt", sep = '\t', header = T, row.names = 1))

分組，這里是兩組，每組5個樣本

group <- rep(c('c1','c2'),each = 5)

由于DESeq包要求接下來的count data必須要整數型，因此我們需要對數據進行取整,然后將數據x和分組信息group讀入到cds對象中

database <- round(as.matrix(x)) cds <- newCountDataSet(database,group)

有生物學重復

cds <- estimateSizeFactors(cds) cds <- estimateDispersions(cds) res <- nbinomTest(cds,"c1","c2")

部分有生物學重復，其實同上

cds <- estimateSizeFactors(cds) cds <- estimateDispersions(cds) res <- nbinomTest(cds,"c1","c2")

沒有生物學重復

cds <- estimateDispersions(cds, method="blind", sharingMode="fit-only" ) res <- nbinomTest(cds,"c1","c2")

查看符合閾值的基因

table(res$padj <0.05) res <- res[order(res$padj),] sum(res$padj<=0.01,na.rm = T) write.csv(res,"路徑")

res.csv結果是這樣的

---

💥💥💥三、DESeq2的詳細使用💥💥💥

參考這篇: DESeq2做差異分析

0. 一些前期準備：

“gene.txt”，是一個基因表達量數據矩陣，包含10列樣本，10個樣本中前5個樣本屬于control組（c），后5個樣本屬于treat組（t）

0.1 構建基因表達矩陣countdata，即讀入數據?read.delim()

?colData的行名要與countData的列名一致！！

gene <- read.delim('C:/Users/wang/Desktop/gene.txt', row.names = 1, sep = '\t', stringsAsFactors = FALSE, check.names = FALSE)

剔除全為0值的行

all <- apply(gene, 1, function(x) all(x==0) ) newdata <- gene[!all,]

指定分組因子順序：差異基因分析需要指定比較分組的先后順序，以確定誰相對于誰的表達量上調/或下調。
···第一種方式是在讀取分組文件后，將分組列轉換為因子類型（factor），并指定因子（分組）順序，因子順序指定對照在前處理在后；
···第二種方式是在后續使用results()獲取差異結果時，指定比較的分組（推薦這種）
注意要保證表達矩陣中的樣本順序和這里的分組順序是一一對應的，前5列為一組，后5列為一組

0.2 構建colData,

?colData的行名要與countData的列名一致！！

colData <- data.frame(group = factor(rep(c('control', 'treat'), each = 5))) colData <- data.frame(row.names=colnames(gene), colData)

兩者的內容，參考這篇(https://www.jianshu.com/p/3a0e1e3e41d0)

1. 構建?DESeqDataSet對象，標準化reads count值，并用于存儲輸入值、中間計算和差異分析的結果

1.1 構建 DESeqDataSet 對象 dds = DESeqDataSet Object

①預處理，將所有樣本基因表達量之和小于1的基因過濾掉（這步？）

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

②差異分析

dds <- DESeqDataSetFromMatrix(countData = gene, colData = colData, design = ~group)

1.2 查看歸一化后的 count 值分布

boxplot(log10(assays(dds)[['cooks']]), range = 0, las = 2) plotDispEsts(dds)

?

cooks距離，詳見http://www.bioconductor.org/packages/release/bioc/vignettes/DESeq2/inst/doc/DESeq2.html

?

但是報錯了

?

我看了一下這個顯示NULL

第二天先運行了1.3的內容后，再運行這里就可以了，不明原因 :-(

boxplot()結果

?

plotDispEsts(dds)

?

1.3?vst標準化，獲取歸一化的基因表達矩陣norm_matrix,?blind = FALSE指定實驗設計不直接用于轉換

vsd <- assay(vst(dds, blind = FALSE)) head(vsd, 10) write.table(vsd, 'norm_matrix.txt', sep = '\t', col.names = NA, quote = FALSE)

vsd

2. 差異基因分析

之后直接運行默認的DESeq2差異分析流程就可以了
函數DESeq()是一個包含因子大小估計、離散度估計、負二項模型擬合、Wald統計等多步在內的過程，結果將返回至DESeqDataSet對象。這步比較耗時，特別是數據量較大時，新舊版DESeq2的運算效率差距極為明顯
通過result()可獲得最終計算的log2倍數變化和校正后p值等信息
①contrast參數用于指定比較的分組順序，即誰相對于誰的表達量上調/或下調
②pAdjustMethod設定p值校正方法
③alpha為顯著性水平，這里0.05為校正后p值小于0.05即為顯著

2.1 標準方法

dds <- DESeq(dds, parallel = FALSE) #若 parallel = TRUE 將啟用多線程模式 suppressMessages(dds) res <- results(dds, contrast = c('group', 'treat', 'control'), pAdjustMethod = 'fdr', alpha = 0.05) write.csv(res, "你的路徑/res.csv")summary(res)plotMA(res, alpha = 0.05, ylim = c(-3, 3))

dds過程

?

suppressMessages(dds)

通過summary()，可以根據預先設定的校正后p值<0.05水平（alpha=0.05，由results()指定），輸出所比較兩組間的上調/下調基因數量。這個結果可供參考，在后續也可以自己根據log2FC和校正后p值自定義作篩選

summary()

?

plotMA()

?

到這兒我發現我和例子中的結果有些差別了，但是還沒找到原因，先過完流程吧 :-(

2.2 an alternate analysis: likelihood ratio test 似然比檢驗

ddsLRT <- DESeq(dds, test = 'LRT', reduced = ~ 1) suppressMessages(ddsLRT) resLRT <- results(ddsLRT, contrast = c('group', 'treat', 'control'), pAdjustMethod = 'fdr', alpha = 0.05) write.csv(resLRT, "你的路徑/ .csv")summary(resLRT)plotMA(resLRT, alpha = 0.05, ylim = c(-3, 3))

差異分析結果保存在res中，可通過as.data.frame()直接轉化為數據框類型
包含了基因id、標準化后的基因表達值的平均值baseMean、log2FoldChange值、顯著性p值pvalue以及校正后p值padj等主要信息

res結果

?

可以先大概看一些差異基因的數目

table(res$padj<0.05)

table()

2.3 可以先按校正和 p 值由小到大排個序，方便查看

deseq_res <- as.data.frame(res[order(res$padj), ])

將行名寫在gene_id列中，這個時候它是最后一列

deseq_res$gene_id <- rownames(deseq_res)

先輸出第7列，再輸出前6列

write.table(deseq_res[c(7, 1:6)], '你的路徑/DESeq2-test.txt', row.names = FALSE, sep = '\t', quote = FALSE)

最后的結果

3 ggplot2對差異基因作圖

3.1 讀進來最后的差異基因結果并進行分類

library(ggplot2) deseq_res <- read.delim('你的路徑 / DESeq2-test.txt', sep = '\t')

|log2FC| >= 1 & FDR p-value < 0.05?定為差異

deseq_res[which(deseq_res$padj %in% NA),'sig'] <- 'no diff' deseq_res[which(deseq_res$log2FoldChange >= 1 & deseq_res$padj < 0.05),'sig'] <- 'up (p.adj < 0.05, log2FC >= 1)' deseq_res[which(deseq_res$log2FoldChange <= -1 & deseq_res$padj < 0.05),'sig'] <- 'down (p.adj < 0.05, log2FC <= -1)' deseq_res[which(abs(deseq_res$log2FoldChange) < 1 | deseq_res$padj >= 0.05),'sig'] <- 'no diff'

也可以獲取上調up /下調down 的差異表達基因（padjust < 0.05，并且|log2(foldchange)|>1)

diff_up = subset(deseq_res,padj < 0.05 & (log2FoldChange > 1)) write.csv(diff_up,file="diff_up.csv",row.names = F)diff_down = subset(deseq_res,padj < 0.05 & (log2FoldChange > 1)) write.csv(diff_down,file="diff_down.csv",row.names = F)

3.2 畫火山圖

①縱軸為-log10(pvalue)，橫坐標為log2FoldChange，差異基因展示為不同顏色

volcano_p <- ggplot(deseq_res, aes(log2FoldChange, -log(padj, 10))) +geom_point(aes(color = sig), alpha = 0.6, size = 1) +scale_color_manual(values = c('blue2', 'gray30', 'red2')) +theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.position = c(0.26, 0.92)) +theme(legend.title = element_blank(), legend.key = element_rect(fill = 'transparent'), legend.background = element_rect(fill = 'transparent')) +geom_vline(xintercept = c(-1, 1), color = 'gray', size = 0.25) +geom_hline(yintercept = -log(0.05, 10), color = 'gray', size = 0.25) +labs(x = 'log2 Fold Change', y = '-log10 p-value', color = NA) +xlim(-5, 5)volcano_p ggsave('你的路徑/volcano_p.png', volcano_p, width = 5, height = 6)

sig映射到color中
背景中fill = 'transparent'，使背景變為透明色
geom_vline和geom_hline為在x軸和y軸添加輔助線
labs在x軸和y軸添加橫縱坐標名稱
xlim限定x軸的顯示范圍
ggsave保存圖片，或者可以直接Export

火山圖

?

②縱軸為log2FoldChange，橫坐標展示為標準化后的基因表達量的平均值 Average log10 baseMean，差異基因用不同顏色顯示

volcano_count <- ggplot(deseq_res, aes(y = log2FoldChange, x = log10(baseMean))) +geom_point(aes(color = sig), alpha = 0.6, size = 1) +scale_color_manual(values = c('blue2', 'gray30', 'red2')) +theme(panel.grid = element_blank(), panel.background = element_rect(color = 'black', fill = 'transparent'), legend.position = c(0.2, 0.9)) +theme(legend.title = element_blank(), legend.key = element_rect(fill = 'transparent'), legend.background = element_rect(fill = 'transparent')) +geom_hline(yintercept = c(-1, 1), color = 'gray', size = 0.25) +labs(y = 'log2 Fold Change', x = 'Average log10 baseMean') +ylim(-5, 5)volcano_countggsave('你的路徑/volcano_count.png', volcano_p, width = 5, height = 6)

火山圖

4 用biomaRt注釋基因

參考這篇

4.1 我們利用useMart()函數選擇“ENSEMBL_MART_ENSEMBL”，并將其賦值給my_mart對象

library('biomaRt') library("curl")my_mart <-useMart("ensembl")

在ensembl數據庫中包含了77個數據集，可用下面這樣的方式查看

datasets <- listDatasets(my_mart) View(datasets)

datasets

4.2 選擇一個數據集datasset，這里選人類的

my_dataset <- useDataset("hsapiens_gene_ensembl",mart = my_mart)

4.3 💥根據ensembl ID獲取基因名、描述或染色體信息

💥💥💥這里前半部分有誤！請一定往下看解決辦法

my_newid <- getBM(attributes = c("ensembl_gene_id","external_gene_name","description","chromosome_name"),filters = "ensembl_gene_id",values = newinput,mart = my_dataset)

?

image.png

這里一直報錯，并且輸出的為內容為0行
找到原因是：EBI數據庫沒有小數點，所以需要進一步替換為整數的形式。需要把小數點去掉！！這個很重要，所以需要加一個步驟

?

①還是將差異文件的行名提取出來

inputdata <- as.data.frame(row.names(deseq_res))

②這里將匹配到的.以及后面的數字連續匹配并替換為空，并重新賦值，一定要是data.frame格式

newinput <- as.data.frame(gsub("\\.\\d*", "", inputdata[,1]))

③getBM()轉換ID

1）attributes參數：用來指定輸出的數據類型，就是你要什么，比如entrezgene，hgnc_id。忘記的話可以用listAttributes(你自定義的dataset)
2）filters參數：用來指定數據的輸入類型，比如你的原始信息是基因的ensembl ID，并且有這些基因的染色體位置信息，那么此處的filter就是ensembl_gene_id和chromosome_name等。
3）values參數：就是你待轉換ID的數據
4）mart參數：此前定義的數據庫，此處就是my_dataset

那么在我這里：
attributes?：我想要輸出"ensembl_gene_id",轉換后的"external_gene_name",轉換后的"description",還可以有"chromosome_name"
filters：我的原始數據"ensembl_gene_id"
mart：之前建立的數據庫

listAttributes(你的dataset)?可以查看可供選擇的attributes

listAttributes(my_dataset)

my_result <- getBM(attributes = c("ensembl_gene_id","external_gene_name","description"),filters = "ensembl_gene_id",values = newinput,mart = my_dataset)

?

ID轉換成功后

這樣就完成了對ensembl_id的轉化和注釋

?

4.4 最后需要把結果文件deseq_res和注釋文件my_result兩者merge起來

merge需要有相同的gene_id
💥但是一定要看看自己文件里的gene_id是不是一致，如果有一個為小數，就要再添加一列取整后的gene_id

①?deseq_res中gene_id有小數點 所以再加一列變成new_deseq_res，新增加的列名為gene_new_id

new_deseq_res <- as.data.frame(deseq_res) new_deseq_res$gene_new_id <- gsub("\\.\\d*", "", deseq_res$gene_id)

② 修改一下列名,把含有小數點的列命名為gene_all_id,取整后的為gene_id，這一步是為了方便merge

colnames(new_deseq_res) <- c('baseMean', 'log2FoldChange','lfcSE','stat','pvalue','padj','gene_all_id','gene_id')

new_deseq_res

③?merge兩個文件，即new_deseq_res和my_resullt，生成final_res文件

by = intersect(names(x), names(y))?為取兩個文件所有列名中列名相同的那列！

final_res <- merge(my_result, new_deseq_res, by = intersect(names(my_result), names(new_deseq_res))) write.table(final_res, 'C:/Users/wang/Desktop/final_result.txt',row.names = FALSE, sep = '\t', quote = FALSE)

結果文件

4.5 還可以找到某個基因所在的通路GO號

參考這篇

① 選出要查找的基因

#舉個例子 entrez = c("673", "837")

② 利用ensembl構建my_mart和my_dataset

my_mart <-useMart("ensembl") #`listDatasets()`可以查看可用的`datasets` datasets <- listDatasets(my_mart) View(datasets)#構建`my_dataset` my_dataset <- useDataset("hsapiens_gene_ensembl",mart = my_mart)

③ 查看可輸出的attributes

listAttributes(my_dataset)

④ 查找GOid

GOid <- getBM(attributes = c('entrezgene_id', 'go_id'),filters = 'entrezgene_id',values = entrez,mart = my_dataset)

結果

4.6 與4.5相反，可以通過所在的通路GO號找到某個基因

getBM(attributes = c('entrezgene_id', 'ensembl_gene_id'),filters = 'go',values = 'GO:0005524',mart = my_dataset)

?

總結

以上是生活随笔為你收集整理的edgeR/limma/DESeq2差异基因分析→ggplot2作火山图→biomaRt转换ID并注释的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。