從我寫的RNA-seq摸索：4. edgeR/limma/DESeq2差異基因分析→ggplot2作火山圖→biomaRt轉(zhuǎn)換ID并注釋中提取出來這部分，方便以后修改補(bǔ)充

參考這篇

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

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

my_mart <-useMart("ensembl")

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

datasets <- listDatasets(my_mart)
View(datasets)

datasets

2 選擇一個(gè)數(shù)據(jù)集datasset，這里選人類的

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

3 ??根據(jù)ensembl ID獲取基因名、描述或染色體信息

??????這里前半部分有誤！請(qǐng)一定往下看解決辦法

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ù)庫(kù)，此處就是my_dataset

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

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 最后需要把結(jié)果文件deseq_res和注釋文件my_result兩者merge起來

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

① deseq_res中gene_id有小數(shù)點(diǎn) 所以再加一列變成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ù)點(diǎn)的列命名為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兩個(gè)文件，即new_deseq_res和my_resullt，生成final_res文件

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

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é)果文件

5 還可以找到某個(gè)基因所在的通路GO號(hào)

參考這篇

① 選出要查找的基因

#舉個(gè)例子
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é)果

6 與5相反，可以通過所在的通路GO號(hào)找到某個(gè)基因

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