生信中级10题

根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)

genes <- c("ENSG00000000003.13","ENSG00000000005.5","ENSG00000000419.11","ENSG00000000457.12","ENSG00000000460.15","ENSG00000000938.11")

library(org.Hs.eg.db)
options(stringsAsFactors = F)
ls("package:org.Hs.eg.db")
ids <- toTable(org.Hs.egENSEMBL2EG)
#等同于ids <- toTable(org.Hs.egENSEMBL)
ids2 <- toTable(org.Hs.egSYMBOL)
str(ids) #查看提取出的数据框变量的特征
str(ids2)
ids$gene_id <- as.numeric(ids$gene_id)  #避免合并后gene_id出现错误
ids2$gene_id <- as.numeric(ids2$gene_id)#避免合并后gene_id出现错误
ids_t <- merge(ids,ids2,by="gene_id")

library(stringr)
genes <- str_match(genes,"([:alnum:]*)[.][:digit:]")[,2]
match_results <- ids_t[ids_t$ensembl_id%in%genes,]
match_results

rm(list=ls())  #一键清除

match_results

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

probe_ids <- c('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")

library(hgu133a.db)
ls("package:hgu133a.db")
ids <- toTable(hgu133aSYMBOL)
symbols_res <- ids[ids$probe_id%in%probe_ids,]
symbols_res

rm(list=ls())  #一键清除

symbols_res

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

• 方法一：

library(CLL)
?CLL  #查看帮助文档
data(sCLLex)
class(sCLLex)  #查看sCLLex是什么对象，决定了怎么提取其中的数据
eSet <- sCLLex
rm(sCLLex)
options(stringsAsFactors = F)
exprSet <- exprs(eSet)
pdata <- pData(eSet)
str(pdata)  #查看数据框变量特征后知道pdata$Disease为因子，但stable不是对照水平
group_list <- relevel(pdata$Disease,ref = "stable")  #用relevel()将stable设定为对照水平

gpl <- eSet@annotation  #用@来取出eSet对象中的注释数据
if (!require(hgu95av2.db)) BiocManager::install("hgu95av2.db")
library(hgu95av2.db)
ls("package:hgu95av2.db")
ids <- toTable(hgu95av2SYMBOL)

expr_df <- as.data.frame(exprSet) #将表达矩阵转为数据框
expr_df$median <- apply(expr_df,1,median)  #将数据框添加median一列
expr_df$probe_id <- rownames(expr_df) #添加probe_id一列
expr_df <- merge(expr_df,ids,by="probe_id")  #将expr_df与ids按照probe_id进行合并

expr_df <- expr_df[order(expr_df$symbol,expr_df$median,decreasing = T),]  #将symbol排序后，再将median从大到小排序
expr_df <- expr_df[!duplicated(expr_df$symbol),]  #根据symbol去重复

rownames(expr_df) <- expr_df$symbol
expr_df <- expr_df[,-c(1,ncol(expr_df)-1,ncol(expr_df))]

identical(colnames(expr_df),rownames(pdata))  #确定表达矩阵中的列名和pdata中的行名是否一致

ex_TP53 <- expr_df[rownames(expr_df)== "TP53",]

library(ggpubr)
library(reshape2)
ex_TP53 <- melt(ex_TP53,variable.name = "variable",value = "value")
str(ex_TP53)
df <- ex_TP53
identical(df$variable,colnames(expr_df))
df$group <- group_list
names(df)[2] <- "expression"
p <- ggboxplot(df, x = "group", y = "expression",
               color = "group", palette =c("#00AFBB", "#E7B800", "#FC4E07"),
               add = "jitter", shape = "group")+
  stat_compare_means(label.y = 4.0,method = "t.test")  
p

ggpubr结果

• 方法二：

suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex  #直接查看数据集（S4对象）中的信息
exprSet=exprs(sCLLex)   
##sCLLex是依赖于CLL这个package的一个对象
samples=sampleNames(sCLLex)
samples  #查看samples
pdata=pData(sCLLex)
pdata
group_list=as.character(pdata[,2]) #pdata[,2]提取后为因子，需转为字符串
dim(exprSet)  #查看基因表达矩阵的维度，即行数及列数（探针数及样本数）
exprSet[1:5,1:5]  #当基因集较大时，查看部分数据可减少内存消耗

#BiocManager::install('hgu95av2.db')  #安装注释包
library(hgu95av2.db)
ls("package:hgu95av2.db") #查看注释包的内容

ids=toTable(hgu95av2SYMBOL)  #toTable可将Bimap对象转为数据框结构
save(ids,exprSet,pdata,file = 'input.Rdata')
length(unique(ids$symbol))  #查看共有多少个symbol，一般有多个探针对应一个symbol
tail(sort(table(ids$symbol)))  ##查看基因对应最多多少个探针
table(sort(table(ids$symbol)))  #显示总的symbol对应的探针数量并table

plot(table(sort(table(ids$symbol))))

table(rownames(exprSet) %in% ids$probe_id)  #注意两个数据集的先后顺序，%in%表示前者中的变量与后者一一比较，返回布尔值
dim(exprSet)
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,]  #取得与注释包交集后的表达矩阵
dim(exprSet)

ids=ids[match(rownames(exprSet),ids$probe_id),]  #match相当于%in%，可以匹配返回逻辑值，也可以将前者数据的顺序匹配给后者
head(ids)
exprSet[1:5,1:5]

identical(ids$probe_id,rownames(exprSet))  #identical判断两个对象是否一致
dat=exprSet
ids$median=apply(dat,1,median) #ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#先对ids$symbol排序，在此基础上按照ids$median中位数从大到小排列的顺序排序，将对应的行赋值为一个新的ids
ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
dat[1:4,1:4]  #保留每个基因ID第一次出现的信息
dim(dat)

df <- dat[rownames(dat)== "TP53",]
identical(colnames(dat),rownames(pdata))
df <- data.frame(expression = df,group = group_list)

library(ggpubr)
p <- ggboxplot(df, x = "group", y = "expression",
               color = "group", palette =c("#00AFBB", "#E7B800", "#FC4E07"),
               add = "jitter", shape = "group")+
  stat_compare_means(label.y = 4.0,method = "t.test") +
  ylab("relative expression of TP53")
p

与第一种方法相比，第二种方法的组别应转为因子，并设定对照stable

• 法三：使用dplyr包

library(dplyr)
data(sCLLex)
class(sCLLex)  #查看sCLLex是什么对象，决定了怎么提取其中的数据
eSet <- sCLLex
rm(sCLLex)
exprSet <- exprs(eSet)
pdata <- pData(eSet)
str(pdata)  #查看数据框变量特征后知道pdata$Disease为因子，但stable不是对照水平
group_list <- relevel(pdata$Disease,ref = "stable")  #用relevel()将stable设定为对照水平

gpl <- eSet@annotation  #用@来取出eSet对象中的注释数据
if (!require(hgu95av2.db)) BiocManager::install("hgu95av2.db")
library(hgu95av2.db)
ls("package:hgu95av2.db")
ids <- toTable(hgu95av2SYMBOL)

dat <- exprSet%>%
  as.data.frame()%>%
  mutate(probe_id=rownames(exprSet),median = apply(exprSet,1,median))
identical(dat$probe_id,rownames(exprSet))
dat <- merge(dat,ids,by = "probe_id")
dat <- dat[order(dat$symbol,dat$median,decreasing= "T"),]
table(sort(table(dat$symbol)))
dat <- dat[!duplicated(dat$symbol),]
rownames(dat) <- dat$symbol
dat <- dat[,-c(1,ncol(dat)-1,ncol(dat))]

identical(colnames(dat),rownames(pdata))
df <- as.numeric(dat[rownames(dat)== "TP53",])
df <- data.frame(expression = df,group = group_list)

library(ggpubr)
p <- ggboxplot(df, x = "group", y = "expression",
               color = "group", palette =c("#00AFBB", "#E7B800", "#FC4E07"),
               add = "jitter", shape = "group")+
  stat_compare_means(label.y = 4.0,method = "t.test") +
  ylab("relative expression of TP53")
p

第三种方法ggpubr作图

找到BRCA1基因在TCGA数据库的乳腺癌数据集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表达情况

因网络问题，暂无解答。提示：使用http://www.cbioportal.org/index.do 定位数据集：http://www.cbioportal.org/datasets

找到TP53基因在TCGA数据库的乳腺癌数据集的表达量分组看其是否影响生存

提示使用：http://www.oncolnc.org/ 数据库

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

genes <- c("ACTR3B", "ANLN","BAG1","mBCL2","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")
rm(list=ls())
library(GEOquery)
library(stringr)
library(AnnoProbe)
library(ggplot2)

gse = "GSE17215"

if(require(AnnoProbe)){
  if(!file.exists(paste0(gse,"_eSet.Rdata"))) geoChina(gse)
  load(paste0(gse,"_eSet.Rdata"))
  eSet <- gset
  rm(gset)
}else{
  eSet <- getGEO(gse, 
                 destdir = '.', 
                 getGPL = F)
}

#(1)提取表达矩阵exp
exp <- exprs(eSet[[1]])
exp[1:4,1:4]
exp = log2(exp+1)
boxplot(exp,cex.axis = 0.75,las = 2,col = rainbow(6))
#(2)提取临床信息
pd <- pData(eSet[[1]])
#(3)提取芯片平台编号
gpl <- eSet[[1]]@annotation

p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]

#(4)分组信息：
group_list <- ifelse(str_detect(pd$title,"paclitaxel"),"paclitaxel","salinomycin")

group_list=factor(group_list,levels = c("paclitaxel","salinomycin"))

#(5)下载注释信息：
if (!require(hthgu133a.db)) BiocManager::install("hthgu133a.db")
library(hthgu133a.db)
ls("package:hthgu133a.db")
ids <- toTable(hthgu133aSYMBOL)

#(6)数据处理：
table(rownames(exp)%in%ids$probe_id)
exp <- exp[rownames(exp)%in%ids$probe_id,]
identical(ids$probe_id,rownames(exp))

exp <- as.data.frame(exp)
exp$probe_id <- rownames(exp)
ids$median <- apply(exp,1,median)
ids <- ids[order(ids$symbol,ids$median,decreasing = T),]
ids <- ids[!duplicated(ids$symbol),]

exp <- merge(exp,ids,by="probe_id")
exp <- exp[match(ids$probe_id,exp$probe_id),]
identical(ids$probe_id,exp$probe_id)
rownames(exp) <- exp$symbol
exp <- exp[,-c(1,ncol(exp)-1,ncol(exp))]

dat <- exp[rownames(exp)%in%genes,]
dat <- as.matrix(dat)

library(RColorBrewer)
library(pheatmap)

identical(rownames(pd),colnames(dat))
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(dat) 

pheatmap(dat,scale = "row",
         cluster_cols = T,
         cluster_rows = T,
         show_rownames = T,
         show_colnames = F,
         color = colorRampPalette(rev(brewer.pal(n = 7, name =
                                                   "RdYlBu")))(100),
         annotation_col = annotation_col,
         fontsize = 8, #调整图片中文本的大小
         legend = T)

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

rm(list=ls())
library(GEOquery)
library(stringr)
library(AnnoProbe)
library(ggplot2)

gse = "GSE24673"

if(require(AnnoProbe)){
  if(!file.exists(paste0(gse,"_eSet.Rdata"))) geoChina(gse)
  load(paste0(gse,"_eSet.Rdata"))
  eSet <- gset
  rm(gset)
}else{
  eSet <- getGEO(gse, 
                 destdir = '.', 
                 getGPL = F)
}

#(1)提取表达矩阵exp
exp <- exprs(eSet[[1]])
exp[1:4,1:4]
#exp = log2(exp+1)
ncol(exp)
boxplot(exp,cex.axis = 0.75,las = 2,col = rainbow(11))
#(2)提取临床信息
pd <- pData(eSet[[1]])
#(3)提取芯片平台编号
gpl <- eSet[[1]]@annotation

p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]

#(4)分组信息：
group_list <- c(rep(c("139","535","984"),each = 3),rep("normal",2))

group_list=factor(group_list,levels = c("normal","139","535","984"))

#(5)下载注释信息：
if (!require(hugene10sttranscriptcluster.db)) BiocManager::install("hugene10sttranscriptcluster.db")
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")
ids <- toTable(hugene10sttranscriptclusterSYMBOL)

#(6)数据处理：
table(rownames(exp)%in%ids$probe_id)
exp <- exp[rownames(exp)%in%ids$probe_id,]
ids <- ids[ids$probe_id%in%rownames(exp),]
identical(ids$probe_id,rownames(exp))

ids$median <- apply(exp,1,median)
ids <- ids[order(ids$symbol,ids$median,decreasing = T),]
ids <- ids[!duplicated(ids$symbol),]

exp <- as.data.frame(exp)
exp$probe_id <- rownames(exp)
exp <- merge(exp,ids,by="probe_id")
exp <- exp[match(ids$probe_id,exp$probe_id),]
identical(ids$probe_id,exp$probe_id)
rownames(exp) <- exp$symbol
exp <- exp[,-c(1,ncol(exp)-1,ncol(exp))]

#画PCA图
dat=as.data.frame(t(exp))
library(FactoMineR)#画主成分分析图需要加载这两个包
library(factoextra) 

#pca的统一操作走起
dat.pca <- PCA(dat, graph = FALSE)
pca_plot <- fviz_pca_ind(dat.pca,
                         geom.ind = "point", # show points only (nbut not "text")
                         col.ind = group_list, # color by groups
                         #palette = c("#00AFBB", "#E7B800"),
                         addEllipses = TRUE, # Concentration ellipses
                         legend.title = "Groups"
)
pca_plot  #浓度椭圆

pca_plot <- fviz_pca_ind(dat.pca,
                         geom.ind = "point", # show points only (nbut not "text")
                         col.ind = group_list, # color by groups
                         #palette = c("#00AFBB", "#E7B800"),
                         addEllipses = TRUE, 
                         ellipse.type = "confidence",
                         legend.title = "Groups"
)
pca_plot  #置信椭圆
#参见：http://www.sthda.com/english/articles/31-principal-component-methods-in-r-practical-guide/112-pca-principal-component-analysis-essentials/

#表达矩阵热图：
library(RColorBrewer)
library(pheatmap)

dat=as.data.frame(t(dat))
identical(rownames(pd),colnames(dat))
annotation_col=data.frame(group=group_list)
rownames(annotation_col)=colnames(dat) 

pheatmap(dat,scale = "row",
         cluster_cols = T,
         cluster_rows = T,
         show_rownames = F,
         show_colnames = F,
         color = colorRampPalette(rev(brewer.pal(n = 7, name =
                                                   "RdYlBu")))(100),
         annotation_col = annotation_col,
         fontsize = 8,#调整图片中文本的大小
         legend = T)

#样本相关性热图：
M <- cor(exp)
pheatmap::pheatmap(M,annotation_col = annotation_col)  #(存疑)

PCA 表达矩阵热图 样本相关性矩阵

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

GPL6244

if (!require(hugene10sttranscriptcluster.db)) BiocManager::install("hugene10sttranscriptcluster.db")
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")
ids <- toTable(hugene10sttranscriptclusterSYMBOL)

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

（注释包：hugene10sttranscriptcluster.db）

rm(list=ls())
#数据下载
options(stringsAsFactors = F)
library(GEOquery)
gse = "GSE42872"
if(require(AnnoProbe)){
  if(!file.exists(paste0(gse,"_eSet.Rdata"))) geoChina(gse)
  load(paste0(gse,"_eSet.Rdata"))
  eSet <- gset
  rm(gset)
}else{
  eSet <- getGEO(gse, 
                 destdir = '.', 
                 getGPL = F)
}
#(1)提取表达矩阵exp
exp <- exprs(eSet[[1]])
exp[1:4,1:4]
#exp = log2(exp+1)
#(2)提取临床信息
pd <- pData(eSet[[1]])
#(3)调整pd的行名顺序与exp列名完全一致
p = identical(rownames(pd),colnames(exp));p
if(!p) exp = exp[,match(rownames(pd),colnames(exp))]
#(4)提取芯片平台编号
gpl <- eSet[[1]]@annotation
save(gse,pd,exp,gpl,file = paste0(gse,"step1output.Rdata"))

#group_list(实验分组)和ids(芯片注释)，每次都需要改
rm(list = ls())  
load(file = "GSE42872step1output.Rdata")
group_list=ifelse(str_detect(pd$title,"Control"),"control","treat")
#设置参考水平，对照在前，处理在后
group_list = factor(group_list,
                    levels = c("control","treat"))
#ids  BioconductorR包
gpl 
#http://www.bio-info-trainee.com/1399.html
if(!require(hugene10sttranscriptcluster.db))BiocManager::install("hugene10sttranscriptcluster.db")
library(hugene10sttranscriptcluster.db)
ls("package:hugene10sttranscriptcluster.db")
ids <- toTable(hugene10sttranscriptclusterSYMBOL)
head(ids)
save(exp,group_list,ids,file = paste0(gse,"step2output.Rdata"))

#数据处理：
table(rownames(exp)%in%ids$probe_id)
exp <- exp[rownames(exp)%in%ids$probe_id,]
ids <- ids[ids$probe_id%in%rownames(exp),]
identical(ids$probe_id,rownames(exp))

ids$median <- apply(exp,1,median) #中位数
ids$sd <- apply(exp,1,sd) #标准差
ids$mad <- apply(exp,1,mad) #中位方差
ids_median <- ids[order(ids$median,decreasing = F),]
median_max <- ids_median[1,]
ids_sd <- ids[order(ids$sd,decreasing = T),]
sd_max <- ids_sd[1,];sd_max
ids_mad <- ids[order(ids$mad,decreasing = T),]
mad_max <- ids_mad[1,];mad_max

所有样本的(平均表达量/sd/mad/)最大的探针及基因名

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

load(file = "GSE42872step1output.Rdata")
load("GSE42872step2output.Rdata")

table(rownames(exp)%in%ids$probe_id)
exp <- exp[rownames(exp)%in%ids$probe_id,]
identical(rownames(exp),ids$probe_id)
exp <- as.data.frame(exp)
exp$median <- apply(exp,1,median)
exp$probe_id <- rownames(exp)
exp_m <- merge(exp,ids,by="probe_id")

exp_m <- exp_m[order(exp_m$symbol,exp_m$median,decreasing = T),]
exp_m <- exp_m[!duplicated(exp_m$symbol),]
rownames(exp_m) <- exp_m$symbol
exp_m <- exp_m[,-c(1,ncol(exp_m)-1,ncol(exp_m))]

library(limma)
design=model.matrix(~group_list)
fit=lmFit(exp_m,design)
fit=eBayes(fit)
deg=topTable(fit,coef=2,number = Inf)
save(exp_m,deg,file="GSE42872limma.Rdata")