今天遇到的数据集合为:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE54236,数据情况如下:
## 魔幻操作,一键清空~
rm(list = ls())
library(AnnoProbe)
library(GEOquery)
library(ggplot2)
library(ggstatsplot)
library(patchwork)
library(reshape2)
library(stringr)
library(limma)
library(tidyverse)
getOption('timeout')
options(timeout=10000)
## 获取并且检查表达量矩阵
## ~~~gse编号需修改~~~
gse_number <- "GSE54236"
dir.create(gse_number)
setwd(gse_number)
getwd()
list.files()
#gset <- geoChina(gse_number)
gset <- getGEO(gse_number, destdir = '.', getGPL = T)
gset[[1]]
a <- gset[[1]]
dat <- exprs(a) # a现在是一个对象,取a这个对象通过看说明书知道要用exprs这个函数
dim(dat) # 看一下dat这个矩阵的维度
dat[1:4,1:4] # 查看dat这个矩阵的1至4行和1至4列,逗号前为行,逗号后为列
range(dat)
table(is.na(dat))
#dat[is.na(dat)] <- 0
dat <- na.omit(dat)
range(dat)
dim(dat)
## 根据生物学背景及研究目的人为分组
## 通过查看说明书知道取对象a里的临床信息用pData
## 挑选一些感兴趣的临床表型。
pd <- pData(a)
colnames(pd)
pd$title
pd$characteristics_ch1
table(pd$characteristics_ch1)
## ~~~分组信息编号需修改~~~
group_list <- ifelse(grepl('Non tumor',pd$title ,ignore.case = T ), "control","case")
group_list <- factor(group_list, levels = c("control","case"))
table(group_list)
## 查看注释平台gpl 获取芯片注释信息
a
a@annotation
gpl <- getGEO(filename=paste0(a@annotation,".soft.gz"))
class(gpl)
gpl_anno <- gpl@dataTable@table
colnames(gpl_anno)
id2name <- gpl_anno[,c("ID" ,"GENE_SYMBOL")]
colnames(id2name) <- c("ID","GENE_SYMBOL")
# 1.过滤掉空的探针
id2name <- na.omit(id2name)
id2name <- id2name[which(id2name$GENE_SYMBOL!=""), ]
# 2.过滤探针一对多
id2name <- id2name[!grepl("\\///",id2name$GENE_SYMBOL), ]
head(id2name)
# 3.多对一取均值
# 合并探针ID 与基因,表达谱对应关系
# 提取表达矩阵
dat <- dat %>%
as.data.frame() %>%
rownames_to_column("ID")
exp <- merge(id2name, dat, by.x="ID", by.y="ID")
# 多对一取均值
exp <- avereps(exp[,-c(1,2)],ID = exp$GENE_SYMBOL) %>%
as.data.frame()
dat <- as.matrix(exp[,pd$geo_accession])
dim(dat)
fivenum(dat['CRP',])
fivenum(dat['GAPDH',])
dat[1:5, 1:6]
save(gse_number, dat, group_list, pd, file = 'step1_output.Rdata')
得到的表达矩阵如下:
进行简单的质控,绘制三张经典的质控图:
## 数据检查
load( file = 'step1_output.Rdata')
source('../scripts/run_check.R')
run_check(dat, group_list, target_gene='AFP', pro=gse_number,path='./')
结果如下,PCA图中 可以看到 肿瘤样本与对照样本全完的混在了一起
## 人类看家基因列表如下(部分基因,最常用的是CRP和CRP)
source('../scripts/run_DEG.R')
group_list
deg <- run_DEG(dat, group_list, target_gene=c('AFP'),pro=gse_number,path='./')
head(deg)
table(deg$g)
肿瘤与正常癌旁样本的差异基因数量非常少:
exp <- dat[deg$g!="stable",]
## 下面是画PCA的必须操作,需要看说明书。
exp=t(exp)#画PCA图时要求是行名时样本名,列名时探针名,因此此时需要转换
exp=as.data.frame(exp)#将matrix转换为data.frame
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")
#~~~主成分分析图p2~~~
dat.pca <- PCA(exp , graph = FALSE)#现在exp最后一列是group_list,需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
pro <- gse_number
this_title <- paste0(pro,'_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=12,hjust = 0.5))
p2
cg <- names(tail(sort(apply(dat,1,sd)),1000))#apply按行('1'是按行取,'2'是按列取)取每一行的方差,从小到大排序,取最大的1000个
exp <- dat[cg,]
## 下面是画PCA的必须操作,需要看说明书。
exp=t(exp)#画PCA图时要求是行名时样本名,列名时探针名,因此此时需要转换
exp=as.data.frame(exp)#将matrix转换为data.frame
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")
#~~~主成分分析图p2~~~
dat.pca <- PCA(exp , graph = FALSE)#现在exp最后一列是group_list,需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
pro <- gse_number
this_title <- paste0(pro,'_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=12,hjust = 0.5))
p2
要是使用差异最大的那部分基因,还是强行可以勉强分开的。
在使用PCA(主成分分析)分析肿瘤样本与癌旁样本时,肿瘤样本和癌旁样本可能在以下情况下分不开:
总结来说,肿瘤样本和癌旁样本在PCA分析中分不开可能是由于样本间差异性不足、样本量不足、数据质量问题、预处理不当、肿瘤微环境的复杂性、肿瘤异质性以及分析方法的选择等多种因素造成的。
老板给我发了另外一个数据集,也是存在肿瘤样本与癌旁样本分不开的情况,数据来自文献《A Candidate Molecular Biomarker Panel for the Detection of Bladder Cancer》,数据编号为:https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE31189,
使用上面的代码,去做做看。