感謝生信技能樹小潔老師

Section1

Illumina數(shù)據(jù)的處理的關(guān)鍵是獲取表達矩陣

#數(shù)據(jù)下載
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) #如果是已經(jīng)取了log，可以把該行代碼注釋掉，避免某個數(shù)據(jù)為0，log2（0）為負(fù)無窮
boxplot(exp)

如果直接使用數(shù)據(jù)下載的代碼下載數(shù)據(jù)，得到的表達矩陣是經(jīng)過標(biāo)準(zhǔn)化處理過的，里面存在負(fù)值。
所以需要采取另一種辦法，按照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的關(guān)聯(lián)包，反正搞了好久也還是裝上去了。。。。
## 首先是從illumina的芯片結(jié)果文件，自己用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) #看一下表達量數(shù)據(jù)分布是否一致

Section2

對數(shù)據(jù)進行分組，由于GSE100748包含未分化BMSC、成脂、成骨、成軟骨等多組細(xì)胞。在此，只關(guān)注成脂0d、7d、21d三個分組，同時以0d作為control。

#(2)提取臨床信息
pd <- pData(eSet)
#在數(shù)據(jù)集中只提取脂肪生成相關(guān)的，并分三組，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)

#匹配關(guān)鍵詞，自行分類
#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"))
#設(shè)置參考水平，指定levels，目的是要把對照組在前，處理組在后
Group = factor(Group,
               levels = c("0d","7d","21d"))
Group

pd用于分組的部分

ids將探針名換成基因名，貌似因為是illumina數(shù)據(jù)的原因，方法一中沒有該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文件列名不統(tǒng)一，活學(xué)活用，有的GPL平臺沒有提供注釋，如GPL16956
  a = getGEO(gpl_number,destdir = ".")
  b = a@dataTable@table
  colnames(b)
  ids2 = b[,c("ID","Symbol")]#有的數(shù)據(jù)不一定寫Symbol，可能寫的是Gene Symbol，如果出現(xiàn)報錯，可以在colname(b)里查看到底哪一列是注釋
  colnames(ids2) = c("probe_id","symbol")
  ids2 = ids2[ids2$symbol!="" & !str_detect(ids2$symbol,"http:///"),]
  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分別對應(yīng)顏色最深
)

其他的pipeline（03以后）沒有什么明顯需要修改的地方了，小潔老師上課時候說火山圖必須兩組兩組之間對比，但你要是不更改，照樣可以畫出火山圖但應(yīng)該是有問題的，變?yōu)榱薱ontrol組和實驗組（另外兩組被自動劃為一組了吧大概）之間的火山圖。（那我看的那篇文章又有問題了哦）。