【原】芯片數(shù)據(jù)分析，so easy?

健明 2021-07-14

展開全文

我最早接觸的高通量數(shù)據(jù)就是RNA-seq,，后來接觸的也基本是高通量測序結(jié)果而不是芯片數(shù)據(jù)，因此我從來沒有分析過一次芯片數(shù)據(jù),，而最近有一個學(xué)員在看生信技能樹在騰訊課堂發(fā)布的課程GEO數(shù)據(jù)庫表達芯片處理之R語言流程（閱讀原文購買）遇到了問題問我請教,，為了解決這個問題，我花了一個晚上時間學(xué)習(xí)這方面的分析,。

注:這篇文章不會介紹R語言的安裝和使用,，也不會介紹GEO數(shù)據(jù)庫的結(jié)構(gòu)，默認(rèn)你都知道,，不知道可以搜索Jimmy的教程,，或者購買視頻學(xué)習(xí)。

數(shù)據(jù)的獲取

數(shù)據(jù)獲取有兩種方式,，R包GEOquery解析和手動下載,。其中前面一種最方便，完成了手動數(shù)據(jù)下載和Bioconductor常見數(shù)據(jù)結(jié)構(gòu)ExpressionSet的構(gòu)造,，關(guān)于這個數(shù)據(jù)結(jié)構(gòu)的具體介紹看Bioconductor的介紹或者視頻,，簡言之,，就是用于存放 實驗信息, 分組信息 和 表達信息, 方便后續(xù)調(diào)用。

library(GEOquery)
gset <- getGEO("GSE13535", GSEMatrix =TRUE, AnnotGPL=TRUE )
show(gset)

一般而言GEOquery解析都是首選,，除非你提供的GSE號還沒被GEOquery記錄或者說網(wǎng)絡(luò)速度感人,，以及你不覺得別人提供的矩陣是你所需要的，你才會決定去手工下載,。分為兩種情況,，一種是下載賽默飛的下機原始數(shù)據(jù)格式CEL，一種是下載單個樣本表達量向量或者含有所有樣本的表達量矩陣,。

先說第一種,，可以直接點擊http下載到tar打包的數(shù)據(jù), 然后解壓縮得到所有的CEL文件

setwd("F:/Project/GEO_project/")library(affy)
affy.data <- ReadAffy()
length(affy.data)
# 13
eset.rma <- rma(affy.data)
exprSet <- exprs(eset.rma)
write.table(exprSet, "expr_rma_matrix.txt", quote=F, sep="\t")

ReadAffy: 讀取當(dāng)前文件下的CEL格式文件，同時第一次還會從bioconductor上下載hugene10stv1用來注釋cel文件,。
rma: 基于robust multi-arrary average(RMA)算法衡量表達量,，從而將AffyBatch對象轉(zhuǎn)換成ExpressionSet
exprs: 獲取ExpressionSet中的表達量矩陣
write.table: 將表達量矩陣信息保存到本地

然后是第二種，以所有樣本的表達矩陣為例,，可以用瀏覽器到ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42589/matrix/下載,，如果你會用Linux的話，可以用wget -4 ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE42nnn/GSE42589/matrix/GSE42589_series_matrix.txt.gz, 才1.7M,。解壓縮這個文件后,，有一個txt文件, 這個txt分為兩個部分。第一個部分是以”!”開頭的樣本的所有信息,，如實驗平臺,、處理、以及分組等信息,。第二個部分則是后面的表達量信息,，

expr.df <- read.table(file = "GSE42589_series_matrix.txt", header =TRUE,
comment.char = "!", row.names=1)

可以從這個角度理解這三種方法：最開始得到的都是CEL文件，CEL文件需要一系列的步驟才能轉(zhuǎn)換成表達矩陣,，例如去除批次效應(yīng),、質(zhì)控和過濾等,，得到的表達矩陣在上傳時會增加元數(shù)據(jù)信息（處理方法,、分組信息），就成為我們下載的GSEXXXX_series_matrix.txt.gz. 通過手工解析加R語言簡單操作得到了R語言中的數(shù)據(jù)框(data.frame),，而GEOquery能夠幫助我們完成下載和解析這兩個步驟,。

三者的優(yōu)先級為：GEOquery > 手工下載表達量矩陣文件 > 手工下載原始的CEL文件。

使用limma進行差異表達分析

limma的核心函數(shù)是lmFit和eBayes,，前者是用于線性擬合,，后者根據(jù)前者的擬合結(jié)果進行統(tǒng)計推斷。

lmFit至少需要兩個輸入,，一個是表達矩陣,，一個是分組對象,。

表達矩陣必須是matrix類數(shù)據(jù)結(jié)構(gòu)，每一列都是存放一個樣本,，每一行是一個探針信息或者是注釋后的基因名,。這里就是向我提問的人出錯的原因，他在讀入數(shù)據(jù)時,，read.table少了參數(shù),，row.names= 1，導(dǎo)致第一列是探針信息,。

# 使用GEOquery
exprSet <- exprs(gset[[1]])
# 基于matrix
expr.df <- read.table(file = "GSE42589_series_matrix.txt", header =TRUE,
comment.char = "!", row.names=1)
# 從cel文件開始
exprSet <- exprs(eset.rma)

試驗設(shè)計矩陣: 沒有試驗設(shè)計矩陣對象,，limma就不知道如何比較。分組數(shù)據(jù)可以手工從之前的matrix.gz整理,，整理到一個excel,，然后用R讀取，或者就是直接從Geoquery的結(jié)果中解析,。

pData <- pData(gset[[1]])
view(pData)

其中title部分告訴了我們分組信息,，2小時和18小時，每個時間段又有vehicle control, PE1.3 embolized, PE2.0 embolized,，也就是2x2的雙因素試驗設(shè)計, 我們可以現(xiàn)在R語言里構(gòu)建實驗設(shè)計的數(shù)據(jù)框,。

sample <- pData$geo_accession
treat_time <- rep(c("2h","18h"),each=11)
treat_type <- rep(rep(c("vehicle_control","PE1.3_embolized","PE2.0_embolized"), c(3,4,4)),
times=2)
design_df <- data.frame(sample, treat_time, treat_type)

根據(jù)Limma的使用手冊的”9.5 Interaction Models: 2 X 2 Factorial Design”進行手續(xù)的分析。這里僅僅展示單個因素的分析過程,，多個因素看文檔依葫蘆畫瓢就行,。

構(gòu)建單因素試驗設(shè)計矩陣，進行線性擬合

TS <- paste(design_df$treat_time, design_df$treat_type, sep=".")
TS
TS <- factor(TS, levels = unique(TS))
design <- model.matrix(~0+TS)
fit <- lmFit(exprSet, design)

然后根據(jù)我們要回答的問題,，來設(shè)置比較對象,。比如不同時間段下控制組哪些基因發(fā)生了差異報答，處理18小時后,，處理組相對于對照組有哪些基因發(fā)生差異表達,，也就是做多次t檢驗。

cont.matrix <- makeContrasts(
vs1  = TS18.vehicle_control-TS2.vehicle_control, # 對照組在前后的差異表達基因
vs2  = TS18.PE2.0_embolized-TS2.PE2.0_embolized, # PE2.0處理前后的差異基因
vs3  = TS18.PE1.3_embolized-TS2.PE1.3_embolized, # PE1.3在處理前后差異基因# 處理18小時候,，PE2.0相對于對照變化的基因再與PE1.3與對照的差異比較
diff = (TS18.PE2.0_embolized-TS18.vehicle_control)-(TS18.PE1.3_embolized-TS18.vehicle_control),
levels = design
)

fit2 <- contrasts.fit(fit, cont.matrix)
results <- decideTests(fit2)

最后的結(jié)果可以用韋恩圖展示vennDiagram(results)

久久国产成人av_抖音国产毛片_a片网站免费观看_A片无码播放手机在线观看,色五月在线观看,亚洲精品m在线观看,女人自慰的免费网址,悠悠在线观看精品视频,一级日本片免费的,亚洲精品久,国产精品成人久久久久久久

【原】芯片數(shù)據(jù)分析，so easy?

數(shù)據(jù)的獲取

使用limma進行差異表達分析

更多分析