最近从geo上下了一个比较大的甲基化数据集(GSE145361),idat文件一共有3778个(grn和red一一对应),本想直接利用构建好的pd文件和idat原始文件直接读入ChAMP,然后一套流程走下来,万万没想到内存超过限制,太可怕了,128G的内存再加上19G的虚拟内存,竟然都不够使,无奈之下只好下载处理好的beta矩阵进行运算,以前都是直接拿原始数据运算,我竟然不知道怎么导入beta矩阵,现在把它记下来:
1.第一步先读入R(由于getGEO下载受限,手动下载的)
beta <- readr::read_delim("./methylation/GSE145361_Vallerga2020_NCOMMS_AvgBeta_Matrix-file.txt", col_names = T, delim = "\t")
运算还可以,比R自带的read.table()要快
2.导入和过滤:由于是直接利用beta矩阵,就不需要用champ.load()或champ.import()导入了,直接利用champ.filter()进行过滤
pd <- readr::read_csv("./methylation/sample_info.csv")
myload <- champ.filter(beta=beta,pd=pd,fixOutlier = F,filterBeads = F,filterDetP = F,autoimpute = F)
ps:后面三个可以不用改T为F,但是fixOutlier一定要改为F,不然始终报错,查看源代码发现与原代码中的which使用有关,我弄清代码之后改为自己手动过滤了
beta <- myload$beta
beta <- na.omit(beta)
beta[beta<=0] <- min(beta[beta>0])
beta[beta>=1] <- max(beta[beta<1])
myload$beta <- beta
3.标准化:beta矩阵只有PBC和BMIQ两种标准化方法,建议BMIQ
时间非常长,耐心等待
myNorm <- champ.norm(beta = myfilter$beta)
翻看以前的学习笔记,发现minfi包也可以直接导入beta矩阵
一般情况下,我们是使用read.metharray.sheet()配合read.metharray.exp()进行导入,没有原始数据的话可以用makeGenomicRatioSetFromMatrix(),这两个包太贴心了~