感谢生信技能树小洁老师

Section1

Illumina数据的处理的关键是获取表达矩阵

#数据下载
rm(list = ls())
options(stringsAsFactors = F)
library(GEOquery)
gse_number = "GSE100748"
eSet <- getGEO(gse_number, 
               destdir = '.', 
               getGPL = F)
class(eSet)
length(eSet)
eSet = eSet[[1]]
#(1)提取表达矩阵exp
exp <- exprs(eSet)
exp[1:4,1:4]
exp = log2(exp+1) #如果是已经取了log，可以把该行代码注释掉，避免某个数据为0，log2（0）为负无穷
boxplot(exp)

如果直接使用数据下载的代码下载数据，得到的表达矩阵是经过标准化处理过的，里面存在负值。
所以需要采取另一种办法，按照jimmy老师的分享，用lumi包来解决这个问题>http://www.bio-info-trainee.com/1944.html

if (!requireNamespace("BiocManager", quietly = TRUE)) install.packages("BiocManager")
if (!requireNamespace("lumi", quietly = TRUE)) BiocManager::install("lumi")
if (!requireNamespace("methylumi", quietly = TRUE)) BiocManager::install("methylumi")
BiocManager::install("GenomicFeatures")
BiocManager::install("methylumi")
BiocManager::install("lumi")
library(GenomicFeatures)
library(methylumi)
library(lumi)#在安装lumi包的时候遇到一堆问题，methylumi是lumi的关联包，反正搞了好久也还是装上去了。。。。
## 首先是从illumina的芯片结果文件，自己用R的lumi包来获取表达矩阵。
fileName <- 'GSE100748_non-normalized.txt.gz' # 从GEO中下载该文件
x.lumi <- lumiR.batch(fileName) ##, sampleInfoFile='sampleInfo.txt')
pData(phenoData(x.lumi))
## Do all the default preprocessing in one step
lumi.N.Q <- lumiExpresso(x.lumi)
### retrieve normalized data
dataMatrix <- exprs(lumi.N.Q) #dataMatrix即表达矩阵
exp = dataMatrix#即可无缝对接pipeline 01
boxplot(exp) #看一下表达量数据分布是否一致

Section2

对数据进行分组，由于GSE100748包含未分化BMSC、成脂、成骨、成软骨等多组细胞。在此，只关注成脂0d、7d、21d三个分组，同时以0d作为control。

#(2)提取临床信息
pd <- pData(eSet)
#在数据集中只提取脂肪生成相关的，并分三组，0d、07d、21d

zerod = pd[pd$`cell type:ch1` %in% "Undifferentiated MSCs",]#直接写0d、7d、21d不行，不能作为变量
sevend = pd[pd$`cell type:ch1` %in% "Adipocytes Day 7",]
twentyoned = pd[pd$`cell type:ch1` %in% "Adipocytes Day 21",]

ad=rbind(zerod,sevend,twentyoned)
pd = ad

#pipeline02 进行分组
 Group(实验分组)和ids(探针注释)
rm(list = ls())  
load(file = "step1output.Rdata")
library(stringr)

#匹配关键词，自行分类
#Group=ifelse(str_detect(pd$source_name_ch1,"control"),"control","RA")
Group=ifelse(str_detect(pd$`cell type:ch1`,"Undifferentiated MSCs"),"0d",
             ifelse(str_detect(pd$`cell type:ch1`,"Adipocytes Day 7"),"7d","21d"))
#设置参考水平，指定levels，目的是要把对照组在前，处理组在后
Group = factor(Group,
               levels = c("0d","7d","21d"))
Group

pd用于分组的部分

ids将探针名换成基因名，貌似因为是illumina数据的原因，方法一中没有该GPL

http://www.bio-info-trainee.com/1399.html
需要用方法二自行注释

方法2 读取GPL平台的soft文件，按列取子集

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

if(T){
  #注：soft文件列名不统一，活学活用，有的GPL平台没有提供注释，如GPL16956
  a = getGEO(gpl_number,destdir = ".")
  b = a@dataTable@table
  colnames(b)
  ids2 = b[,c("ID","Symbol")]#有的数据不一定写Symbol，可能写的是Gene Symbol，如果出现报错，可以在colname(b)里查看到底哪一列是注释
  colnames(ids2) = c("probe_id","symbol")
  ids2 = ids2[ids2$symbol!="" & !str_detect(ids2$symbol,"///"),]
  ids = ids2 
  }

热图

pheatmap(n2,
         show_colnames =F,
         show_rownames = F,
         annotation_col=annotation_col,
         color = colorRampPalette(colors = c("blue","white","red"))(100),#热图的颜色
         breaks = seq(-4,4,length.out = 100)#注释条的取值范围，-4、4分别对应颜色最深
)

其他的pipeline（03以后）没有什么明显需要修改的地方了，小洁老师上课时候说火山图必须两组两组之间对比，但你要是不更改，照样可以画出火山图但应该是有问题的，变为了control组和实验组（另外两组被自动划为一组了吧大概）之间的火山图。（那我看的那篇文章又有问题了哦）。