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

三者區(qū)別

三者區(qū)別

差異分析流程:

1 初始數(shù)據(jù)

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

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

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

4 多重假設(shè)得FDR值：

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

5 差異基因篩選：pvalue、padj

??????一龄砰、 edgeR的使用??????

因為目前沒有合適的數(shù)據(jù)盟猖，所以數(shù)據(jù)來源于這里 參考這篇：劉堯科學網(wǎng)博客

0. 前期工作

1. 用到的gene.txt文件，內(nèi)容如下
   
   gene.txt文件內(nèi)容
   
   c1表示為一組换棚，c2表示為另一組式镐，.后為第幾個樣本
   讀取數(shù)據(jù)并設(shè)置分組，要保證樣本名稱和分組名稱的順序是一一對應的固蚤。

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

1. 構(gòu)建DGEList對象

根據(jù)基因表達量矩陣以及樣本分組信息構(gòu)建DGEList對象

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

2. 過濾低表達的基因

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

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

CPM 按照基因或轉(zhuǎn)錄本長度歸一化后的表達即 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]

實際的數(shù)據(jù)分析中踢匣，還需多加嘗試告匠，選擇一個合適的過濾條件

3. 標準化

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

TMM法

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

dgelist_norm

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

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

可以考慮一下出現(xiàn)較大偏差的原因

4. 估算離散值

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

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

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

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

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

一定要記住是誰較誰

5. 差異基因分析

前面都是準備工作，現(xiàn)在可以開始正式分析了悠鞍！

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

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

使用edgeR包中的函數(shù)glmFit()和glmLRT()實現(xiàn)么翰，其中glmFit()用于將每個基因的read count值擬合到模型中牺汤，glmLRT()用于對給定系數(shù)進行統(tǒng)計檢驗。

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

topTags(lrt)

write.csv(topTags(lrt, n = nrow(dgelist$counts)), 'glmLRT.csv', quote = FALSE) #輸出主要結(jié)果
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() 的結(jié)果

decideTestsDGE()可用于統(tǒng)計差異基因數(shù)量硬鞍。屏幕輸出了其默認值（供參考慧瘤，大多數(shù)情況下我們還是優(yōu)先根據(jù)Fold Change值以及p值等手動去篩選，而不會在意這個程序自己判斷的數(shù)值）
-1表示下調(diào)基因數(shù)量固该，1表示上調(diào)基因數(shù)量锅减，0表示無差異基因數(shù)量。注意伐坏，對于這里的示例數(shù)據(jù)怔匣，基因表達量上調(diào)/下調(diào)均為“分組c2”相較于“分組c1”而言的。

輸出的 glmLRT.csv

logFC即log2轉(zhuǎn)化后的 Fold Change值桦沉，但是要注意每瞒，這里不是簡單的將基因的read count值直接對比，而是分別計算了基因在兩組中的CPM值纯露，然后據(jù)此計算的logFC
logCPM是log2轉(zhuǎn)化后的CPM值
LR剿骨，似然比統(tǒng)計
PValue，差異表達的p值
FDR埠褪，F(xiàn)DR校正后的p值

最后結(jié)果浓利，也可以畫火山圖

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

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

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

topTags(lrt)

write.csv(topTags(lrt, n = nrow(dgelist$counts)), 'glmQLFTest.csv', quote = FALSE)  #輸出主要結(jié)果
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()直接執(zhí)行兩組負二項分布count之間基因均值差異的精確檢驗汇在。

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

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()函數(shù)無法在此處適用，這里使用的作圖函數(shù)plotSmear()是edgeR包中的方法
圖中縱軸為log2 Fold Change值脏答；橫軸為log2 CPM值糕殉，反映了基因表達量信息；紅色的點表示差異基因（未使用顏色進一步區(qū)分上調(diào)/下調(diào)）殖告，黑色的點為無差異基因阿蝶。

結(jié)果是這樣

4) voom線性建模（limma）

limma包可以說是處理RNA-seq數(shù)據(jù)上的“老大”了，功能強大自然無需多說黄绩。因此也很容易得知羡洁，limma包中同樣提供了多種差異基因分析的方法，其中最常用的就是voom方法（請允許我這么稱呼它）
我們?nèi)钥梢曰谏衔那皫撞将@得的預處理結(jié)果（DGEList對象爽丹、標準化數(shù)據(jù)筑煮、估算的離散值等），繼續(xù)使用limma包voom方法來完成后續(xù)的差異基因分析

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

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) #輸出主要結(jié)果

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必須要整數(shù)型，因此我們需要對數(shù)據(jù)進行取整,然后將數(shù)據(jù)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結(jié)果是這樣的

??????三贿条、DESeq2的詳細使用??????

參考這篇: DESeq2做差異分析

0. 一些前期準備：

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

0.1 構(gòu)建基因表達矩陣countdata钓葫，即讀入數(shù)據(jù) 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,]

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

0.2 構(gòu)建colData,

?colData的行名要與countData的列名一致Ｑ纪ⅰ枣抱！

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

兩者的內(nèi)容，參考這篇(http://www.reibang.com/p/3a0e1e3e41d0)

1. 構(gòu)建 DESeqDataSet對象辆床，標準化reads count值佳晶，并用于存儲輸入值、中間計算和差異分析的結(jié)果

1.1 構(gòu)建 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

boxplot()結(jié)果

plotDispEsts(dds)

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

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差異分析流程就可以了
函數(shù)DESeq()是一個包含因子大小估計吱型、離散度估計逸贾、負二項模型擬合、Wald統(tǒng)計等多步在內(nèi)的過程津滞，結(jié)果將返回至DESeqDataSet對象铝侵。這步比較耗時，特別是數(shù)據(jù)量較大時触徐，新舊版DESeq2的運算效率差距極為明顯
通過result()可獲得最終計算的log2倍數(shù)變化和校正后p值等信息
①contrast參數(shù)用于指定比較的分組順序咪鲜，即誰相對于誰的表達量上調(diào)/或下調(diào)
②pAdjustMethod設(shè)定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()

plotMA()

到這兒我發(fā)現(xiàn)我和例子中的結(jié)果有些差別了览祖，但是還沒找到原因，先過完流程吧 :-(

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))

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

res結(jié)果

可以先大概看一些差異基因的數(shù)目

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)

最后的結(jié)果

3 ggplot2對差異基因作圖

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

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'

也可以獲取上調(diào)up /下調(diào)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'丝里，使背景變?yōu)橥该魃?br> 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_count

ggsave('你的路徑/volcano_count.png', volcano_p, width = 5, height = 6)

火山圖

4 用biomaRt注釋基因

參考這篇

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

library('biomaRt')
library("curl")

my_mart <-useMart("ensembl")

在ensembl數(shù)據(jù)庫中包含了77個數(shù)據(jù)集营密，可用下面這樣的方式查看

datasets <- listDatasets(my_mart)
View(datasets)

datasets

4.2 選擇一個數(shù)據(jù)集datasset械媒，這里選人類的

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

4.3 ??根據(jù)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

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

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

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

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

③getBM()轉(zhuǎn)換ID

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

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

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轉(zhuǎn)換成功后

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

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

① deseq_res中gene_id有小數(shù)點 所以再加一列變成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)

② 修改一下列名,把含有小數(shù)點的列命名為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)

結(jié)果文件

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

參考這篇

① 選出要查找的基因

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

② 利用ensembl構(gòu)建my_mart和my_dataset

my_mart <-useMart("ensembl") 

#`listDatasets()`可以查看可用的`datasets`
datasets <- listDatasets(my_mart)
View(datasets)

#構(gòu)建`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)

結(jié)果

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

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