TCGA 甲基化芯片数据下载和整理

测试开头

今天是生信星球陪你的第611天

大神一句话，菜鸟跑半年。我不是大神，但我可以缩短你走弯路的半年~

就像歌儿唱的那样，如果你不知道该往哪儿走，就留在这学点生信好不好~

这里有豆豆和花花的学习历程，从新手到进阶，生信路上有你有我！

花花写于2020-4-30，是努力薅代码的一天！

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。（实在不想再次下载信号值矩阵，将就看）

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