TCGA甲基化芯片数据质控和过滤

在step1中，我们获得了TCGA中OSCC 的32个病人的T-N配对样本和对应的临床信息,并将其组成了一个名为my_Load的ChAMP对象。

rm(list = ls())
library(ChAMP)
library(stringr)
load('./Rdata/step1_myLoad.Rdata')
myLoad$beta[1:4,1:4]
#>            TCGA-CV-6959-01 TCGA-CV-6436-11 TCGA-CV-5966-11 TCGA-CV-6939-11
#> cg00000029       0.2999564       0.2905091       0.3686591       0.3632641
#> cg00000108       0.4633994       0.4575858       0.4726964       0.4645730
#> cg00000109       0.4633994       0.4575858       0.4726964       0.4645730
#> cg00000165       0.6190992       0.1511787       0.2370243       0.1444160
myLoad$pd[1:4,1:4]
#>           sampleID anatomic_neoplasm_subdivision      patient group_list
#> 1: TCGA-CV-6959-01                   Oral Tongue TCGA-CV-6959      Tumor
#> 2: TCGA-CV-6436-11                   Oral Tongue TCGA-CV-6436     Normal
#> 3: TCGA-CV-5966-11                   Oral Cavity TCGA-CV-5966     Normal
#> 4: TCGA-CV-6939-11                   Oral Tongue TCGA-CV-6939     Normal

2.样本和探针过滤

2.1 归一化

做后续差异分析之前，需要对信号值矩阵进行归一化。这一步骤消耗计算资源较多，配置不够可能会跑很久或者会中断。

norm_file = "./Rdata/step2_champ_myNorm.Rdata"
if(!file.exists(norm_file)){
  myNorm <- champ.norm(beta=myLoad$beta,arraytype="450K",cores=8)
  save(myNorm,file = norm_file)
}
load(norm_file)

# 归一化过程产生了缺失值,需要将有NA的样本和它们的配对样本一起删掉
num.na <- apply(myNorm,2,function(x)(sum(is.na(x))))
table(num.na)
#> num.na
#>      0 258616 260092 264579 
#>     61      1      1      1
names(num.na) = colnames(myNorm)
dt = names(num.na[num.na>0])
dn = str_replace(dt,"-01","-11")
keep = setdiff(colnames(myNorm),c(dt,dn))
myNorm = myNorm[,keep]
pd = myLoad$pd
pd = pd[pd$sampleID %in% keep,]
identical(pd$sampleID,colnames(myNorm))
#> [1] TRUE

删除缺失值样本后，还剩58个（29对）样本。

2.2 数据质控三张图

# 主成分分析
library(FactoMineR)
library(factoextra) 
dat <- t(myNorm)

group_list=pd$group_list
table(group_list)
#> group_list
#> Normal  Tumor 
#>     29     29

dat.pca <- PCA(dat, graph = FALSE) 
fviz_pca_ind(dat.pca,
             geom.ind = "point", 
             col.ind = group_list, 
             addEllipses = TRUE, 
             legend.title = "Groups")

# 热图
cg=names(tail(sort(apply(myNorm,1,sd)),1000))
library(pheatmap)
ac=data.frame(group=group_list)
rownames(ac)=colnames(myNorm)  
pheatmap(myNorm[cg,],show_colnames =F,show_rownames = F,
         annotation_col=ac)

dev.off()
#> null device 
#>           1

# 相关关系矩阵热图
pheatmap::pheatmap(cor(myNorm[cg,]),
                   annotation_col = ac,
                   show_rownames = F,
                   show_colnames = F)

2.3 剔除聚类失败的样本

图中看出三个样本异常,删掉它们和它们的配对样本。

pn = c("TCGA-CV-5971-01","TCGA-CV-6953-11","TCGA-CV-6955-11")
drop = str_sub(colnames(myNorm),1,12) %in% str_sub(pn,1,12)
table(drop)
#> drop
#> FALSE  TRUE 
#>    52     6
myNorm = myNorm[,!drop]
dim(myNorm)
#> [1] 412481     52

pd = pd[!(pd$patient %in% str_sub(pn,1,12)),]
identical(pd$sampleID,colnames(myNorm))
#> [1] TRUE

save(pd,myNorm,file = "./Rdata/step2_filtered_pd_myNorm.Rdata")

根据top1000sd的热图和相关性热图，会发现三个样本是异常的，因此又剔除3对，剩下26对（52个）样本，用于下一步的差异分析。我试了一下这三个样本不删除的话，后面做差异甲基化位点的热图也是聚类不成功的，删掉会好些。