需求
和這一篇的需求一致仇哆，不同公司的芯片原始數(shù)據(jù)處理方式也不同醒叁，最終都是為了得到表達(dá)矩陣，今天介紹的是illumina beadchip 芯片消痛，用神奇的limma搞定它且叁。
找不到某些GEO數(shù)據(jù)的表達(dá)矩陣腫么辦
illumina beadchip芯片原始數(shù)據(jù)處理,參考limma的userguide。
1.幫助文檔里的示例
rm(list=ls())
library(limma)
x <- read.ilmn(files="probe profile.txt",
other.columns="Detection")
## Reading file probe profile.txt ... ...
x$E[1:4,1:4]
## 1 2 3 4
## ILMN_1762337 52.34406 46.10429 54.01238 47.65873
## ILMN_2055271 69.91481 73.86729 58.64280 72.36581
## ILMN_1736007 57.47208 53.70050 53.39117 49.43674
## ILMN_2383229 53.60817 57.50813 48.22960 48.18757
讀取原始文件秩伞，這么簡(jiǎn)單就拿到表達(dá)矩陣了逞带，原示例中有一個(gè)control probe profile.txt，文檔說(shuō)了它不是必須的纱新，我直接去掉咯掰担。畫(huà)個(gè)箱線圖
boxplot(log2(x$E),range=0,ylab="log2 intensity")
image.png
做背景校正、標(biāo)準(zhǔn)化
y <- neqc(x)
## Note: inferring mean and variance of negative control probe intensities from the detection p-values.
聽(tīng)起來(lái)很厲害怒炸，實(shí)際上一個(gè)函數(shù)搞定。
探針過(guò)濾
原始文件里的信息毡代，被拆分成了一個(gè)表達(dá)矩陣和一個(gè)Detection P值矩陣阅羹，P值可以用于過(guò)濾樣本。
示例數(shù)據(jù)總共有12個(gè)樣本教寂，過(guò)濾標(biāo)準(zhǔn)是至少在3個(gè)樣本里P值<0.05捏鱼。
x$other$Detection[1:4,1:4]
## 1 2 3 4
## ILMN_1762337 0.55849580 0.6754875 0.13698630 0.601388900
## ILMN_2055271 0.03064067 0.0000000 0.04931507 0.002777778
## ILMN_1736007 0.27715880 0.2924791 0.15342470 0.486111100
## ILMN_2383229 0.47353760 0.1866295 0.36575340 0.563888900
dim(y)
## [1] 48803 12
expressed <- rowSums(y$other$Detection < 0.05) >= 3 ;table(expressed)
## expressed
## FALSE TRUE
## 24112 24691
y <- y[expressed,]
dim(y)
## [1] 24691 12
這時(shí)的表達(dá)矩陣，也就能對(duì)接常規(guī)芯片數(shù)據(jù)的其他分析了酪耕，比如差異分析导梆。
exp = y$E
exp[1:4,1:4]
## 1 2 3 4
## ILMN_2055271 5.085517 5.294789 5.047373 5.274919
## ILMN_1653355 5.803398 5.901535 5.378863 5.567170
## ILMN_1787689 5.164395 4.929249 5.401202 5.006137
## ILMN_1745607 10.508010 9.922797 5.689654 5.317240
2.GEO數(shù)據(jù)實(shí)例
去GSE16997的頁(yè)面下載它的補(bǔ)充文件“GSE16997_raw.txt”。
讀取、背景校正和標(biāo)準(zhǔn)化
rm(list = ls())
x <- read.ilmn(files="GSE16997_raw.txt",
expr="Sample",
probeid="ID_REF",
other.columns="Detection Pval")
## Reading file GSE16997_raw.txt ... ...
y <- neqc(x,detection.p="Detection Pval")
## Note: inferring mean and variance of negative control probe intensities from the detection p-values.
探針過(guò)濾
x$other$Detection[1:4,1:4]
## 1 2 3 4
## ILMN_1762337 0.55849580 0.6754875 0.13698630 0.601388900
## ILMN_2055271 0.03064067 0.0000000 0.04931507 0.002777778
## ILMN_1736007 0.27715880 0.2924791 0.15342470 0.486111100
## ILMN_2383229 0.47353760 0.1866295 0.36575340 0.563888900
dim(y)
## [1] 48803 12
expressed <- rowSums(y$other$`Detection Pval` < 0.05) >= 3 ;table(expressed)
## expressed
## FALSE TRUE
## 24112 24691
y <- y[expressed,]
dim(y)
## [1] 24691 12
y$E[1:4,1:4]
## 1 2 3 4
## ILMN_2055271 5.085517 5.294789 5.047373 5.274919
## ILMN_1653355 5.803398 5.901535 5.378863 5.567170
## ILMN_1787689 5.164395 4.929249 5.401202 5.006137
## ILMN_1745607 10.508010 9.922797 5.689654 5.317240
