Aglient的芯片在科研界也是一大寵兒，通常根據(jù)其染色分為單通道和多通道兩種。最為奇葩的是Aglient芯片的許多表達(dá)矩陣下載后發(fā)現(xiàn)有空值、負(fù)值，因此就要求我們從原始數(shù)據(jù)開始著手。下面就一起學(xué)習(xí)下吧。

核心函數(shù)：

read.maimages(raw_datas, source = "agilent", green.only = T, other.columns = "gIsWellAboveBG")

1.單通道芯片

以下以GSE23558為例，是《aglient芯片原始數(shù)據(jù)處理》的學(xué)習(xí)筆記。

1.1 數(shù)據(jù)下載及讀取

rm(list = ls())
library(tidyverse)
library(limma)
library(GEOquery)
library(AnnoProbe)
gse="GSE23558"
#setwd(gse)
#geoChina(gse)
load("D:/jianshu/microarry-analysis/GSE23558/GSE23558_eSet.Rdata") #提取原始數(shù)據(jù)

# 分組信息
pd <- pData(gset[[1]])
raw_dir <- "D:/jianshu/microarry-analysis/GSE23558/GSE23558_RAW"
raw_datas <- paste0(raw_dir,"/",list.files(raw_dir))
raw_order <- str_extract(raw_datas,"GSM\\d*")
pd <- pd[match(raw_order,rownames(pd)),]
pd <- pd %>% 
  select(geo_accession,`tissue:ch1`)
colnames(pd) <- c("id","type")
pd$type <- case_when(pd$type=="Oral Tumor"~"tumor",
                     T~"normal")
pd$type <- factor(pd$type,levels = c("normal","tumor"))
group_list <- pd$type
names(group_list) <- pd$id

#原始數(shù)據(jù)讀取
data.raw <- read.maimages(raw_datas,
                          source = "agilent",
                          green.only = T,
                          other.columns = "gIsWellAboveBG")

## Read D:/jianshu/microarry-analysis/GSE23558/GSE23558_RAW/GSM577914.txt.gz 
## Read D:/jianshu/microarry-analysis/GSE23558/GSE23558_RAW/GSM577915.txt.gz 
## Read D:/jianshu/microarry-analysis/GSE23558/GSE23558_RAW/GSM577916.txt.gz 
......
## Read D:/jianshu/microarry-analysis/GSE23558/GSE23558_RAW/GSM577943.txt.gz 
## Read D:/jianshu/microarry-analysis/GSE23558/GSE23558_RAW/GSM577944.txt.gz 
## Read D:/jianshu/microarry-analysis/GSE23558/GSE23558_RAW/GSM577945.txt.gz

1.2 背景矯正和標(biāo)準(zhǔn)化

data.bg <- backgroundCorrect(data.raw,method = "normexp")

## Array 1 corrected
## Array 2 corrected
## Array 3 corrected
.....
## Array 29 corrected
## Array 30 corrected
## Array 31 corrected
## Array 32 corrected

data.norm <- normalizeBetweenArrays(data.bg,method = "quantile")

1.3 基因過濾

去掉對照探針、未匹配到genesymbol探針、去表達(dá)探針（至少在一般樣本中高于背景）、重復(fù)探針。

ctrl <- data.norm$genes$ControlType==1L
Nosymbol <- is.na(data.norm$genes$GeneName)
IsExpr <- rowSums(data.norm$other$gIsWellAboveBG>0)>= nrow(pd)/2
Isdup <- duplicated(data.norm$genes$GeneName)
data.filt <- data.norm[!ctrl&!Nosymbol&IsExpr&!Isdup,]
dim(data.filt)

## [1] 20650    32

過濾后剩余2萬零650個(gè)探針。

1.4 表達(dá)矩陣

data.exp <- data.filt@.Data[[1]]
library(RColorBrewer)
colors <- brewer.pal(12,"Set3")
boxplot(data.exp,col=colors,las=3)

image.png

colnames(data.exp) <- str_extract(colnames(data.exp),"GSM\\d*")

1.5 獲得基因名

anno <- data.filt$genes
nrow(anno);nrow(data.exp)

## [1] 20650

## [1] 20650

rownames(data.exp)=anno$GeneName
ids <- unique(anno$GeneName)
data.exp <- data.exp[!duplicated(anno$GeneName),]

其實(shí)，整個(gè)過程相當(dāng)于對作者上傳的標(biāo)準(zhǔn)化矩陣進(jìn)行了修復(fù)。

1.6 差異分析

design <- model.matrix(~group_list)
fit <- lmFit(data.exp,design)
fit1 <- eBayes(fit,trend = T,robust=T)
summary(decideTests(fit))

##        (Intercept) group_listtumor
## Down             0            2090
## NotSig           0           16887
## Up           20650            1673

options(digits = 4)
deg <- topTable(fit1,coef = 2,n=dim(data.exp)[1])
boxplot(data.exp[rownames(deg)[1],]~group_list)

image.png

2.雙通道芯片

rm(list = ls())
gse="GSE29609"
library(limma)
library(AnnoProbe)
#setwd(gse)
#geoChina(gse)
load("D:/jianshu/microarry-analysis/GSE29609_eSet.Rdata")
pd <- Biobase::pData(gset[[1]])
raw_dir <- "D:/jianshu/microarry-analysis/GSE29609/GSE29609_RAW"

raw_datas <- paste0(raw_dir,"/",list.files(raw_dir,pattern = "GSM\\d*"))
raw_order <- str_extract(raw_datas,"GSM\\d*")
pd <- pd[match(raw_order,rownames(pd)),]

#原始數(shù)據(jù)讀取
data.raw <- read.maimages(raw_datas,
                          source = "agilent",
                          green.only = F,
                          other.columns = "gIsWellAboveBG")

## Read D:/jianshu/microarry-analysis/GSE29609/GSE29609_RAW/GSM733579_US22502565_251239115211_S01_A01_GE2_44k_1005.txt.gz 
## Read D:/jianshu/microarry-analysis/GSE29609/GSE29609_RAW/GSM733580_US22502565_251239125482_S01_A01_GE2_44k_1005.txt.gz 
## Read D:/jianshu/microarry-analysis/GSE29609/GSE29609_RAW/GSM733581_US22502565_251239144561_S01_A01_GE2_44k_1005.txt.gz 
......
## Read D:/jianshu/microarry-analysis/GSE29609/GSE29609_RAW/GSM733615_US22502565_251239125485_S01_A01_GE2_44k_1005.txt.gz 
## Read D:/jianshu/microarry-analysis/GSE29609/GSE29609_RAW/GSM733616_US22502565_251239115213_S01_A01_GE2_44k_1005.txt.gz 
## Read D:/jianshu/microarry-analysis/GSE29609/GSE29609_RAW/GSM733617_US22502565_251239144552_S01_A01_GE2_44k_1005.txt.gz

2.1 背景矯正、標(biāo)準(zhǔn)化

data.bg <- backgroundCorrect(data.raw,method = "normexp")

## Array 1 corrected
## Array 2 corrected
## Array 3 corrected
......
## Array 37 corrected
## Array 38 corrected
## Array 39 corrected
## Array 1 corrected
## Array 2 corrected
## Array 3 corrected
......
## Array 37 corrected
## Array 38 corrected
## Array 39 corrected

data.norm <- normalizeBetweenArrays(data.bg,method = "quantile")

ctrl <- data.norm$genes$ControlType==1L
Nosymbol <- is.na(data.norm$genes$GeneName)
#IsExpr <- rowSums(data.norm$other$gIsWellAboveBG>0)>= nrow(pd)/2
Isdup <- duplicated(data.norm$genes$GeneName)
data.filt <- data.norm[!ctrl&!Nosymbol&!Isdup,]
dim(data.filt)

## [1] 31036    39

data.exp <- data.filt@.Data[[4]]
library(RColorBrewer)
colors <- brewer.pal(12,"Set3")
boxplot(data.exp,col=colors,las=3)

image.png

疑問：雙通道芯片，我是按照單通道的處理的，不知道是否正確？希望得到大神的指導(dǎo)，萬分感謝。

參考鏈接：

aglient芯片原始數(shù)據(jù)處理

雙通道芯片MA和density