搜索
查看: 1103|回复: 1

[basic] R处理芯片数据

[复制链接]

29

主题

29

帖子

149

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
149
发表于 2018-9-20 18:13:52 | 显示全部楼层 |阅读模式
1) 熟悉CEL file

从 NCBI GEO (http://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE24460)下载GSE24460. 将得到一个 GSE24460_RAW.tar
[Python] 纯文本查看 复制代码
if("affy" %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");biocLite("affy")}
suppressMessages(library(affy))
ls('package:affy')
myData <- ReadAffy(celfile.path="F:/R/bioinformation")
#myData1 <- ReadAffy(filenames = "F:/R/bioinformation with R cookbook/GSM602658_MCF71.CEL") 

--------如果不是从CEL文件读取,而是有多个独立的临床、实验、表达矩阵等文件,则需根据这些文件构建新的ExpressionSet对象---------------

每一部分组合成 ExpressionSet 对象,都扮演各自的角色。 exprs object 是表达量, phenotypic data 是样本临床信息 ( sex, age, treatment ...), annotated package 提供基本数据操作工具 。
[Python] 纯文本查看 复制代码
###################构建ExpressionSet对象(包含临床、实验、表达矩阵等多种信息)#####################
######利用自带数据集演示#######
suppressMessages(library(Biobase))
DIR <- system.file("extdata", package="Biobase")
exprsLoc <- file.path(DIR, "exprsData.txt")
pDataLoc <- file.path(DIR, "pData.txt")
exprs <- as.matrix(read.csv(exprsLoc, header = TRUE, sep = "\t", row.names = 1, as.is = TRUE))#读取表达矩阵
class(exprs)
dim(exprs)
pData <- read.table(pDataLoc, row.names = 1, header = TRUE, sep = "\t") #读取临床信息
pData <- new("AnnotatedDataFrame", data = pData)   #构建pData对象
exData <- new("MIAME", name="ABCabc", lab="XYZ Lab", contact="abc@xyz", title="", abstract="", url="www.xyz") #编译实验信息
exampleSet <- new("ExpressionSet", exprs = exprs, phenoData = pData, experimentData = exData, annotation = "hgu133a2")#构建ExpressionSet对象
str(exampleSet)
validObject(exampleSet) #检验构建的ExpressionSet对象的有效性

2)了解affBatch对象结构
[Python] 纯文本查看 复制代码
myData
str(myData)
pData(myData)   #临床信息
phenoData(myData)
exprs(myData) #获取表达矩阵
annotation(myData)# 获取注释信息
probeNames(myData) #获取探针名称
sampleNames(myData) #获取样本名称

3)Checking the quality of data(质控)

与质量相关的问题可能源于:
1、杂交,因为芯片上的荧光不均匀会导致不同的强度分布,非特异性结合或其他生物/技术原因可能在数据中产生背景噪声。
2、不合适的实验设计可能会影响整个数据集。 使用此类数据将导致数据分析期间的错误或不确定的推断。
因此,必须在开始数据分析之前确保数据质量。 这是通过寻找边远数组,数组内的分布,批处理效果等来实现。 有各种分析和诊断图可用于计算
[Python] 纯文本查看 复制代码
if("arrayQualityMetrics" %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");biocLite("arrayQualityMetrics")}
suppressMessages(library(arrayQualityMetrics))
if("arrayQualityMetrics" %in% rownames(installed.packages()) == FALSE) {source("http://bioconductor.org/biocLite.R");biocLite("arrayQualityMetrics")}
suppressMessages(library(arrayQualityMetrics))
arrayQualityMetrics(myData, outdir="microarray")  #质控
browseURL(file.path("microarray", "index.html"))
MAplot(myData, pairs=TRUE, plot.method="smoothScatter") #MAplot图
plotDensity.AffyBatch(myData) #密度图
boxplot(myData)               #箱型图
rnaDeg <- AffyRNAdeg(myData)  #查看RNA降解
plotAffyRNAdeg(rnaDeg)
summaryAffyRNAdeg(rnaDeg) #获取RNA降解情况

4)Data 标准化
标准化用来降低技术影响,产生可比较数据。因为有多种方法可以标准化数据,这里讲vsn, loess, quantile三种方法。
[Python] 纯文本查看 复制代码
library(vsn)
ls('package:vsn')
myData <- ReadAffy(celfile.path="F:/R/bioinformation with R cookbook") #读取所有CEL文件,返回的是AffyBatch object
#myData.VSN <- normalize.AffyBatch.vsn(myData)  #vsn包提示找不到该函数
myData.loess <- normalize.AffyBatch.loess(myData)
boxplot(myData.loess)
myData.quantile <- normalize.AffyBatch.quantiles(myData)
boxplot(myData.quantile)



5)移除批次效应

[Python] 纯文本查看 复制代码
source("http://bioconductor.org/biocLite.R")
biocLite("bladderbatch")
library(sva) # contains batch removing utilities
library(bladderbatch) # The data to be used
data(bladderdata)
pheno <- pData(bladderEset) #提取表达矩阵 and pheno data
edata <- exprs(bladderEset) #提取样本临床信息
pheno
myData <- bladderEset[,sampleNames(bladderEset)[1:8]]  #提取同样本不同批次子集
arrayQualityMetrics(myData, outdir="qc_be")  #质控
mod1 <- model.matrix(~as.factor(cancer), data=pData(myData))[,c(1,3)]  #model矩阵
batch <- pData(myData)$batch  #样本批次信息
edata <- exprs(myData)        #表达矩阵
combat_edata <- ComBat(dat=edata, batch=batch, mod=mod1, par.prior=TRUE)#经验贝叶斯方法去merge批次效应
myData2 <- myData
exprs(myData2) <- combat_edata
arrayQualityMetrics(myData2, outdir="qc_nbe")


6) PCA
[Python] 纯文本查看 复制代码
myData <- ReadAffy(celfile.path="F:/R/bioinformation with R cookbook") #读取所有CEL文件,返回的是AffyBatch object
myData.pca <- exprs(myData)
myPca <- prcomp(myData.pca, scale=TRUE) #prcomp函数计算主成分
summary(myPca)
colors <- c("green","cyan","violet","magenta") 
pairs(myPca$x, col=colors)

7)差异表达基因
[Python] 纯文本查看 复制代码
library(affy) # Package for affy data handling
library(antiProfilesData) # Package containing input data
library(affyPLM) # Normalization package for eSet
library(limma) # limma analysis package
data(apColonData)
myData <- apColonData[, sampleNames(apColonData)[1:16]]
myData_quantile <- normalize.ExpressionSet.quantiles(myData)
design <- model.matrix(~0 + pData(myData)$Status)
fit <- lmFit(myData_quantile,design)
fit
fitE <- eBayes(fit)
tested <- topTable(fitE, adjust="fdr", sort.by="B", number=Inf)
DE <- x[tested$adj.P.Val<0.01,]
dim[DE]
DE <- x[tested$adj.P.Val< 0.01 & abs(x$logFC) >2,]






本帖子中包含更多资源

您需要 登录 才可以下载或查看,没有帐号?立即注册

x
回复

使用道具 举报

0

主题

1

帖子

43

积分

新手上路

Rank: 1

积分
43
发表于 2019-2-12 22:32:29 | 显示全部楼层
请问一下:
7)差异表达基因
第13行代码x是代表什么?前面一直没有出现x,运行到这一步就出错了
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-4-26 10:44 , Processed in 0.031623 second(s), 24 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.