花了10多万的队列就只为了这一张图吗？

生信技能树

发布于 2024-12-27 19:56:59

7900

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

前面已经给出了3个GEO芯片数据挖掘分析点，详见：

现在继续完成老板分配的任务，《100个GEO芯片数据分析》，真的是信息量很大啊。又遇到了一个有意思的芯片数据，数据如下：

https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE44861，

作者一下子做了111个样本的芯片测序，那可是在2013年提交进GEO的，说明芯片测序可能在更早的时候，那个时候单个表达量芯片样品起码得两三千块钱人民币，这个队列保守估计仅仅是芯片费用十几万了。

样本表型包括：111 colon tissues from tumors and adjacent noncancerous tissues

然后按照芯片的标准分析，如下：

首先是获取样本分组：

## 魔幻操作，一键清空~
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 <- "GSE44861"
dir.create(gse_number)
setwd(gse_number)
getwd()
list.files() 

#gset <- geoChina(gse_number)
gset <- getGEO(gse_number, destdir = '.', getGPL = T)
gset[[1]]

#######
## 根据生物学背景及研究目的人为分组
## 通过查看说明书知道取对象a里的临床信息用pData
## 挑选一些感兴趣的临床表型。
a <- gset[[1]]
pd <- pData(a)
colnames(pd)
pd$title
pd$characteristics_ch1.3
table(pd$characteristics_ch1.3)
table(pd$`tissue:ch1`)
table(pd$`rs5995355 genotype:ch1`)
table(pd$`tissue:ch1`,pd$`rs5995355 genotype:ch1`)

## ～～～分组信息编号需修改～～～
group_list <- ifelse(grepl('N',pd$title ,ignore.case = T ), "normal","tumor")
group_list <- factor(group_list, levels = c("normal","tumor"))
table(group_list)

# normal  tumor 
#    55     56

共有55个癌旁正常样本，56个肿瘤样本。（非常标准的实验设计，理论上癌症组织跟癌旁肯定是差异很大），前面我们就讨论过：

正常组织与癌旁组织可以一视同仁吗

然后是提取芯片表达矩阵，探针id转换基因symbol：

#######################################################
a <- gset[[1]]
dat <- exprs(a) # a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)        # 看一下dat这个矩阵的维度
dat[1:4,1:4]    # 查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
range(dat)

## 查看注释平台gpl 获取芯片注释信息
a
gpl <- getGEO(filename="GPL3921.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')

吭哧吭哧一顿出图：

经典三张图：PCA，sd top1000 基因热图、样本相关性热图

## 数据检查
load( file = 'step1_output.Rdata')
source('../scripts/run_check.R')
# ERα基因的symbol是ESR1
run_check(dat, group_list, target_gene='NCF4', pro=gse_number,path='./')

可以看到这个数据集的生物学分组不是很好，两个分组（癌症和癌旁）里面的病人样品在PCA图根本分不开呀！

接着我对肿瘤样本 vs 癌旁正常样本做了差异分析：

## 人类看家基因列表如下（部分基因，最常用的是CRP和CRP）
# ERα基因的symbol是ESR1
source('../scripts/run_DEG.R')
group_list
deg <- run_DEG(dat, group_list, target_gene=c('NCF4'),pro=gse_number,path='./')
head(deg)
table(deg$g)

差异结果如下：