現(xiàn)在繼續(xù)完成老板分配的任務(wù)，《100個GEO芯片數(shù)據(jù)分析》，真的是信息量很大啊。又遇到了一個有意思的芯片數(shù)據(jù)，數(shù)據(jù)如下：

https://www.ncbi.nlm./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]]

#######
## 根據(jù)生物學(xué)背景及研究目的人為分組
## 通過查看說明書知道取對象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個腫瘤樣本。（非常標準的實驗設(shè)計，理論上癌癥組織跟癌旁肯定是差異很大），前面我們就討論過：

• 正常組織與癌旁組織可以一視同仁嗎

然后是提取芯片表達矩陣，探針id轉(zhuǎn)換基因symbol：

#######################################################
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)

## 查看注釋平臺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.多對一取均值
# 合并探針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')

吭哧吭哧一頓出圖：

經(jīng)典三張圖：PCA，sd top1000 基因熱圖、樣本相關(guān)性熱圖

## 數(shù)據(jù)檢查
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='./')

可以看到這個數(shù)據(jù)集的生物學(xué)分組不是很好，兩個分組（癌癥和癌旁）里面的病人樣品在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)

差異結(jié)果如下：

然后老板看完我的結(jié)果，進行了靈魂發(fā)問：

你這個分組有問題呀！你跟作者的文獻中的結(jié)果對比了嗎？作者做了什么分析嗎？可以對的上嗎？

我：我說作者沒有沒有做差異分析

老板：這么大的數(shù)據(jù)隊列，作者為什么沒有做差異分析呢？

我：愣在原地不知所措，哈哈哈哈哈哈哈哈哈。。。（其實我想反駁但是不敢發(fā)話）

來看看，作者文獻里面的結(jié)果分析（全部結(jié)果只有這一樣簡單的圖）：

靈魂發(fā)問：你說為什么呢？（畢竟這個數(shù)據(jù)花了差不多10萬塊呢?。?/span>