前往小程序,Get更优阅读体验!
立即前往
首页
学习
活动
专区
工具
TVP
发布
社区首页 >专栏 >aglient芯片原始数据处理

aglient芯片原始数据处理

作者头像
生信技能树
发布2020-06-11 15:47:40
3.5K1
发布2020-06-11 15:47:40
举报
文章被收录于专栏:生信技能树

导读

我多次在学徒作业强调了 3大基因芯片产商里面,就Agilent公司的芯片比较难搞,比如Agilent芯片表达矩阵处理(学徒作业) 以及 oligo包可以处理agilent芯片吗,这个作业难度非常高,不过我们生信技能树优秀讲师:小洁在繁重的授课压力下抽空整理了相关数据处理经验分享给大家,下面看她的表演:

本文讲的是aglient芯片原始数据的处理,参考资料是limma 的userguide文档。GEO数据库下载的表达矩阵不符合预期,比如是空的,或者是有负值的,那我们就处理一下它的原始数据。aglient的芯片应用也很广泛,举个OSCC的栗子:GSE23558,跟着学习学习。

1.下载和读取数据

1.1获取临床信息数据

从前,提到GEO数据下载,我们只有GEOquery,神功盖世,但是死于网速。后来就有了中国人寄几的GEO镜像,AnnoProbe包。还没有正式发表,就已经初露锋芒了,因简单易学,下载迅速,在我们的粉丝圈子里很受欢迎。

代码语言:javascript
复制
rm(list = ls())
library(stringr)
library(AnnoProbe)
library(GEOquery)
library(limma)
gse = "GSE23558"
geoChina(gse)
代码语言:javascript
复制
## you can also use getGEO from GEOquery, by 
## getGEO("GSE23558", destdir=".", AnnotGPL = F, getGPL = F)
代码语言:javascript
复制
## $GSE23558_series_matrix.txt.gz
## ExpressionSet (storageMode: lockedEnvironment)
## assayData: 41000 features, 32 samples 
##   element names: exprs 
## protocolData: none
## phenoData
##   sampleNames: GSM577914 GSM577915 ... GSM577945 (32 total)
##   varLabels: title geo_accession ... tissue:ch1 (42 total)
##   varMetadata: labelDescription
## featureData: none
## experimentData: use 'experimentData(object)'
##   pubMedIds: 22072328
## 28433800 
## Annotation: GPL6480

提供一个GSE编号就可以下载啦。因为表达矩阵是处理过的,我们不要,所以只提取临床信息表格,从中获得分组信息。

代码语言:javascript
复制
load("GSE23558_eSet.Rdata")
pd <- pData(gset[[1]])

调整pd的行名与文件读取的顺序一致,并定义分组信息。

代码语言:javascript
复制
raw_dir = "rawdata/GSE23558_RAW/"
raw_datas = paste0(raw_dir,"/",dir(raw_dir))

#调整pd与rawdata的顺序一致
raw_order = str_extract(raw_datas,"GSM\\d*")
pd = pd[match(raw_order,rownames(pd)),]

group_list <- ifelse(stringr::str_detect(pd$title,"HealthyControl"),"normal","tumor")
group_list <- factor(group_list, levels=c("normal","tumor"))
1.2 读取原始数据

这个原始数据下载,在GEO主页,可能对大家的网络有一点点要求,可以参考:下载GEO数据太慢?快用axel

代码语言:javascript
复制
x <- read.maimages(raw_datas,
                   source="agilent", 
                   green.only=TRUE,
                   other.columns = "gIsWellAboveBG")
代码语言:javascript
复制
## Read rawdata/GSE23558_RAW//GSM577914.txt 
## Read rawdata/GSE23558_RAW//GSM577915.txt 
## Read rawdata/GSE23558_RAW//GSM577916.txt 
## Read rawdata/GSE23558_RAW//GSM577917.txt 
## Read rawdata/GSE23558_RAW//GSM577918.txt 
## Read rawdata/GSE23558_RAW//GSM577919.txt 
## Read rawdata/GSE23558_RAW//GSM577920.txt 
## Read rawdata/GSE23558_RAW//GSM577921.txt 
## Read rawdata/GSE23558_RAW//GSM577922.txt 
## Read rawdata/GSE23558_RAW//GSM577923.txt 
## Read rawdata/GSE23558_RAW//GSM577924.txt 
## Read rawdata/GSE23558_RAW//GSM577925.txt 
## Read rawdata/GSE23558_RAW//GSM577926.txt 
## Read rawdata/GSE23558_RAW//GSM577927.txt 
## Read rawdata/GSE23558_RAW//GSM577928.txt 
## Read rawdata/GSE23558_RAW//GSM577929.txt 
## Read rawdata/GSE23558_RAW//GSM577930.txt 
## Read rawdata/GSE23558_RAW//GSM577931.txt 
## Read rawdata/GSE23558_RAW//GSM577932.txt 
## Read rawdata/GSE23558_RAW//GSM577933.txt 
## Read rawdata/GSE23558_RAW//GSM577934.txt 
## Read rawdata/GSE23558_RAW//GSM577935.txt 
## Read rawdata/GSE23558_RAW//GSM577936.txt 
## Read rawdata/GSE23558_RAW//GSM577937.txt 
## Read rawdata/GSE23558_RAW//GSM577938.txt 
## Read rawdata/GSE23558_RAW//GSM577939.txt 
## Read rawdata/GSE23558_RAW//GSM577940.txt 
## Read rawdata/GSE23558_RAW//GSM577941.txt 
## Read rawdata/GSE23558_RAW//GSM577942.txt 
## Read rawdata/GSE23558_RAW//GSM577943.txt 
## Read rawdata/GSE23558_RAW//GSM577944.txt 
## Read rawdata/GSE23558_RAW//GSM577945.txt
代码语言:javascript
复制
dim(x)
代码语言:javascript
复制
## [1] 45015    32
2.背景校正和标准化
代码语言:javascript
复制
y <- backgroundCorrect(x, method="normexp")
代码语言:javascript
复制
## Array 1 corrected
## Array 2 corrected
## Array 3 corrected
## Array 4 corrected
## Array 5 corrected
## Array 6 corrected
## Array 7 corrected
## Array 8 corrected
## Array 9 corrected
## Array 10 corrected
## Array 11 corrected
## Array 12 corrected
## Array 13 corrected
## Array 14 corrected
## Array 15 corrected
## Array 16 corrected
## Array 17 corrected
## Array 18 corrected
## Array 19 corrected
## Array 20 corrected
## Array 21 corrected
## Array 22 corrected
## Array 23 corrected
## Array 24 corrected
## Array 25 corrected
## Array 26 corrected
## Array 27 corrected
## Array 28 corrected
## Array 29 corrected
## Array 30 corrected
## Array 31 corrected
## Array 32 corrected
代码语言:javascript
复制
y <- normalizeBetweenArrays(y, method="quantile")
class(y)
代码语言:javascript
复制
## [1] "EList"
## attr(,"package")
## [1] "limma"

3. 基因过滤

  • 去除对照探针
  • 去除匹配不到genesymbol的探针
  • 去除不表达的探针,去除的标准是:至少在一半样本中高于背景,通过y(other)gIsWellAboveBG来判断。
  • 我自己加上了一个,测到多次的基因,只保留一个探针。
代码语言:javascript
复制
Control <- y$genes$ControlType==1L;table(Control)
代码语言:javascript
复制
## Control
## FALSE  TRUE 
## 43529  1486
代码语言:javascript
复制
NoSymbol <- is.na(y$genes$GeneName);table(NoSymbol)
代码语言:javascript
复制
## NoSymbol
## FALSE 
## 45015
代码语言:javascript
复制
IsExpr <- rowSums(y$other$gIsWellAboveBG > 0) >= 16;table(IsExpr)
代码语言:javascript
复制
## IsExpr
## FALSE  TRUE 
## 13088 31927
代码语言:javascript
复制
Isdup <- duplicated(y$genes$GeneName);table(Isdup)
代码语言:javascript
复制
## Isdup
## FALSE  TRUE 
## 30328 14687
代码语言:javascript
复制
yfilt <- y[!Control & !NoSymbol & IsExpr & !Isdup, ]
dim(yfilt)
代码语言:javascript
复制
## [1] 20650    32

可以看到,过滤后剩下了2万多个探针。

4.得到表达矩阵

代码语言:javascript
复制
exp = yfilt@.Data[[1]]
boxplot(exp)
代码语言:javascript
复制
exp[1:2,1:2]
代码语言:javascript
复制
##      rawdata/GSE23558_RAW//GSM577914 rawdata/GSE23558_RAW//GSM577915
## [1,]                        9.284154                       11.473334
## [2,]                        7.341236                        7.474406

得到的表达矩阵没问题,但行名和列名均有问题。行名应该是探针名,列名是样本名,调整一下。

4.1获得样本名
代码语言:javascript
复制
colnames(exp) = str_extract(colnames(exp),"GSM\\d*")
exp[1:2,1:2]
代码语言:javascript
复制
##      GSM577914 GSM577915
## [1,]  9.284154 11.473334
## [2,]  7.341236  7.474406
4.2 获得基因名

limma文档里写的是用了注释R包,在本例的原文件是里有探针注释的,这里直接使用。

可以直接将exp的行名转为基因名。行名不能重复,所以先去重

代码语言:javascript
复制
anno = yfilt$genes
nrow(anno);nrow(exp)
代码语言:javascript
复制
## [1] 20650
代码语言:javascript
复制
## [1] 20650
代码语言:javascript
复制
rownames(exp)=rownames(anno)
ids = unique(anno$GeneName)
exp = exp[!duplicated(anno$GeneName),]
rownames(exp) = ids
exp[1:4,1:4]
代码语言:javascript
复制
##             GSM577914 GSM577915 GSM577916 GSM577917
## APOBEC3B     9.284154 11.473334 10.439071 11.661000
## A_32_P77178  7.341236  7.474406  7.310818  7.397149
## ATP11B       9.963452  8.915621 10.193873  9.321954
## DNAJA1      13.469790 13.201078 12.827357 13.389431

至此,得到了标准的表达矩阵。后面要做什么就看你啦,这就相当于修复了一下数据库里那个被标准化过的表达矩阵。

5.差异分析

代码语言:javascript
复制
design <- model.matrix(~group_list)
fit <- lmFit(exp,design)
fit <- eBayes(fit,trend=TRUE,robust=TRUE)
summary(decideTests(fit))
代码语言:javascript
复制
##        (Intercept) group_listtumor
## Down             0            2102
## NotSig           0           16928
## Up           20650            1620
代码语言:javascript
复制
deg = topTable(fit,coef=2,n=dim(y)[1])
boxplot(exp[rownames(deg)[1],]~group_list)

这里直接走limma的简易流程,可以画差异最显著的那个基因表达量看看,可以看到差异是超级明显了!

代码语言:javascript
复制
save(exp,group_list,deg,gse,file = paste0(gse,"deg.Rdata"))
本文参与 腾讯云自媒体同步曝光计划,分享自微信公众号。
原始发表:2020-06-10,如有侵权请联系 cloudcommunity@tencent.com 删除

本文分享自 生信技能树 微信公众号,前往查看

如有侵权,请联系 cloudcommunity@tencent.com 删除。

本文参与 腾讯云自媒体同步曝光计划  ,欢迎热爱写作的你一起参与!

评论
登录后参与评论
0 条评论
热度
最新
推荐阅读
目录
  • 1.下载和读取数据
    • 1.1获取临床信息数据
      • 1.2 读取原始数据
        • 2.背景校正和标准化
        • 3. 基因过滤
        • 4.得到表达矩阵
          • 4.1获得样本名
            • 4.2 获得基因名
            • 5.差异分析
            相关产品与服务
            数据库
            云数据库为企业提供了完善的关系型数据库、非关系型数据库、分析型数据库和数据库生态工具。您可以通过产品选择和组合搭建,轻松实现高可靠、高可用性、高性能等数据库需求。云数据库服务也可大幅减少您的运维工作量,更专注于业务发展,让企业一站式享受数据上云及分布式架构的技术红利!
            领券
            问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档