0.介绍
本文的数据出自2017年的一篇文章：Seven-CpG-based prognostic signature coupled with gene expression predicts survival of oral squamous cell carcinoma。参考链接：TCGA数据库的癌症甲基化芯片数据重分析
这也是B站的生信技能树甲基化芯片课程示例数据。经过我的研究和整理，成了今天的推文,这也将作为我以后开新课的重要基础，看到就是赚到~
前面写的甲基化学习记录，有几张思维导图：
甲基化学习记录
因为当时电脑配置跑不了甲基化数据，代码就没有学成，中途放弃。现在找补回来！
1.输入数据和R包
从UCSC的xena网站上可以找到TCGA-HNSC的甲基化数据和病人的临床信息。将其存放于raw_data.（这里使用的是GDC-TCGA下载的信号值矩阵和从TCGA下载的临床信息和生存信息），信号值矩阵很大，下下来才发现GDC里的临床信息是没有完整的site信息的，遂，转战TCGA。（实在不想再次下载信号值矩阵，将就看）
image.png
library(data.table)
library(impute)
library(ChAMP)
library(stringr)
library(tibble)
options(stringsAsFactors = F)
if(!dir.exists("raw_data"))dir.create("raw_data")
if(!dir.exists("Rdata"))dir.create("Rdata")
if(!dir.exists("figure"))dir.create("figure")
1.1 临床信息表格
TCGA里没有专门的OSCC，而是合并为HNSC了。所以需要整个下载下来，然后筛选OSCC对应的样本。
这里就需要用到一个很重要的符号：%in%。点击充电%in%很简单
pd <- fread("./raw_data/HNSC_clinicalMatrix")
#colnames(pd)
table(pd$anatomic_neoplasm_subdivision)
#>
#> Alveolar Ridge Base of tongue Buccal Mucosa Floor of mouth Hard Palate
#> 18 30 23 69 8
#> Hypopharynx Larynx Lip Oral Cavity Oral Tongue
#> 10 140 3 89 159
#> Oropharynx Tonsil
#> 9 46
oscc <- c("Oral Cavity","Oral Tongue","Buccal Mucosa","Lip","Alveolar Ridge","Hard Palate","Floor of mouth")
table(pd$anatomic_neoplasm_subdivision %in% oscc)
#>
#> FALSE TRUE
#> 235 369
dim(pd);pd <- pd[pd$anatomic_neoplasm_subdivision %in% oscc,];dim(pd)
#> [1] 604 132
#> [1] 369 132
可见，604个样本中有369个是属于OSCC的。已经挑选出了这些样本对应的临床信息。进行简化和整理：
1.样本ID保留至15位，为了和甲基化信号值矩阵的样本ID保持一致 2.完整的pd有130多列，从中挑出有用的列 3.生成group_list(Tumor-normal分组)和patient(病人ID)列 4.选出tumor-normal配对的样本
pd <- pd[,c("sampleID","anatomic_neoplasm_subdivision")]
pd$sampleID = str_sub(pd$sampleID,1,15)
pd$patient = str_sub(pd$sampleID,1,12)
pd$group_list = ifelse(as.numeric(str_sub(pd$sampleID,14,15))<10,"Tumor","Normal")
table(pd$group_list)
#>
#> Normal Tumor
#> 48 321
tp = pd$patient[pd$group_list =="Normal"]
np = pd$patient[pd$group_list =="Tumor"]
okp = intersect(tp,np)
pd$patient <- str_sub(pd$sampleID,1,12)
pd <- pd[pd$patient %in% okp,]
rownames(pd) <- pd$sampleID
dim(pd)
#> [1] 96 4
到此，得到了48对病人的临床信息。接下来要找他们对应的甲基化信号值数据。
1.2 信号值矩阵
读取进来，并筛选有配对样本甲基化数据的病人。
a = data.table::fread("raw_data/TCGA-HNSC.methylation450.tsv.gz", data.table = F)
a[1:4,1:4]
#> Composite Element REF TCGA-CQ-6227-01A TCGA-H7-7774-01A TCGA-CV-6943-01A
#> 1 cg00000029 0.2533996 0.5590821 0.2670461
#> 2 cg00000108 NA NA NA
#> 3 cg00000109 NA NA NA
#> 4 cg00000165 0.4846212 0.6797669 0.4371214
a = column_to_rownames(a,"Composite Element REF")
colnames(a)= str_sub(colnames(a),1,15)
# 列名(样本)筛选，第一个条件是有对应的配对pd信息
k1 = str_sub(colnames(a),1,12) %in% pd$patient
table(k1)
#> k1
#> FALSE TRUE
#> 500 80
a = a[,k1]
# 列名筛选，第二个条件是有配对的甲基化表达量
patient = str_sub(colnames(a),1,12)
group_list = ifelse(as.numeric(str_sub(colnames(a),14,15))<10,"Tumor","Normal")
table(group_list)
#> group_list
#> Normal Tumor
#> 32 48
tp = patient[group_list =="Normal"]
np = patient[group_list =="Tumor"]
okp = intersect(tp,np);length(okp)
#> [1] 32
a = a[,patient %in% okp];dim(a)
#> [1] 485577 64
# 对应的修改pd
pd = pd[match(str_sub(colnames(a),1,15),pd$sampleID),]
identical(str_sub(colnames(a),1,15),pd$sampleID)
#> [1] TRUE
2.整理为ChAMP的对象
beta=as.matrix(a)
# beta信号值矩阵里面不能有NA值
beta=impute.knn(beta)
sum(is.na(beta))
beta=beta$data
beta=beta+0.00001
myLoad=champ.filter(beta = beta ,pd = pd) #这一步已经自动完成了过滤
dim(myLoad$beta)
save(myLoad,file = './Rdata/step1_myLoad.Rdata')
网友评论