ccRCC数据分析-GSE66270-GPL570
数据集介绍
GEO链接：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66270
芯片平台：GPL570, [HG-U133_Plus_2] Affymetrix Human Genome U133 Plus 2.0 Array]
样品列表：
14 normal tissue samples and 14 tumor
GSM1618387 primary ccRCC, patient 1
GSM1618388 adj. normal kidney sample, patient 1
GSM1618389 primary ccRCC, patient 2
GSM1618390 adj. normal kidney sample, patient 2
GSM1618391 primary ccRCC, patient 3
GSM1618392 adj. normal kidney sample, patient 3
GSM1618393 primary ccRCC, patient 4
GSM1618394 adj. normal kidney sample, patient 4
GSM1618395 primary ccRCC, patient 5
GSM1618396 adj. normal kidney sample, patient 5
GSM1618397 primary ccRCC, patient 6
GSM1618398 adj. normal kidney sample, patient 6
GSM1618399 primary ccRCC, patient 7
GSM1618400 adj. normal kidney sample, patient 7
GSM1618401 primary ccRCC, patient 8
GSM1618402 adj. normal kidney sample, patient 8
GSM1618403 primary ccRCC, patient 9
GSM1618404 adj. normal kidney sample, patient 9
GSM1618405 primary ccRCC, patient 10
GSM1618406 adj. normal kidney sample, patient 10
GSM1618407 primary ccRCC, patient 11
GSM1618408 adj. normal kidney sample, patient 11
GSM1618409 primary ccRCC, patient 12
GSM1618410 adj. normal kidney sample, patient 12
GSM1618411 primary ccRCC, patient 13
GSM1618412 adj. normal kidney sample, patient 13
GSM1618413 primary ccRCC, patient 14
GSM1618414 adj. normal kidney sample, patient 14
文章链接是： Integrated microRNA and mRNA Signature Associated with the Transition from the Locally Confined to the Metastasized Clear Cell Renal Cell Carcinoma Exemplified by miR-146-5p. PLoS One 2016;11(2):e0148746. PMID: 26859141 https://pubmed.ncbi.nlm.nih.gov/26859141/ Cooperative Effect of miR-141-3p and miR-145-5p in the Regulation of Targets in Clear Cell Renal Cell Carcinoma. PLoS One 2016;11(6):e0157801. PMID: 27336447
核心步骤
R包加载
rm(list = ls())
library(AnnoProbe)
library(GEOquery)
library(ggplot2)
library(ggstatsplot)
library(reshape2)
library(patchwork)
library(stringr)
2. 获取并且检查表达量矩阵
本数据集的关键是：当发现表达矩阵有很多负值（经过归一化的矩阵）时，如何根据cel文件获得count矩阵
# 获取表达量矩阵
gse_number <- 'GSE66270'
gset <- geoChina(gse_number)
a=gset[[1]]
dat=exprs(a)
dim(dat)
# 检查，判断需不需要取log
dat[1:4,1:4]
boxplot(dat[,1:4],las=2)
#通过观察boxplot可以很容易地看出数据是经过归一化的矩阵（图非常地一目了然），而我们最好是拿到count自己分析，所以需要去下载raw_data（cel文件）到工作目录，然后自己处理，代码见下一个chunk。
在 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE66270下载芯片原始测序数据-cel文件到工作目录，解压得到download文件，然后运行如下代码：
rm(list = ls())
gse_number = "GSE66270"
gpl_number = "GPL570"
if(!require(oligo)) BiocManager::install(c( 'oligo' ),ask = F,update = F)
if(!require(pd.hg.u133.plus.2)) BiocManager::install(c( 'pd.hg.u133.plus.2' ),ask = F,update = F)
#目前数据已经下载并解压一次，得到工作目录下的 Download 文件夹
dir='./download'
od=getwd()
setwd(dir)
celFiles <- list.celfiles(listGzipped = T)
celFiles
affyRaw <- read.celfiles( celFiles )
setwd(od)
eset <- rma(affyRaw)
eset
save(eset,celFiles,file = "GSE66270_exp.Rdata")
write.exprs(eset,file="data.txt")
#接下来是常规流程
获取并且检查表达量矩阵
a=eset
dat=exprs(a) #a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
dat[1:4,1:4]
library(limma)
dat=normalizeBetweenArrays(dat)
dat[1:4,1:4]
#重命名样本名称
colnames_dat = substring(colnames(dat),1,18);colnames_dat
colnames(dat)=colnames_dat
dat[1:4,1:4]
# 画图，使用ggplot需宽数据变长数据
class(dat)
data <- as.data.frame(dat)
data <- melt(data)
head(data)
gse_number = "GSE66270"
gpl_number = "GPL570"
title <- paste (gse_number, "/", gpl_number, sep ="")
p1 <- ggplot( data,aes(x=variable,y=value))+
geom_boxplot()+
theme_ggstatsplot()+
theme(panel.grid = element_blank(),
axis.text=element_text(size=10,face = 'bold'),
axis.text.x=element_text(angle=90),
plot.title = element_text(hjust = 0.5,size =15))+
xlab('')+
ylab('')+
ggtitle(title)
p1
可以看到，处理前后我们的表达量矩阵的表达量范围箱线图如下所示：
根据生物学背景及研究目的人为分组
#生成样本分组
group_list_pair = substring(celFiles,12,18);group_list_pair
group_list = ifelse(str_detect(group_list_pair,"NN"),"control","ccRCC");group_list
group_list = factor(group_list,levels = c("control","ccRCC"))
table(group_list)
probe_id 和symbol的转换并对应至表达矩阵
获取芯片注释信息
代码如下：
gpl_number='GPL570'
ids = idmap(gpl_number)#极简函数，一步注释
可以看到此芯片的探针与基因ID或者symbol的对应关系如下所示:
> head(ids)
probe_id symbol
## 193731 1053_at RFC2
## 193732 117_at HSPA6
## 193733 121_at PAX8
## 193734 1255_g_at GUCA1A
## 193735 1316_at THRA
## 193736 1320_at PTPN21
探针基因ID对应以及去冗余
代码如下：
library(tidyr)
library(dplyr)
#接下来，使探针与基因symbol一一对应
ids=as.data.frame(ids)
table(rownames(dat) %in% ids$probe_id)
dat=dat[rownames(dat) %in% ids$probe_id,]
ids=ids[match(rownames(dat),ids$probe_id),]
ids$probe_id=as.character(ids$probe_id)
rownames(dat)=ids$probe_id
ids=ids[ids$probe_id %in% rownames(dat),]
dat=dat[ids$probe_id,]
#下一步是：基因symbol去冗余-按照表达量最大值筛选
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 ，保留每个基因最大表达量结果s
#获得去冗余之后的dat/exp
dat=dat[ids$probe_id,] #新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat
#把ids的symbol这一列中的每一行给dat作为dat的行名
rownames(dat)=ids$symbol
dat[1:4,1:4]
#保存为R数据文件：step1-output.Rdata
save(dat,exp,ids,group_list,group_list_pair,gpl_number,file = "step1-output.Rdata")
最后得到了表达量矩阵如下所示：
> dat[1:4,1:4] #保留每个基因ID第一次出现的信息
GSM1618387_J862_NC GSM1618388_J862_NN GSM1618389_J976_NC GSM1618390_J976_NN
ZZZ3 8.821700 7.918982 7.908951 8.045114
ZZEF1 7.709334 7.582121 7.699571 7.611648
ZYX 8.730555 8.150840 8.061261 7.763805
ZYG11B 8.965619 9.786744 9.041259 9.098561
以及最简单的2分组，如下所示：
> table(group_list)
group_list
ccRCC control
14 14
三、标准步骤之质控
得到标准的3张图，包括主成分分析，高变基因的表达量热图，样品相关性热图
## 下面是画PCA的必须操作，需要看说明书。
gse_number = "GSE66270"
exp <- dat
exp=t(exp)#画PCA图时要求是行名时样本名，列名时探针名，因此此时需要转换
exp=as.data.frame(exp)#将matrix转换为data.frame
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")
dat.pca <- PCA(exp , graph = FALSE)#现在exp最后一列是group_list，需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
# 画图，主成分分析图p2
this_title <- paste0(gse_number,'_PCA')
p2 <- fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = group_list, # color by groups
palette = "Dark2",
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups")+
ggtitle(this_title)+
theme_ggstatsplot()+
theme(plot.title = element_text(size=15,hjust = 0.5))
p2
# 下面是1000_sd热图
library(pheatmap)
cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
n=t(scale(t(dat[cg,])))
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
ac=data.frame(Group=group_list)
rownames(ac)=colnames(n)
# 画图，高变基因的表达量热图p3
p3 <- pheatmap::pheatmap(n,
show_colnames =F,
show_rownames = F,
main = gse_number,
annotation_col=ac,
breaks = seq(-3,3,length.out = 100))#因为已经手动设置了表达量最大值，所以，可以不用设置break
p3
# 画图，样品相关性热图p4
colD=data.frame(Group=group_list)
exprSet=t(exp)
rownames(colD)=colnames(exprSet)#问题-exprSet设置成转置后的exp
p4 <- pheatmap::pheatmap(cor(exprSet),#热图对样本-列 操作
annotation_col = colD,
show_rownames = F,
main = gse_number
)
p4
出图如下:
标准步骤之limma差异分析-配对样本进行配对分析
options(digits = 4) #设置全局的数字有效位数为4
#生成配对信息，注意pairinfo中数字的最大值为样本对数，而数字的排列与exp样本排列顺序一致
m=vector()
for (i in 1:14) {
m=c(m,rep(i,times=2))
}
pairinfo = factor(m)
group = factor(group_list,levels = c("control","ccRCC"))
#差异分析
library(limma)
design <- model.matrix(~group+pairinfo)
exprSet=as.data.frame(dat)
dim(exprSet)
library(tidyr)
a=drop_na(exprSet)
dim(a)
fit <- lmFit(exprSet,design)
dim(fit)
fit2 <- eBayes(fit)
deg = topTable(fit2,coef=2,number=Inf,adjust='fdr')
head(deg)
boxplot(dat[rownames(deg)[1],]~group)
deg[1,]
差异分析结果前10行如下所示：
> deg[1:10,]
logFC AveExpr t P.Value adj.P.Val B
NDUFA4L2 6.725 8.916 51.72 7.568e-19 1.528e-14 32.55
TMPRSS2 -3.093 7.344 -43.68 1.035e-17 1.045e-13 30.39
EGLN3 4.414 8.502 41.61 2.187e-17 1.139e-13 29.75
ERBB4 -5.156 6.270 -41.53 2.257e-17 1.139e-13 29.72
CLCNKB -4.746 8.538 -39.97 4.071e-17 1.644e-13 29.21
IGFBP3 3.701 11.317 38.03 8.756e-17 2.746e-13 28.52
RALYL -4.688 6.257 -37.83 9.522e-17 2.746e-13 28.45
SPAG4 3.905 7.301 36.95 1.369e-16 3.456e-13 28.12
AQP2 -4.766 9.204 -35.91 2.121e-16 4.412e-13 27.73
ANGPTL4 3.971 9.077 35.84 2.185e-16 4.412e-13 27.70
有了差异分析就可以进行标准的可视化，包括火山图和上下调的差异基因热图
代码如下：
nrDEG=deg
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(P.Value))#简单画图看一下
df=nrDEG
df$v= -log10(P.Value) #df新增加一列'v',作为新的绘图参数，值为-log10(P.Value)
#设定上下调基因
df$g=ifelse(df$P.Value>0.05,'stable',
ifelse( df$logFC >2,'up',
ifelse( df$logFC < -2,'down','stable') )
)
#统计上下调基因数量
table(df$g)
#给绘制火山图用的数据新增一列symbol
df$name=rownames(df)
head(df)
logFC_t = 2
#设置可循环使用的plot标题
this_tile <- paste0('Cutoff for logFC is ',round(logFC_t,3),
'\nThe number of up gene is ',nrow(df[df$g == 'up',]) ,
'\nThe number of down gene is ',nrow(df[df$g == 'down',])
)
#画图，火山图p5
p5 <- ggplot(data = df,
aes(x = logFC,
y = -log10(P.Value))) +
geom_point(alpha=0.6, size=1.5,
aes(color=g)) +
ylab("-log10(Pvalue)")+
scale_color_manual(values=c("#34bfb5", "#828586","#ff6633"))+
geom_vline(xintercept= 0,lty=4,col="grey",lwd=0.8) +
xlim(-3, 3)+
theme_classic()+
ggtitle(this_tile )+
theme(plot.title = element_text(size=12,hjust = 0.5),
legend.title = element_blank(),
)
p5
#热图
library(pheatmap)
x=deg$logFC
names(x)=rownames(deg)
cg=c(names(head(sort(x),100)),
names(tail(sort(x),100)))#对x进行从小到大排列，取前100及后100，并取其对应的探针名，作为向量赋值给cg
n=t(scale(t(dat[cg,])))
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(group=group_list)
rownames(ac)=colnames(n)
# 画图，上下调的差异基因热图p6
p6 <- pheatmap(n,show_colnames =F,
show_rownames = F,
cluster_cols = T,
main = gse_number,
annotation_col=ac)
p6
出图如下：
标准步骤之生物学功能数据库注释
我们这里不根据任何武断的阈值来区分统计学显著的上下调基因，而是直接根据基因的变化情况排序进行gsea分析，而且仅仅是展示kegg这个生物学功能数据库的注释情况！
symbol 和 ENTREZID转换
gsea分析需要基因的ENTREZID，根据物种进行转换
# 加ENTREZID列，用于富集分析（symbol转entrezid，然后inner_join）
deg$symbol=rownames(deg)
library(org.Hs.eg.db)
library(clusterProfiler)
s2e <- bitr(deg$symbol,
fromType = "SYMBOL",
toType = "ENTREZID",
OrgDb = org.Hs.eg.db)#人
#bitr()用于SYMBOL转ENTREZID
#其他物种http://bioconductor.org/packages/release/BiocViews.html#___OrgDb
dim(deg)
dim(s2e)
setdiff(deg$symbol,s2e$SYMBOL)
DEG <- inner_join(deg,s2e,by=c("symbol"="SYMBOL"))
gsea富集
library(dplyr)
library(ggplot2)
geneList=DEG$logFC
names(geneList)=DEG$ENTREZID
geneList=sort(geneList,decreasing = T)
head(geneList)
library(clusterProfiler)
kk_gse <- gseKEGG(geneList = geneList,
organism = 'hsa',#按需替换
#nPerm = 1000,
minGSSize = 10,
pvalueCutoff = 0.9,
verbose = FALSE)
tmp=kk_gse@result
dim(tmp)
kk=DOSE::setReadable(kk_gse, OrgDb='org.Hs.eg.db',keyType='ENTREZID')#按需替换
#DOSE::setReadable():mapping geneID to gene Symbol
tmp=kk@result
dim(tmp)
pro='comp1'
write.csv(kk@result,paste0(pro,'_kegg.gsea.csv'))
save(kk,file = 'gsea_kk.Rdata')
富集可视化
上面的kk这个变量就存储了kegg这个生物学功能数据库的gsea分析结果，我们进行简单可视化，代码如下：
# 展现前6个上调通路和6个下调通路
down_k <- kk_gse[tail(order(kk_gse$enrichmentScore,decreasing = F)),];down_k$group=-1
up_k <- kk_gse[head(order(kk_gse$enrichmentScore,decreasing = F)),];up_k$group=1
dat=rbind(up_k,down_k)
colnames(dat)
dat$pvalue = -log10(dat$pvalue)
dat$pvalue=dat$pvalue*dat$group
dat=dat[order(dat$pvalue,decreasing = F),]
# gsea分析结果p7
p7<- ggplot(dat, aes(x=reorder(Description,order(pvalue, decreasing = F)), y=pvalue, fill=group)) +
geom_bar(stat="identity") +
scale_fill_gradient(low="#34bfb5",high="#ff6633",guide = FALSE) +
scale_x_discrete(name ="Pathway names") +
scale_y_continuous(name ="log10P-value") +
coord_flip() +
theme_ggstatsplot()+
theme(plot.title = element_text(size = 15,hjust = 0.5),
axis.text = element_text(size = 12,face = 'bold'),
panel.grid = element_blank())+
ggtitle("Pathway Enrichment")
p7
#具体看上面条形图里面的每个通路的gsea分布情况p8
library(enrichplot)
p8 <- gseaplot2(kk, geneSetID = rownames(down_k))+
gseaplot2(kk, geneSetID = rownames(up_k))
p8
如果你也有类似的数据分析需求，却苦于不会写代码，可以考虑找我们的工程师帮忙哦！
转录组产品线
明码标价之转录组常规测序服务（仅需799每个样品） 明码标价之普通转录组上游分析 明码标价之转录组下游分析仅需800元 明码标价之转录组测序数据的可变剪切 明码标价之RNA-Seq数据的内含子保留分析