前言
27K的数据是很老的芯片数据，但是客户有需求就要找方法分析，主流的DNA甲基化芯片R包minfi和champ都只支持450K和850K的芯片。所以在bioconductor中搜索到了methylumi这个包，可以从idat读数据，经过质控得到beta值矩阵，之后用limma做差异分析。
可以参考这篇文章[https://www.ncbi.nlm.nih.gov/pubmed/32375395]
450K和850K的芯片分析就简单多了，生信技能树也写好了pipeline可以参考[https://github.com/jmzeng1314/methy_array]
从TCGA下载
library(TCGAbiolinks)
m1 <- GDCquery(
project = "TCGA-KIRC",
data.category = "DNA Methylation",
data.type = "Masked Intensities",
data.format = "idat",
platform = "Illumina Human Methylation 27"
)
GDCdownload(m1)
移动文件到KIRC_27K_idat文件夹下
mv GDCdata/TCGA-KIRC/harmonized/DNA_Methylation/Masked_Intensities/*idat KIRC_27K_idat
制作重命名shell脚本
metadata.cart.2023-02-09.json这个文件是在TCGA官网选中样本和数据类型后下载的样本信息，里边包含了样本名和文件名。
library(jsonlite)
library(magrittr)
library(data.table)
j=jsonlite::read_json('metadata.cart.2023-02-09.json')
tb <- map_dfr(j,~tibble(file_name=.x$file_name,new_name=.x$associated_entities[[1]]$entity_submitter_id))
tb %<>% arrange(new_name)
colors = str_split(tb$file_name,'_',simplify = T)[,3]
n=rep(1:(826/2),2) %>% sort()
ns = sprintf("%03d",n)
sample_name = tb$new_name
tb$new_name <- paste0(tb$new_name,'_','R',ns,'C',ns,'_',colors)
tb$mv = "mv"
fwrite(tb[,c(3,1,2)],"KIRC_27K_idat/rename.sh",sep=' ',col.names=F )
重命名
cd KIRC_27K_idat
bash rename.sh
重命名脚本的前几行
mv 348750ad-930b-4a62-98fc-165a8216cf42_noid_Grn.idat TCGA-A3-3306-01A-01D-0859-05_R001C001_Grn.idat
mv 348750ad-930b-4a62-98fc-165a8216cf42_noid_Red.idat TCGA-A3-3306-01A-01D-0859-05_R001C001_Red.idat
mv 63c1410d-9b54-47a8-bb8f-08a030dacab0_noid_Grn.idat TCGA-A3-3306-11A-01D-0859-05_R002C002_Grn.idat
mv 63c1410d-9b54-47a8-bb8f-08a030dacab0_noid_Red.idat TCGA-A3-3306-11A-01D-0859-05_R002C002_Red.idat
methylumi读数据
idatPath必须是绝对路径
idatPath <- "~/Project/20230203_DNAmeth/data/KIRC_27K_idat"
mset27k <- methylumIDAT(getBarcodes(path=idatPath), idatPath=idatPath)
sampleNames(mset27k) <- unique(sample_name)
标准化前的质控图
qc.probe.plot(mset27k, by.type=TRUE)
标准化后的质控图
mset27k_pp <- stripOOB(normalizeMethyLumiSet(methylumi.bgcorr(mset27k)))
qc.probe.plot(mset27k_pp, by.type=TRUE)
limma差异分析
跟mRNA芯片的差异分析一样，最后的deg_df可以用于绘制火山图
beta=betas(mset27k_pp)
beta %<>% as.data.frame() %>% dplyr::select(matches('^.{13}[01]1A'))
group = ifelse(str_detect(colnames(beta),'^.{13}01'),'Tumor','Normal')
design <- model.matrix(~ 0 + factor(group))
colnames(design) <- levels(factor(group))
rownames(design) <- colnames(beta)
contrasts <- paste0("Tumor", "-", "Normal")
contrast.matrix <- makeContrasts(contrasts = contrasts, levels = design)
fit <- lmFit(beta, design)
fit <- contrasts.fit(fit, contrast.matrix)
fit <- eBayes(fit, 0.01)
deg_df <- topTable(fit, adjust = "fdr", sort.by = "B", number = nrow(beta)) %>% na.omit()
Reference
https://www.bioconductor.org/packages/release/bioc/vignettes/methylumi/inst/doc/methylumi.pdf
https://www.bioconductor.org/packages/release/bioc/manuals/methylumi/man/methylumi.pdf
https://www.bioconductor.org/packages/release/bioc/vignettes/methylumi/inst/doc/methylumi450k.pdf
https://mp.weixin.qq.com/s/BIxtWJAO8AXHbNDItS7PFQ
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC7291297/
网友评论