0.介紹
本文的數(shù)據(jù)出自2017年的一篇文章：Seven-CpG-based prognostic signature coupled with gene expression predicts survival of oral squamous cell carcinoma。參考鏈接?：?TCGA數(shù)據(jù)庫(kù)的癌癥甲基化芯片數(shù)據(jù)重分析
這也是B站的生信技能樹(shù)甲基化芯片課程示例數(shù)據(jù)。經(jīng)過(guò)我的研究和整理，成了今天的推文,這也將作為我以后開(kāi)新課的重要基礎(chǔ)，看到就是賺到~
前面寫(xiě)的甲基化學(xué)習(xí)記錄，有幾張思維導(dǎo)圖：
甲基化學(xué)習(xí)記錄
因?yàn)楫?dāng)時(shí)電腦配置跑不了甲基化數(shù)據(jù)，代碼就沒(méi)有學(xué)成，中途放棄。現(xiàn)在找補(bǔ)回來(lái)！
1.輸入數(shù)據(jù)和R包
從UCSC的xena網(wǎng)站上可以找到TCGA-HNSC的甲基化數(shù)據(jù)和病人的臨床信息。將其存放于raw_data.（這里使用的是GDC-TCGA下載的信號(hào)值矩陣和從TCGA下載的臨床信息和生存信息），信號(hào)值矩陣很大，下下來(lái)才發(fā)現(xiàn)GDC里的臨床信息是沒(méi)有完整的site信息的，遂，轉(zhuǎn)戰(zhàn)TCGA。（實(shí)在不想再次下載信號(hào)值矩陣，將就看）
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里沒(méi)有專(zhuān)門(mén)的OSCC，而是合并為HNSC了。所以需要整個(gè)下載下來(lái)，然后篩選OSCC對(duì)應(yīng)的樣本。
這里就需要用到一個(gè)很重要的符號(hào)：%in%。點(diǎn)擊充電%in%很簡(jiǎn)單
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
可見(jiàn)，604個(gè)樣本中有369個(gè)是屬于OSCC的。已經(jīng)挑選出了這些樣本對(duì)應(yīng)的臨床信息。進(jìn)行簡(jiǎn)化和整理：
1.樣本ID保留至15位，為了和甲基化信號(hào)值矩陣的樣本ID保持一致 2.完整的pd有130多列，從中挑出有用的列 3.生成group_list(Tumor-normal分組)和patient(病人ID)列 4.選出tumor-normal配對(duì)的樣本
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對(duì)病人的臨床信息。接下來(lái)要找他們對(duì)應(yīng)的甲基化信號(hào)值數(shù)據(jù)。
1.2 信號(hào)值矩陣
讀取進(jìn)來(lái)，并篩選有配對(duì)樣本甲基化數(shù)據(jù)的病人。
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)
# 列名(樣本)篩選，第一個(gè)條件是有對(duì)應(yīng)的配對(duì)pd信息
k1 = str_sub(colnames(a),1,12) %in% pd$patient
table(k1)
#> k1
#> FALSE TRUE
#> 500 80
a = a[,k1]
# 列名篩選，第二個(gè)條件是有配對(duì)的甲基化表達(dá)量
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
# 對(duì)應(yīng)的修改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的對(duì)象
beta=as.matrix(a)
# beta信號(hào)值矩陣?yán)锩娌荒苡蠳A值
beta=impute.knn(beta)
sum(is.na(beta))
beta=beta$data
beta=beta+0.00001
myLoad=champ.filter(beta = beta ,pd = pd) #這一步已經(jīng)自動(dòng)完成了過(guò)濾
dim(myLoad$beta)
save(myLoad,file = './Rdata/step1_myLoad.Rdata')
