事情是这样的：曾老板分配了一个任务，完成100个GEO数据库的基因芯片数据分析，一开始我可能觉得这个事情比较简单，心想就我都干了这么多年了。然后吭哧吭哧在完成了几十个数据集后（已经被批评了速度太慢，距离任务下来过的时间比较久，后面还是得踏实下来好好干！），老板喊我去对结果。
就眼看老板一顿操作猛如虎，迅雷不及掩耳，就随意一看，就挑出了8个数据集的结果是有问题的！
我可能跑数据分析结果都没有走心去看，就这，我可能都得学个好几年。
下面来看第一个数据：
从cel文件开始，看看矩阵相关性会不会变化。
https://ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE17351
我从头梳理了一下这个数据，
这次从cel芯片数据开始分析
，在这个过程中，我带入了思考，遇到了以下几个问题：
-
问题1：数据是非常经典的GPL570平台，众所周知有2w个左右的基因，但是这个数据的matrix表达矩阵存在很多NA值，处理完后基因水平的表达矩阵只剩下14080个基因
-
问题2：CEL文件处理的表达与matrix表达矩阵的样本相关性分析 会有什么变化吗？
首先从cel开始处理：
library(AnnoProbe)
library(GEOquery)
library(ggplot2)
library(ggstatsplot)
library(patchwork)
library(reshape2)
library(stringr)
library(limma)
library(tidyverse)
# BiocManager::install("oligo", configure.args="--disable-threading", force = TRUE)
library(oligo)
packageVersion("oligo")
library(pd.hg.u133.plus.2)
getOption('timeout')
options(timeout=10000)
## 1.获取并且检查表达量矩阵
## ～～～gse编号需修改～～～
gse_number "GSE17351"
dir.create(gse_number)
setwd(gse_number)
getwd()
list.files()
## 3.表达矩阵处理: 使用cel文件
# https://ftp.ncbi.nlm.nih.gov/geo/series/GSE33nnn/GSE33479/suppl/GSE33479_RAW.tar
url "https://ftp.ncbi.nlm.nih.gov/geo/series/GSE17nnn/GSE17351/suppl/GSE17351_RAW.tar"
download.file(url = url, destfile = "GSE17351_RAW.tar")
## 读取芯片原始cell文件
# options(BioC_mirror="https://mirrors.westlake.edu.cn/bioconductor")
configure.args "--disable-threading")
options(configure.args = configure.args)
celFiles "GSE17351_RAW/",listGzipped=T, full.name=TRUE)
celFiles
rawData eset ## 一步完成背景校正、分位数校正标准化和RMA(robust multichip average) 表达整合
e colnames(e) ".CEL.gz","",colnames(e))
dim(e)
e[1:4,1:4]
range(e)
拿到表达矩阵就走后面的分析了，这里遇到一个报错：
使用oligo包处理的时候报错
：
eset Background correcting
Error in basicRMA(pms, pnVec, normalize, background) :
ERROR; return code from pthread_create() is 22
解决办法：
kimi告诉我：针对您在使用
rma
函数时遇到的
ERROR; return code from pthread_create() is 22
错误，根据搜索结果，这是一个与多线程创建相关的问题。这个错误通常发生在尝试创建新线程时，特别是与
pthread_create()
函数有关。
将oligo包更新到了：‘1.70.0’版本，并在代码运行前设置：
configure.args "--disable-threading")
options(configure.args = configure.args)
cel得到表达的样本相关性分析
#~~~~~corplot~~~~~
exprSet=dat
pheatmap::pheatmap(cor(exprSet))
# 组内的样本的相似性应该是要高于组间的！
colD=data.frame(Group=group_list)
rownames(colD)=colnames(exprSet)
pheatmap::pheatmap(cor(exprSet),
annotation_col = colD,
show_rownames = F)
#～～～样品相关性热图p4～～～
p4 color = colorRampPalette(c("#34bfb5", "white","#ff6633" ))(100),
annotation_col = colD,
show_rownames = T,
show_colnames = T,
display_numbers=T,
main = paste(pro,"correlation",sep = "_"),
)
p4
ggsave('qc_corplot.pdf',plot = p4,width = 5,height = 5)
上面的表达量矩阵是我自己读取cel文件后使用oligo包的rma算法拿到的，
结果如下：log2-transformed RMA signal intensities
matrix矩阵样本相关性分析
其实作者本身也有自己的处理方式，可以直接下载作者的表达量矩阵文件，然后看 matrix中看数据集，里面的处理方法：
log2-transformed GC-RMA signal intensities
结果：
前后对比这两个相关性热图可以看出来：两次相关性分析有一些不一样，cel处理的相关性变高了，但是组内与组间样本相关性都高了，样本间的聚类关系并有没改变。
处理matrix矩阵时的插曲：这么多NA值吗？
## 1.获取并且检查表达量矩阵
## ～～～gse编号需修改～～～
gse_number "GSE17351"
dir.create(gse_number)
setwd(gse_number)
getwd()
list.files()
#gset
gset '.', getGPL = T)
gset[[1]]
## 2.根据生物学背景及研究目的人为分组
## 通过查看说明书知道取对象a里的临床信息用pData
## 挑选一些感兴趣的临床表型。
a pd colnames(pd)
pd$title
pd$characteristics_ch1
table(pd$characteristics_ch1)
## ～～～分组信息编号需修改～～～
group_list 'Normal',pd$title ,ignore.case = T ), "control","case")
group_list "control","case"))
table(group_list)
## 3.表达矩阵处理
dat # a现在是一个对象，取a这个对象通过看说明书知道要用exprs这个函数
dim(dat) # 看一下dat这个矩阵的维度
dat[1:4,1:4] # 查看dat这个矩阵的1至4行和1至4列，逗号前为行，逗号后为列
range(dat)
dat range(dat)
可以看到很多NA行：
我第一次遇到这种矩阵，非常经典的平台，为什么有这么多NA行探针，处理完后基因只剩下：
