【原】在什么情況下癌癥樣本與正常樣本是分不開的呢？

健明 2025-01-02 發(fā)布于廣東

展開全文

在老板給出的100個GEO芯片數(shù)據(jù)分析中，我們前面已經(jīng)介紹過四個有意思的數(shù)據(jù)集了，大家可以先去讀一下，詳見：

正常組織與癌旁組織可以一視同仁嗎？
2萬個基因少一半也不影響最后的差異分析富集結(jié)果啊？
火山圖展示差異分析上下調(diào)基因數(shù)量，不平衡無處不在嗎？
花了10多萬的隊列就只為了這一張圖嗎？

今天遇到的數(shù)據(jù)集合為：https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE54236，數(shù)據(jù)情況如下：

首先還是標(biāo)準(zhǔn)的芯片數(shù)據(jù)分析得到表達矩陣：

## 魔幻操作，一鍵清空~
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現(xiàn)在是一個對象，取a這個對象通過看說明書知道要用exprs這個函數(shù)
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)

## 根據(jù)生物學(xué)背景及研究目的人為分組
## 通過查看說明書知道取對象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.多對一取均值
# 合并探針I(yè)D 與基因，表達譜對應(yīng)關(guān)系
# 提取表達矩陣
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')

得到的表達矩陣如下：

進行簡單的質(zhì)控，繪制三張經(jīng)典的質(zhì)控圖：

## 數(shù)據(jù)檢查
load( file = 'step1_output.Rdata')
source('../scripts/run_check.R')
run_check(dat, group_list, target_gene='AFP', pro=gse_number,path='./')

結(jié)果如下，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)

腫瘤與正常癌旁樣本的差異基因數(shù)量非常少：

如果只用差異基因做PCA分析：

exp <- dat[deg$g!="stable",]
## 下面是畫PCA的必須操作，需要看說明書。
exp=t(exp)#畫PCA圖時要求是行名時樣本名，列名時探針名，因此此時需要轉(zhuǎn)換
exp=as.data.frame(exp)#將matrix轉(zhuǎn)換為data.frame 

library("FactoMineR")#畫主成分分析圖需要加載這兩個包
library("factoextra")  
#～～～主成分分析圖p2～～～
dat.pca <- PCA(exp , graph = FALSE)#現(xiàn)在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

如果使用sd值指標(biāo)挑選：

cg <- names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列?。┤∶恳恍械姆讲?，從小到大排序，取最大的1000個

exp <- dat[cg,]
## 下面是畫PCA的必須操作，需要看說明書。
exp=t(exp)#畫PCA圖時要求是行名時樣本名，列名時探針名，因此此時需要轉(zhuǎn)換
exp=as.data.frame(exp)#將matrix轉(zhuǎn)換為data.frame 

library("FactoMineR")#畫主成分分析圖需要加載這兩個包
library("factoextra")  
#～～～主成分分析圖p2～～～
dat.pca <- PCA(exp , graph = FALSE)#現(xiàn)在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可能無法有效地將它們區(qū)分開來。這可能是因為腫瘤和癌旁組織在分子層面的異質(zhì)性較低，或者腫瘤微環(huán)境與正常組織的差異不大。
樣本量不足：如果分析的樣本數(shù)量較少，PCA可能無法捕捉到足夠的變異性來區(qū)分腫瘤和癌旁樣本。樣本量的不足可能導(dǎo)致分析結(jié)果的不穩(wěn)定和不準(zhǔn)確。
數(shù)據(jù)質(zhì)量問題：如果輸入PCA的數(shù)據(jù)質(zhì)量不高，例如存在大量的噪聲或者技術(shù)偏差，這可能會掩蓋腫瘤和癌旁樣本之間的真實差異，導(dǎo)致PCA分析無法有效分離樣本。
預(yù)處理不當(dāng)：在進行PCA之前，數(shù)據(jù)需要適當(dāng)?shù)念A(yù)處理，包括標(biāo)準(zhǔn)化、去除批次效應(yīng)等。如果預(yù)處理步驟不當(dāng)，可能會影響PCA分析的結(jié)果，使得腫瘤和癌旁樣本難以區(qū)分。
腫瘤微環(huán)境的復(fù)雜性：腫瘤微環(huán)境可能包含多種細胞類型，包括免疫細胞、基質(zhì)細胞等，這些細胞的存在可能使得腫瘤樣本的基因表達模式與癌旁樣本相似，從而難以通過PCA區(qū)分。
腫瘤異質(zhì)性：腫瘤本身的異質(zhì)性也可能導(dǎo)致PCA分析中腫瘤樣本與癌旁樣本難以區(qū)分。不同腫瘤樣本之間可能存在顯著的基因表達差異，這使得它們在PCA分析中難以與癌旁樣本區(qū)分開來。
分析方法的選擇：PCA是一種線性降維技術(shù)，對于非線性關(guān)系可能不夠敏感。如果腫瘤和癌旁樣本之間的差異是非線性的，可能需要考慮使用其他更復(fù)雜的分析方法。

總結(jié)來說，腫瘤樣本和癌旁樣本在PCA分析中分不開可能是由于樣本間差異性不足、樣本量不足、數(shù)據(jù)質(zhì)量問題、預(yù)處理不當(dāng)、腫瘤微環(huán)境的復(fù)雜性、腫瘤異質(zhì)性以及分析方法的選擇等多種因素造成的。

學(xué)徒作業(yè)

老板給我發(fā)了另外一個數(shù)據(jù)集，也是存在腫瘤樣本與癌旁樣本分不開的情況，數(shù)據(jù)來自文獻《A Candidate Molecular Biomarker Panel for the Detection of Bladder Cancer》，數(shù)據(jù)編號為：https://www.ncbi.nlm./geo/query/acc.cgi?acc=GSE31189，

使用上面的代碼，去做做看。