1. clusterProfiler包，bitr()转换ID函数
clusterProfiler 是Y叔写的一个功能强大的R包，可以用来做各种富集分析，如GO、KEGG、DO（Disease Ontology analysis）、Reactome pathway analysis、GSEA富集分析等。还具有非常优秀的富集分析结果可视化功能。
1.1 bitr()的使用方法：
bitr(geneID, fromType, toType, OrgDb, drop = TRUE)
geneID：一个含有基因名的向量
orgDb：人类的注释包是org.Hs.eg.db，小鼠是org.Mm.eg.db
fromType：输入的基因名的类型
toType：需要转换成的类型，可以是多种类型，用大写character类型向量表示
1.2 查看org.Hs.eg.db中可以被选择/使用的类型
keytypes(org.Hs.eg.db)
[1] "ACCNUM" "ALIAS" "ENSEMBL" "ENSEMBLPROT" "ENSEMBLTRANS"
[6] "ENTREZID" "ENZYME" "EVIDENCE" "EVIDENCEALL" "GENENAME"
[11] "GO" "GOALL" "IPI" "MAP" "OMIM"
[16] "ONTOLOGY" "ONTOLOGYALL" "PATH" "PFAM" "PMID"
[21] "PROSITE" "REFSEQ" "SYMBOL" "UCSCKG" "UNIGENE"
[26] "UNIPROT"
1.3 示例
library(clusterProfiler)
library(org.Hs.eg.db)
s2e <- bitr(deg$symbol,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)
deg <- inner_join(deg,s2e,by=c("symbol"="SYMBOL")) #将转换后的结果加入矩阵，变成一列
需要注意的是，一些数据比如从TCGA上下载的数据，行名是带有版本号的，如下 ：
View(exp)
a = rownames(exp)
head(a)
# [1] "ENSG00000000003.13"
# [2] "ENSG00000000419.11"
# [3] "ENSG00000000457.12"
# [4] "ENSG00000000460.15"
# [5] "ENSG00000000938.11"
# [6] "ENSG00000000971.14"
这样的ENTREZID是不能被bitr()识别的，需要按点号分隔开，保留点号前面的部分，才能被bitr()识别。
library(clusterProfiler)
library(org.Hs.eg.db)
library(stringr)
a = str_split(a,"\\.",simplify = T)[,1] #不带\\直接写"."是正则表达式任意字符的意思，需要\\来转义。
id = bitr(a,
fromType = "ENSEMBL",
toType = "SYMBOL",
OrgDb = "org.Hs.eg.db")
2. 探针ID转化为gene symbol
2.1 获得探针和gene symbol的对应关系
- 方法一：根据gpl_number，从http://www.bio-info-trainee.com/1399.html中寻找对应的bioc_package，下载对应的包。
比如最常用的GPL570，对应hgu133plus2
if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db") #注意加上.db
library(hgu133plus2.db)
ls("package:hgu133plus2.db")
# [1] "hgu133plus2"
# [2] "hgu133plus2_dbconn"
# [3] "hgu133plus2_dbfile"
# [4] "hgu133plus2_dbInfo"
# [5] "hgu133plus2_dbschema"
# [6] "hgu133plus2.db"
# [7] "hgu133plus2ACCNUM"
# [8] "hgu133plus2ALIAS2PROBE"
# [9] "hgu133plus2CHR"
# [10] "hgu133plus2CHRLENGTHS"
# [11] "hgu133plus2CHRLOC"
# [12] "hgu133plus2CHRLOCEND"
# [13] "hgu133plus2ENSEMBL"
# [14] "hgu133plus2ENSEMBL2PROBE"
# [15] "hgu133plus2ENTREZID"
# [16] "hgu133plus2ENZYME"
# [17] "hgu133plus2ENZYME2PROBE"
# [18] "hgu133plus2GENENAME"
# [19] "hgu133plus2GO"
# [20] "hgu133plus2GO2ALLPROBES"
# [21] "hgu133plus2GO2PROBE"
# [22] "hgu133plus2MAP"
# [23] "hgu133plus2MAPCOUNTS"
# [24] "hgu133plus2OMIM"
# [25] "hgu133plus2ORGANISM"
# [26] "hgu133plus2ORGPKG"
# [27] "hgu133plus2PATH"
# [28] "hgu133plus2PATH2PROBE"
# [29] "hgu133plus2PFAM"
# [30] "hgu133plus2PMID"
# [31] "hgu133plus2PMID2PROBE"
# [32] "hgu133plus2PROSITE"
# [33] "hgu133plus2REFSEQ"
# [34] "hgu133plus2SYMBOL"
# [35] "hgu133plus2UNIGENE"
# [36] "hgu133plus2UNIPROT"
ids <- toTable(hgu133plus2SYMBOL)
head(ids)
# probe_id symbol
# 1 1053_at RFC2
# 2 117_at HSPA6
# 3 121_at PAX8
# 4 1255_g_at GUCA1A
# 5 1316_at THRA
# 6 1320_at PTPN21
- 方法二：读取GPL平台的soft文件，按列取子集
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570
a = getGEO(gpl_number,destdir = ".")
b = a@dataTable@table
colnames(b)
ids2 = b[,c("ID","Gene Symbol")]
colnames(ids2) = c("probe_id","symbol")
ids2 = ids2[ids2$symbol!="" & !str_detect(ids2$symbol,"///"),]
- 方法三：官网下载，文件读取
参考：http://www.affymetrix.com/support/technical/byproduct.affx?product=hg-u133-plus - 方法四：自主注释
参考：https://mp.weixin.qq.com/s/mrtjpN8yDKUdCSvSUuUwcA
2.2 加上探针注释
#为deg数据框添加几列（deg是做完差异分析后，行为probe_id，列为logFC等差异分析结果的矩阵）
#1.加probe_id列，把行名变成一列
library(dplyr)
deg <- mutate(deg,probe_id=rownames(deg))
head(deg)
#2.加上探针注释
#多个探针对应一个基因：按照基因去重复
table(!duplicated(ids$probe_id))
table(!duplicated(ids$symbol)) #结果显示很多基因对应了多个探针
#按symbol列去重，常见标准有3个：最大值/平均值/随机去重
#随机去重，另两个见zz.去重方式.R
ids = ids[!duplicated(ids$symbol),] #保留重复值中第一个出现的
deg <- inner_join(deg,ids,by="probe_id")
head(deg)
nrow(deg)
3. gtf文件获得id转换信息
对于TCGA上下载的数据，其使用的geneid是ensembl id，做完差异分析后需要把结果中的每个ensembl id转换为对应的symbol和类型(mRNA/lncRNA或其它)。而gtf文件中就包含了我们需要的gene symbol和ensembl id的对应关系。
思路：
1. 找到TCGA数据对应的参考基因组注释版本。
2. 下载该版本的参考基因组注释文件，提取ensembl id 与symbol的对应关系及每个基因的gene type信息。
3. 可以将symbol和gene type 用merge添加到差异分析结果中，也可以在差异分析前先转换矩阵的行名。
3.1 找参考基因组版本
可以看到，使用的参考基因组版本是genecode的v22。（版本很多，这个是14年的版本了）
在gtf文件里并不是直接分出了lncRNA，需要找gtf文件里对biotype的说明，不看不知道，一看发现这是一个很长的表格。
进入https://www.gencodegenes.org，点击Documentation下的Biotypes，
其中对lncRNA的说明是：Generic long non-coding RNA biotype that replaced the following biotypes: 3prime_overlapping_ncRNA, antisense, bidirectional_promoter_lncRNA, lincRNA, macro_lncRNA, non_coding, processed_transcript, sense_intronic and sense_overlapping.
所以需要将genetype里这些类型对应的行挑出来，就是lncRNA了。 然后与表达矩阵行名进行匹配替换，就可以分别得到mRNA和lncRNA的矩阵了。
3.2 转换
#step1:读取并探索gtf文件----
#BiocManager::install("rtracklayer")
library(rtracklayer)
gtf = rtracklayer::import("gencode.v22.annotation.gtf")
class(gtf)
gtf = as.data.frame(gtf);dim(gtf)
colnames(gtf)
table(gtf$type)
#step2:先筛选出gene对应的行
gtf_gene = gtf[gtf$type=="gene",]
save(gtf_gene,file = "gtf_gene.Rdata")
table(rownames(deg) %in% gtf_gene$gene_id) #deg是需要进行id转换的矩阵
# FALSE TRUE
# 3 30345
an = gtf_gene[,c("gene_name","gene_id","gene_type")] #取出gtf文件中这三列
head(an)
# gene_name gene_id gene_type
# 1 DDX11L1 ENSG00000223972.5 transcribed_unprocessed_pseudogene
# 13 WASH7P ENSG00000227232.5 unprocessed_pseudogene
# 26 MIR6859-3 ENSG00000278267.1 miRNA
# 29 RP11-34P13.3 ENSG00000243485.3 lincRNA
# 37 MIR1302-9 ENSG00000274890.1 miRNA
# 40 FAM138A ENSG00000237613.2 lincRNA
deg = merge(deg,an,by.x = "row.names",by.y = "gene_id")
# mRNA和lncRNA总共有多少个？
lnc = c("3prime_overlapping_ncRNA", "antisense", "bidirectional_promoter_lncRNA", "lincRNA", "macro_lncRNA", "non_coding", "processed_transcript", "sense_intronic" , "sense_overlapping")
k1 = gtf_gene$gene_type %in% lnc;table(k1) #gtf中共14826个lnc
# k1
# FALSE TRUE
# 45657 14826
k2 = gtf_gene$gene_type == "protein_coding";table(k2) #gtf中共19814个mRNA
# k2
# FALSE TRUE
# 40669 19814
# deg中有多少mRNA和lncRNA?
k3 = deg$gene_type %in% lnc;table(k3) #deg中共7501个lnc
# k3
# FALSE TRUE
# 22844 7501
k4 = deg$gene_type =="protein_coding";table(k4) #deg中共17464个mRNA
# k4
# FALSE TRUE
# 12881 17464
# 差异的mRNA和lncRNA 各有多少
k5 = deg$change !="NOT"
table(k3&k5)
# FALSE TRUE #有396个差异性lnc
# 29949 396
table(k4&k5)
# FALSE TRUE #有1084个差异性mRNA
# 29261 1084
表达矩阵的行名id转换
load("gtf_gene.Rdata")
an = gtf_gene[,c("gene_name","gene_id","gene_type")]
exp = exp[rownames(exp) %in% an$gene_id,] #从表达矩阵中去除不存在于gtf中的gene名
an = an[match(rownames(exp),an$gene_id),] #将an的顺序调整成和表达矩阵行名一致
identical(an$gene_id,rownames(exp)) #检查是否一致
# [1] TRUE
k = !duplicated(an$gene_name);table(k) #因为要用gene_symbol做行名，而行名是不能重复的，所以要先去重
# k
# FALSE TRUE
# 193 30152
an = an[k,]
exp = exp[k,]
rownames(exp) = an$gene_name
常用的转换需求总结：
- GO分析：需要将gene symbol转化为ENTREZID
- 分析芯片数据：需要将探针ID转化为gene symbol
- 将TCGA上下载的矩阵中的ENTREZID转化为gene symbol（1和3实质上是一样的），可以直接转换矩阵，也可以在做完差异分析后，把gene symbol作为一列连接在差异表达矩阵上。
代码来自2021生信技能树数据挖掘课
网友评论