在step1中，我們獲得了TCGA中OSCC 的32個(gè)病人的T-N配對樣本和對應(yīng)的臨床信息,并將其組成了一個(gè)名為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 歸一化

做后續(xù)差異分析之前，需要對信號(hào)值矩陣進(jìn)行歸一化。這一步驟消耗計(jì)算資源較多，配置不夠可能會(huì)跑很久或者會(huì)中斷。

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)

# 歸一化過程產(chǎn)生了缺失值,需要將有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個(gè)（29對）樣本。

2.2 數(shù)據(jù)質(zhì)控三張圖

# 主成分分析
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

# 相關(guān)關(guān)系矩陣熱圖
pheatmap::pheatmap(cor(myNorm[cg,]),
                   annotation_col = ac,
                   show_rownames = F,
                   show_colnames = F)

2.3 剔除聚類失敗的樣本

圖中看出三個(gè)樣本異常,刪掉它們和它們的配對樣本。

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")

根據(jù)top1000sd的熱圖和相關(guān)性熱圖，會(huì)發(fā)現(xiàn)三個(gè)樣本是異常的，因此又剔除3對，剩下26對（52個(gè)）樣本，用于下一步的差異分析。我試了一下這三個(gè)樣本不刪除的話，后面做差異甲基化位點(diǎn)的熱圖也是聚類不成功的，刪掉會(huì)好些。