R语言基因表达数据转化log r语言基因id转symbol

转载

fjfdh 2023-10-09 13:17:40

文章标签 R语言基因表达数据转化log hg bc 数据集 文章分类 R语言后端开发

1.请根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)
包中自带函数toTable可以将各种命名方式转换为数据框
其中每种命名方式都和共同的gene_id对应,可以通过gene_id对各个命名数据框进行联结操作。

> head(toTable(org.Hs.egENSEMBL))
  gene_id      ensembl_id
1       1 ENSG00000121410
2       2 ENSG00000175899
3       3 ENSG00000256069

options(stringAsfactor = F)
id <- read.table(col.names = c("ensembl_id_all"),
                   text = "ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11
")

library(org.Hs.eg.db)
library(stringr)
id$ensembl_id <- str_split(id$ensembl_id_all,"[.]",simplify = T)[,1]
g2s <- toTable(org.Hs.egSYMBOL)
g2e <- toTable(org.Hs.egENSEMBL)
id <- merge(id,g2e,by="ensembl_id",all.x = T)
id <- merge(id,g2s,by="gene_id",all.x = T)
id <- id[order(id$ensembl_id_all),]  # 排序

for (i in 1:nrow(id)) {
  print(paste(id$ensembl_id_all[i],'->',id$symbol[i]))
}

[1] "ENSG00000000003.13 -> TSPAN6"
[1] "ENSG00000000005.5 -> TNMD"
[1] "ENSG00000000419.11 -> DPM1"
[1] "ENSG00000000457.12 -> SCYL3"
[1] "ENSG00000000460.15 -> C1orf112"
[1] "ENSG00000000938.11 -> FGR"

自己的方法

# 读入数据，用col.names指定列名
> id <- read.table(col.names = c("id"),text = "
+                  ENSG00000000003.13
+ ENSG00000000005.5
+ ENSG00000000419.11
+ ENSG00000000457.12
+ ENSG00000000460.15
+ ENSG00000000938.11
+                  ")
> name_list <- select(org.Hs.eg.db,keys = str_sub(as.vector(id$id), 1, 15),columns = c("SYMBOL","ENTREZID","GENENAME"),keytype = "ENSEMBL")
'select()' returned 1:1 mapping between keys and columns
> cat(str_c("\n",name_list$ENSEMBL, " -> ", name_list$SYMBOL))

ENSG00000000003 -> TSPAN6 
ENSG00000000005 -> TNMD 
ENSG00000000419 -> DPM1 
ENSG00000000457 -> SCYL3 
ENSG00000000460 -> C1orf112 
ENSG00000000938 -> FGR

查看所有的命名方式

> keytypes(org.Hs.eg.db)
 [1] "ACCNUM"       "ALIAS"        "ENSEMBL"      "ENSEMBLPROT"  "ENSEMBLTRANS" "ENTREZID"     "ENZYME"      
 [8] "EVIDENCE"     "EVIDENCEALL"  "GENENAME"     "GO"           "GOALL"        "IPI"          "MAP"         
[15] "OMIM"         "ONTOLOGY"     "ONTOLOGYALL"  "PATH"         "PFAM"         "PMID"         "PROSITE"     
[22] "REFSEQ"       "SYMBOL"       "UCSCKG"       "UNIGENE"      "UNIPROT"

2.根据R包hgu133a.db找到下面探针对应的基因名(symbol)

library(hgu133a.db)
Pid <- read.table(col.names = c("probe_id"),
                  text = "1053_at
117_at
                  121_at
                  1255_g_at
                  1316_at
                  1320_at
                  1405_i_at
                  1431_at
                  1438_at
                  1487_at
                  1494_f_at
                  1598_g_at
                  160020_at
                  1729_at
                  177_at
                  ")

Pid <- merge(Pid,toTable(hgu133aSYMBOL),by = "probe_id", x.all = T)

3.找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图

suppressPackageStartupMessages(library(CLL)) # 加载没有信息提示
data(sCLLex) # 载入CLL数据集

#获取表达和样本信息
exprSet=exprs(sCLLex) 
pd <- pData(sCLLex)

#加载注释包
library(hgu95av2.db)
ids <- toTable(hgu95av2SYMBOL)

# 找到基因对应探针
tp53choose <- ids[,2]=="TP53"
TP53_probe_id<- ids[tp53choose,][1]

# 画图
> boxplot(exprSet[TP53_probe_id[1,1],]~pd$Disease)
> boxplot(exprSet["1939_at",]~pd$Disease)
> boxplot(exprSet[TP53_probe_id$probe_id[1],]~pd$Disease)

R语言基因表达数据转化log r语言基因id转symbol_hg

6.下载数据集GSE17215的表达矩阵并且提取下面的基因画热图

rm(list = ls()) 
options(stringsAsFactors = F)

# 使用GEOquery包下载GEO数据（需提前下载安装此包）
# 注意下载文件大小是否正确
# 感觉中科大的biocondater好用
library(GEOquery)
f='GSE17215_eSet.Rdata'
if(!file.exists(f)){
  gset <- getGEO('GSE17215', destdir=".",
                 AnnotGPL = F,     ## 注释文件
                 getGPL = F)       ## 平台文件
  save(gset,file=f)   ## 保存到本地
}
# 先保存后加载
load('GSE17215_eSet.Rdata')
# 因为这个GEO数据集只有一个GPL平台，所以下载到的是一个含有一个元素的list            
class(gset)
length(gset)
class(gset[[1]])#ExpressionSet
dim(gset[[1]])
head(exprs(gset[[1]]))


# not list易于操作
a=gset[[1]]
# 获取表达矩阵
dat=exprs(a)
dim(dat)
# 加载注释包
library(hgu133a.db)
# 看一下包里有什么
ls("package:hgu133a.db")#查下包里有哪些函数，这里为了查看有哪些基因命名方式
# 获取探针与基因的对应关系
ids=toTable(hgu133aSYMBOL)
head(ids)
# 筛选ids中包含的probe_id，将两个表行名设为同序
dat=dat[ids$probe_id,]
dat[1:4,1:4] 
# 取dat表达量的中位数，生成ids的median列
ids$median=apply(dat,1,median)
# 因为只保留最大的，将ids按照symbol和median排序（降序），为了去除重复探针时，保留最大检测量的探针，一个基因可能对应多种探针。
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
# 去除重复的symbol值，保留第一位，即表达量最大的一位（dim一下看看前后对比）
ids=ids[!duplicated(ids$symbol),]
# 同上
dat=dat[ids$probe_id,]
# 先同序，再直接改为symbol（dim一下瞧瞧两者行数）
rownames(dat)=ids$symbol
dat[1:4,1:4]  
dim(dat)
# 好，进入正题
ng='ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T'
# 把基因加上“”，使之与rownames(dat)格式一致
ng=strsplit(ng,' ')[[1]]
#ng=str_split(ng,' ')[[1]]
# 看一下，多少基因存在于dat中，可得41个TRUE&9个FALSE
table(ng %in%  rownames(dat))
# 清洗掉不存在的ng，注意这一步存在排序（连同下一步理解）
ng=ng[ng %in%  rownames(dat)]
# 再改行名
dat=dat[ng,]
# 画图
dat=log2(dat)
pheatmap::pheatmap(dat,scale = 'row')

R语言基因表达数据转化log r语言基因id转symbol_数据集_02

自己的做法

# 下载数据
f = 'GSE17215_eSet.Rdata'
options("BioC_mirror"<- "https://mirrors.ustc.edu.cn/bioc/")
library(GEOquery)
if (!exists(f)) {
  gset <- getGEO('GSE17215',
                 AnnotGPL = F,
                 getGPL = F)
  save(gset,file = f)
}
load(f)
#class(gset)
#length(gset)
#class(gset[[1]])   ExpressionSet
eSet <- gset[[1]]
dat <- exprs(eSet)
# 转换探针数据
library(hgu133a.db)
#ls("package:hgu133a.db")#看下包中函数，同时也可得到命名方式
ids <- toTable(hgu133aSYMBOL)
dat <- dat[ids$probe_id,]
ids$median <- apply(dat,1,median)
ids <- ids[order(ids$symbol,ids$median,decreasing = T),]
ids <- ids[!duplicated(ids$symbol),]
dat <- dat[ids$probe_id,]
rownames(dat) <- ids$symbol

#热图
gn <- "ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T"
library(stringr)
gn <- str_split(gn,' ')[[1]]
table(gn %in% rownames(dat))
gn <- gn[gn %in% rownames(dat)]
dat <- dat[gn,]
library(pheatmap)
#dat <- log2(dat + 1)
#pheatmap(dat,scale = "row")

dat2 <- t(scale(t(dat)))
hc=hclust(dist(t(dat2)))
clus = cutree(hc, 2)
group_list <- as.factor(clus)
ac=data.frame(g=group_list)
pheatmap(dat2,annotation_col = ac)

> ac
          g
GSM431121 1
GSM431122 1
GSM431123 2
GSM431124 2
GSM431125 2
GSM431126 2

R语言基因表达数据转化log r语言基因id转symbol_数据集_03

7.下载数据集GSE24673的表达矩阵计算样本的相关性并且绘制热图，需要标记上样本分组信息解

rm(list = ls())  
options(stringsAsFactors = F)
library(GEOquery)
GSE_name = 'GSE24673'
options( 'download.file.method.GEOquery' = 'libcurl' ) #windows系统
gset <- getGEO( GSE_name, getGPL = F)

# 保存并载入数据
save( gset, file = 'GSE24673_gset.Rdata' )
load(file ="GSE24673_gset.Rdata")
a <- gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# 我们通过观察样本信息（pd）来得知11个样本的性质（分组信息），建立小分组
group_list=c('rbc','rbc','rbc',
             'rbn','rbn','rbn',
             'rbc','rbc','rbc',
             'normal','normal')
dat[1:4,1:4]
# 相关性分析并作图
M=cor(dat)
pheatmap::pheatmap(M)
# 标记分组信息，相关性分析并作图
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)

R语言基因表达数据转化log r语言基因id转symbol_R语言基因表达数据转化log_04

自己方法

options()$BioC_mirror 
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
rm(list = ls())  
options(stringsAsFactors = F)
library(GEOquery)
f='GSE24673_gset.Rdata'
if (!exists(f)) {
  gset <- getGEO('GSE24673',
                 AnnotGPL = F,
                 getGPL = F)
  save(gset,file = f)
}
load(f)
a <- gset[[1]]
dat <- exprs(a)
pd <- pData(a)

# 我们通过观察样本信息（pd）来得知11个样本的性质（分组信息），建立小分组
group_list=c('rb139','rb139','rb139',
             'rb535','rb535','rb535',
             'rb984','rb984','rb984',
             'normal','normal')
# 用列的mean代替缺失值
dat2=apply(dat,2,function(x){
  x[is.na(x)]=mean(x,na.rm = T)
  return(x)
})

# 相关性分析并作图
M=cor(dat2)

tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)

R语言基因表达数据转化log r语言基因id转symbol_R语言基因表达数据转化log_05

8.找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 对应的R的bioconductor注释包，并且安装它！

参考：http://www.bio-info-trainee.com/1399.html

#BiocManager::install("请输入自己找到的R包",ask = F,update = F)
if (!require(hugene10sttranscriptcluster.db)) {
  BiocManager::install("hugene10sttranscriptcluster.db")
}
library(hugene10sttranscriptcluster.db)

9下载数据集GSE42872的表达矩阵，并且分别挑选出所有样本的(平均表达量/sd/mad/)最大的探针，并且找到它们对应的基因。

rm(list = ls())
options(stringsAsFactors = F)
options()$BioC_mirror 
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
f='GSE42872_eSet.Rdata'
library(GEOquery)
if (!exists(f)) {
  gset <- getGEO('GSE42872',
                 AnnotGPL = F,
                 getGPL = F)
  save(gset,file = f)
}
load(f)
a <- gset[[1]]
dat <- exprs(a)
pd <- pData(a)

boxplot(dat,col = 2:7,las = 2)

mea <- sort(apply(dat,1,mean),decreasing = T)[1]
s <- sort(apply(dat,1,sd),decreasing = T)[1]
ma <- sort(apply(dat,1,mad),decreasing = T)[1]
names(mea)

library('hugene10sttranscriptcluster.db')
ls('package:hugene10sttranscriptcluster.db')
ids <- toTable(hugene10sttranscriptclusterSYMBOL)
head(ids)
ids[ids$probe_id == names(s),]

10.下载数据集GSE42872的表达矩阵，并且根据分组使用limma做差异分析，得到差异结果矩阵

rm(list = ls()) 
options(stringsAsFactors = F)
options()$BioC_mirror 
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
gse <- 'GSE42872'
f <- paste0(gse,'_eSet.Rdata')
library(GEOquery)
if (!file.exists(f)) {
  gset <- getGEO(gse,
                 AnnotGPL = F,
                 getGPL = F)
  save(gset,file = f)
}
load(f)

#length(gset)
a <- gset[[1]]
dat <- exprs(a)
pd <- pData(a)
boxplot(dat,col = 2:(ncol(dat) + 1), las = 2, cex.axis =0.7)

#Step2: 构建design 样本分组矩阵

group <- factor(rep(c("Control","Vemurafenib"),each = 3))
design <- model.matrix(~0+group)
rownames(design) = colnames(dat)
colnames(design) <- levels(group)


library(limma)
fit <- lmFit(dat,design)

contrast.matrix<-makeContrasts(paste0(unique(group),collapse = "-"),levels = design)
contrast.matrix 
contrast.matrix <- makeContrasts(Control-Vemurafenib,levels = design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
DEG <- topTable(fit2,coef = 1,n=Inf)
nrDEG <- na.omit(DEG)

#Visualization
p_thred <- 0.05
logFC_thred <- 2
p4draw <- function(DEG,p_thred,logFC_thred) {
  pp <- data.frame(symbols=rownames(DEG),
                   logFC=DEG$logFC,
                   P.Value=DEG$P.Value)
  pp$groups = ifelse(pp$P.Value > p_thred, "stable", 
                     ifelse(pp$logFC > logFC_thred, "up", 
                            ifelse(pp$logFC < -logFC_thred, "down", 
                                   "stable")))
  return(pp)
}
pp <- p4draw(nrDEG,p_thred,logFC_thred)
table(pp$groups)
#  down stable     up 
#    92  33144     61 

library(ggplot2)
drawVolcano <- function(pp){
  g = ggplot(data = pp, aes(x = logFC, y = -log10(P.Value), color = groups)) + 
    geom_point(alpha = 0.4, size = 1.75) + 
    geom_vline(xintercept = c(-2,2),lty=3,lwd=0.8) +
    geom_hline(yintercept=-log10(0.05),lty=3,lwd=0.8) +
    xlim(c(-7,7)) +
    theme_set(theme_set(theme_bw(base_size = 20))) + 
    labs(x="log2 fold change",y="-log10 p-value",title="GSE42872_DEG_Volcano") +
    theme(plot.title = element_text(size = 15,hjust = 0.5)) +
    scale_colour_manual(values = c("blue","black", "red")) +
    theme(legend.title = element_blank())
  return(g)
}
drawVolcano(pp)

R语言基因表达数据转化log r语言基因id转symbol_R语言基因表达数据转化log_06

本文章为转载内容，我们尊重原作者对文章享有的著作权。如有内容错误或侵权问题，欢迎原作者联系我们进行内容更正或删除文章。

上一篇：spring boot使用自定义配置项限制系统的总注册用户数 springboot自定义自动配置

下一篇：python 挂单 python下单脚本

提问和评论都可以，用心的回复会被更多人看到评论

发布评论

相关文章

官方博客	全部文章	热门标签	班级博客
了解我们	网站地图	意见反馈

鸿蒙开发者社区	51CTO学堂
51CTO	软考资讯