R包minfi处理DNA甲基化芯片数据

minfi：Analyze Illumina Infinium DNA methylation arrays

1. 数据读入

if (!require("BiocManager", quietly = TRUE))install.packages("BiocManager")BiocManager::install("minfi",force = TRUE)
BiocManager::install("minfiData")library(minfi)library(minfiData)### 1. 数据读入
# RGsetEx
# ## RGsetEx: RGChannelSet, 622,399 features
# class(RGsetEx) #"RGChannelSet"# MsetEx <- preprocessRaw(RGsetEx)
# class(MsetEx)  # "MethylSet"
# ## MsetEx: MethylSet, 485,512 features# GMsetEx <- mapToGenome(MsetEx)
# class(GMsetEx) #"GenomicMethylSet"
# ## GMsetEx: GenomicMethylSet, 485,512 features## ----baseDir------------------------------------------------------------------
baseDir <- system.file("extdata", package = "minfiData")
list.files(baseDir)
list.files(file.path(baseDir, "5723646052"))
## 读Illumina methylation sample sheet(csv结尾的文件),含有样本信息以及文件位置
targets <- read.metharray.sheet(baseDir) # targets 为data.frame# Usage
# read.metharray.sheet(base, pattern = "csv$", ignore.case = TRUE,
#                      recursive = TRUE, verbose = TRUE)# sub(baseDir, "", targets$Basename) # 替换
## 通过样本信息（targets为data.frame）读取甲基化芯片数据idat
RGset <- read.metharray.exp(targets = targets)
class(RGset) # "RGChannelSet"
# 只读取某一文件夹下面的idat文件
RGset2 <- read.metharray.exp(file.path(baseDir, "5723646052"))
# 读取baseDir下所有子文件夹下面的idat文件
RGset3 <- read.metharray.exp(baseDir, recursive = TRUE)

2. 质控

### 2. 质控
## 甲基化β值的豆状密度图。
par(mar=c(5,6,4,2))  # c(bottom, left, top, right)
densityBeanPlot(RGset, sampNames=names, sampGroups=group)
## β值密度曲线
densityPlot(RGset,sampGroups=group)# 检测失败的位点
detP<- detectionP(RGset, type = "m+u")
failed <- detP>0.01
colMeans(failed) # Fraction of failed positions per sample
sum(rowMeans(failed)>0.5) # How many positions failed in >50% of samples

3.数据类型转化

### 3.数据类型转化
# 产生没有标准的MethylSet类型
Mset <- preprocessRaw(RGset)
# quantile normalization
gmSet <- preprocessQuantile(Mset)# GMset <- mapToGenome(Mset)
# getMeth(GMset)
# getUnmeth(GMset)## 对于RGset，Mset，GMset， gmSet对象 都有这些函数
pd <- pData(gmSet) # 表型数据
getBeta(gmSet)
group <- pd$Sample_Group # 样品分组
annotation(gmSet) # 注释信息
getAnnotatio(gmSet)
ids <- rownames(gmSet)  # 探针id名
names <- colnames(gmSet)  # 样品名
beta <- getBeta(gmSet) #β值
#assays(RGset)[[1]]  # assays(RGset)$Green
#assays(RGset)[2] # assays(RGset)$Red

4. 读入TCGA数据和GEO数据

### 4. 读入TCGA数据和GEO数据
data_file <- "/Users/zhengxueming/test/test_READ.methylation450.tsv"
gmSet_tcga <- readTCGA(data_file)
class(gmSet_tcga) #"GenomicRatioSet"
?readGEORawFile
?read.metharray
?read.metharray.exp
?read.metharray.sheet

5.甲基化差异化分析

### 5.甲基化差异化分析
## 检测差异甲基化位点
# 要用normalization后的gmSet数据
pd <- pData(gmSet)
beta <- getBeta(gmSet)
group <- pd$Sample_Group
dmp <- dmpFinder(beta, pheno = group, type = "categorical")
class(dump) # "data.frame"
# Usage
# dmpFinder(dat, pheno, type = c("categorical", "continuous"),
#           qCutoff = 1, shrinkVar = FALSE)
# Continuous phenotypes are tested with linear regression,
# while an F-test is used for categorical phenotypes.head(dmp)
## 选择统计学差异的位点
# sum(dmp$qval < 0.05, na.rm=TRUE)
dmp_sig <- dmp[which(dmp$qval < 0.05),]
head(dmp_sig)##检测差异甲基化区域
designMatrix <- model.matrix(~ gmSet$status)
#运行时间较长，并行处理
require(doParallel)
registerDoParallel(cores = 2)
bumps <- bumphunter(gmSet, design = designMatrix, B = 10,type = "Beta", cutoff = 0.25)class(bumps) #"list"
names(bumps)
dim(bumps$table)
head(bumps$coef)
length(bumps$coef)
head(bumps$pvaluesMarginal)
bumps$pvaluesMarginal
head(bumps$fitted)bumps_table <- bumps$table
dim(bumps_table)
head(bumps_table)
summary(bumps_table$fwer)
bumps_table_sig <- bumps_table[which(bumps_table$p.value < 0.01),]
head(bumps_table_sig)
dim(bumps_table_sig)bumps_table_sig[which(bumps_table_sig$chr=="chr7"),]# dat <- dummyData()
# # Enable parallelization
# require(doParallel)
# registerDoParallel(cores = 2)
# # Find bumps
# dmrs <- bumphunter(dat$mat, design=dat$design, chr=dat$chr, pos=dat$pos,
#                     cluster=dat$cluster, coef=2, cutoff= 0.28, nullMethod="bootstrap",
#                     smooth=TRUE, B=250, verbose=TRUE,
#                     smoothFunction=loessByCluster)
# head(dmrs$table[,1:4], n =3)

参考：

https://www.bioconductor.org/packages/release/bioc/manuals/minfi/man/minfi.pdf

R包minfi处理DNA甲基化芯片数据相关推荐

国家生物信息中心开发DNA甲基化芯片数据标准化方法—GMQN
过去十年来,由于DNA甲基化芯片技术的不断发展以及测序成本的快速下降,DNA甲基化芯片数据呈现爆发式增长.这些数据是表观基因组关联研究(Epigenome-Wide Association Studi ...
r语言怎么把txt数据变成一个Rdata格式_甲基化芯片数据下载如何读入到R里面
数据是一切的开始,前面我们介绍了一些背景知识,主要是理解什么是DNA甲基化,为什么要检测它,以及芯片和测序两个方向的DNA甲基化检测技术.具体介绍在:甲基化的一些基础知识,也了解了甲基化芯片的一般分析 ...
DNA甲基化测序数据的分析流程及相关软件总结
目前检测DNA甲基化的方法众多,主要可以分为以下几类(如表1所示): 图片来源(凡时财等,中国科学: 生命科学,2015) <更多精彩,可关注微信公众号:AIPuFuBio,和大型免费综合生物信 ...
甲基化系列 2. 甲基化芯片数据介绍与下载（GEO）
点击关注,桓峰基因桓峰基因的教程不但教您怎么使用,还会定期分析一些相关的文章,学会教程只是基础,但是如果把分析结果整合到文章里面才是目的,觉得我们这些教程还不错,并且您按照我们的教程分析出来不错的结 ...
R包cgdsr下载MSKCC癌症基因组数据
#install.packages("cgdsr",repos="http://cran.us.r-project.org") library(cgdsr) # ...
芯片数据分析笔记【05】 | 处理芯片数据的R包
芯片数据分析笔记[01] | 基因芯片的基本原理芯片数据分析笔记[02] | 芯片数据库芯片数据分析笔记[03] | GEO数据库使用教程及在线数据分析工具芯片数据分析笔记[04] | Arra ...
keil551的芯片包不能用_r语言中使用Bioconductor 分析芯片数据
原文链接: r语言中使用Bioconductor 分析芯片数据tecdat.cn 介绍芯片数据分析流程有些复杂,但使用 R 和 Bioconductor 包进行分析就简单多了.本教程将一步一步的展 ...
amt630a芯片中文资料_甲基化芯片学习记录
今天是生信星球陪你的第539天大神一句话,菜鸟跑半年.我不是大神,但我可以缩短你走弯路的半年~ 就像歌儿唱的那样,如果你不知道该往哪儿走,就留在这学点生信好不好~ 这里有豆豆和花花的学习历程,从新 ...
DNA甲基化数据分析专题
欢迎关注"生信修炼手册"! DNA 甲基化作为重要的表观遗传学的标记,能够调控基因的表达,在生长发育和疾病相关研究领域都有着重要意义.测定甲基化的手段有很多,芯片作为一种成熟的手段 ...
DNA甲基化与表观遗传学数据挖掘与分析学习会（10月26-27 上海）
两天一夜高强度训练,理论与实战相结合赠送高清视频供学员复习使用会务背景表观遗传学调控多种生命活动及疾病发生发展,已成为近年来生命科学领域的研究热点,促使了生物学多领域研究的突破性进展.随着实验和 ...

R包minfi处理DNA甲基化芯片数据

minfi：Analyze Illumina Infinium DNA methylation arrays

1. 数据读入

2. 质控

3.数据类型转化

4. 读入TCGA数据和GEO数据

5.甲基化差异化分析

R包minfi处理DNA甲基化芯片数据相关推荐

最新文章

热门文章