TCGA甲基化數(shù)據(jù)下載以及相關(guān)臨床信息整理

之前的文章里已經(jīng)按照教程進(jìn)行了TCGA數(shù)據(jù)庫的一些練習(xí)（RNA-seq、芯片、生存分析）拴鸵，現(xiàn)在學(xué)習(xí)TCGA甲基化數(shù)據(jù)的分析過程肄程。

參考文章：
1.甲基化的一些基礎(chǔ)知識
2.Make Decision: DNA甲基化檢測方法，哪一款適合你？
3.一文了解 MethylationEPIC 850K 甲基化芯片
4.DNA甲基化研究方法速遞

摘自參考文章1：

目前5個甲基化數(shù)據(jù)庫：
（1）Roadmap：NIH的“表觀組學(xué)線圖計(jì)劃” (Roadmap Epigenomics Mapping Consortium)，包含367個人類主要組織和細(xì)胞類型的DNA甲基化圖譜
（2）BLUEPRINT：歐洲“血液表觀基因組項(xiàng)目”（BLUEPRINT of Haematopoietic Epigenomes），包含與人類復(fù)雜疾病相關(guān)的82個不同血液細(xì)胞的DNA甲基化圖譜达布。
（3）ENCODE：“DNA元件百科全書”計(jì)劃（The Encyclopedia of DNA Elements），包含來自世界各國32個研究機(jī)構(gòu)對206個人類不同的細(xì)胞系和組織進(jìn)行了DNA甲基化水平的測定逾冬。
（4）ICGC：“國際癌癥基因組聯(lián)盟”（The International Cancer Genome Consortium）黍聂，旨在從基因組、表觀基因組和轉(zhuǎn)錄組等多維數(shù)據(jù)層面研究癌癥的發(fā)生和發(fā)展身腻，包含27種常見癌癥的9000多個樣本的DNA甲基化數(shù)據(jù)产还，
（5）TCGA：包含34種癌癥類型的10000多個樣本的DNA甲基化數(shù)據(jù)，并且保留了癌癥患者詳細(xì)的臨床數(shù)據(jù)資料嘀趟，為生存分析提供了大量的數(shù)據(jù)資源脐区。

練習(xí)的文獻(xiàn)來自：Seven-CpG-based prognostic signature coupled with gene expression predicts survival of oral squamous cell carcinoma。
目標(biāo)：下載TCGA數(shù)據(jù)庫中口腔癌的甲基化芯片信號值矩陣她按，然后挑選有N-T配對的32個病人的數(shù)據(jù)進(jìn)行差異分析牛隅。

這些甲基化結(jié)果使用的平臺是illumina公司的Human methylation 450炕柔，可檢測人全基因組45萬個甲基化位點(diǎn)，全面覆蓋了96%的CpG島媒佣，并根據(jù)需求加入了CpG島以外的CpG位點(diǎn)匕累、人類干細(xì)胞非CpG甲基化位點(diǎn)、正常組織與腫瘤（多種癌癥）組織差異甲基化位點(diǎn)默伍、編碼區(qū)以外的CpG島哩罪、miRNA啟動子區(qū)域和已通過GWAS的疾病相關(guān)區(qū)域的位點(diǎn)，同時覆蓋了Human Methylation27 BeadChip的90%的位點(diǎn)（Human Methylation27是世界首款DNA甲基化芯片）巡验。

練習(xí)之前需要提示的是：

這個TCGA的甲基化數(shù)據(jù)很大，如果你下載完整的文件一共大概有80多個G（TCGAbiolink包下載）碘耳，你需要一個大的儲存空間显设，或者整一個移動硬盤。下載也是很費(fèi)時辛辨，網(wǎng)速不好的話捕捂，過夜下載很正常。而且內(nèi)存要求非常高斗搞，我嘗試了提高系統(tǒng)虛擬內(nèi)存指攒，然而還是不行。于是我就換了一種方法下載僻焚。(見下面）
另外：網(wǎng)上的幾篇教程里面的甲基化數(shù)據(jù)臨床樣品信息似乎和我下載的不太一樣允悦，于是我就自己摸索著整理了一下，結(jié)果得到29對樣品虑啤，與文獻(xiàn)里和教程里的不一樣（32對）隙弛，不過沒有太大關(guān)系，練習(xí)的目的旨在熟悉處理過程狞山。如果有嚴(yán)重強(qiáng)迫癥的童鞋可以移步其他教程全闷，本篇文章絕大多數(shù)代碼是自己寫的，有些繁瑣萍启，但至少是自己研究了2天半的結(jié)果总珠。

（一）數(shù)據(jù)下載

（1）甲基化信號矩陣下載：

去網(wǎng)站：https://xenabrowser.net/datapages/，選擇：

出現(xiàn)新頁面后勘纯，點(diǎn)擊：

再點(diǎn)擊這里：

下載后解壓局服。得到一個名為HumanMethylation450的文件。

（2）下載臨床數(shù)據(jù)

去TCGA網(wǎng)站上選擇好TCGA-HNSC的甲基化數(shù)據(jù)屡律，彈出的頁面應(yīng)該是：

https://portal.gdc.cancer.gov/projects/TCGA-HNSC

把上面的“Biospecimen”和“clinical”兩個部分的“tsv”文件都下載下來：

兩個文件解壓縮后腌逢，會解壓出多個tsv文件，其中一個叫"sample.tsv”和“clinical.tsv”兩個文件的就是我們需要的臨床文件超埋。

（二）數(shù)據(jù)預(yù)處理

（1）臨床數(shù)據(jù)預(yù)處理

首先處理臨床數(shù)據(jù)搏讶，因?yàn)槲覀儾⒉恍枰械臄?shù)據(jù)佳鳖，我們先把樣品進(jìn)行過濾，讀取的是“sample.tsv”文件媒惕，里面有我們要的“腫瘤/正诚捣裕”信息：

#先讀取sample.tsv文件
> pd.all <- read.delim("sample.tsv", header = T, stringsAsFactors = F)
> dim(pd.all) #你會發(fā)現(xiàn)雖然下載頁面介紹只有528個病人，但是這里的樣品信息有1573個妒蔚，說明每個病人的取樣不止1次穿挨，后面要刪掉
[1] 1573   37
> colnames(pd.all) #所有的臨床信息，我們并不需要所有的
#提取部分臨床信息
> pd <- pd.all[,c("case_id","case_submitter_id","sample_submitter_id","sample_type")]
> table(pd$sample_type)  #看一下樣品的類型肴盏，這里有4種科盛，只取后面的兩種，即只取實(shí)體組織
Blood Derived Normal           Metastatic        Primary Tumor  Solid Tissue Normal 
                 511                    2                  978                   82 
#這里看到原位腫瘤有978個樣品菜皂，遠(yuǎn)遠(yuǎn)超過了528個病人的數(shù)量贞绵，再一次驗(yàn)證了存在同一個病人多次取樣的情況
> tissue = c("Primary Tumor","Solid Tissue Normal") #取后面兩種樣品類型
> pd_tissue <- pd[pd$sample_type %in% tissue,] #從所有的臨床信息里提取我們要的兩種類型的信息
> table(pd_tissue$sample_type)
        Primary Tumor       Solid Tissue Normal 
                978                  82

截至目前為止，提取出了實(shí)體組織的樣品恍飘。下面要對樣品進(jìn)行過濾榨崩，因?yàn)橛行┠[瘤樣品的取樣多于1次：

> table(pd_tissue$case_submitter_id) #查看case_submitter_id（例如：TCGA-CQ-A4CE）分別出現(xiàn)了幾次
#結(jié)果里有1,2,3。
#首先我們要先去掉1的那些樣品章母，因?yàn)檫@些樣品一定不是配對的母蛛。
#刪除在case_submitter_id（例如：TCGA-CQ-A4CE）里只出現(xiàn)過一次的行
#代碼參考：https://blog.csdn.net/dingchenxixi/article/details/50865277
> deleteuniquelines <- function(x) {# x為輸入的數(shù)據(jù)框
  stand.col <- x$case_submitter_id 
  count <- table(stand.col) #table函數(shù)可以得到每個上述列每個數(shù)所出現(xiàn)的頻數(shù)
  if (all(count < 2)) stop("no repeated records") 
  else {
    ind <- sapply(stand.col, function(t) ifelse(count[as.character(t)] > 1, TRUE, FALSE))
  }
  return(x[ind, ])
}
> pd_tissue_filtered = deleteuniquelines(pd_tissue)
> dim(pd_tissue_filtered) #過濾完不配對的樣品，還剩982個乳怎，但這里面不是全配對的彩郊，有的是同一個腫瘤樣品取了兩次
#[1] 982   4
#分別提取正常和腫瘤組織樣品
> nt = pd_tissue_filtered[pd_tissue_filtered$sample_type =="Solid Tissue Normal",]
> tt = pd_tissue_filtered[pd_tissue_filtered$sample_type =="Primary Tumor",]
#對于正常組織，由于都只取了一次樣品蚪缀，所以不進(jìn)行過濾
#對于腫瘤樣品：只取tt里sample_submitter_id編碼最后一位是"A"的樣品焦辅，因?yàn)锽是福爾馬林固定石蠟包埋組織
#從B以后就不太好了，如果你table一下會發(fā)現(xiàn)還有Z椿胯，所以只取A的樣品
> tt <- tt[substr(tt$sample_submitter_id,16,16) =="A",]
#取和正常對照匹配的腫瘤樣品
> tt <- tt[tt$case_submitter_id %in% nt$case_submitter_id,]
> dim(tt)
[1] 82  4
#合并正常和腫瘤樣品
> paired_tissue <- rbind(nt,tt)
> dim(paired_tissue)
[1] 164   4

到目前為止筷登，我們提取出了成對的樣品（在臨床信息里）。一共是82對哩盲。下面我們需要對腫瘤發(fā)生的位置進(jìn)行篩選前方，只取發(fā)生在口腔里的腫瘤組織和其對照，需要讀取“clinical.tsv”文件：

#讀取“clinical.tsv”文件廉油，提取腫瘤位置信息
> tumor_site <-read.delim("clinical.tsv",header = T,stringsAsFactors = F)
#這里需要注意的是惠险，臨床樣品里有兩列分別是“取樣位置”和“腫瘤發(fā)生位置”，你要取的是“腫瘤發(fā)生的位置”
> tumor_site <- tumor_site[,c("case_id","case_submitter_id","tissue_or_organ_of_origin")]
> library(stringr)
> tumor_site[, c("anatomic_neoplasm_subdivision", "takeout")] <- str_split_fixed(tumor_site$tissue_or_organ_of_origin, ",", 2)
> tumor_site <- tumor_site[,c(1,2,4)]
> table(tumor_site$anatomic_neoplasm_subdivision) #看一下所有的腫瘤位置信息
#這里的分類信息和有些教程里的不一樣抒线，我覺得無所謂
     Anterior floor of mouth               Base of tongue             Border of tongue 
                           4                           48                            2 
                Cheek mucosa               Floor of mouth                          Gum 
                          38                          108                           16 
                 Hard palate                  Hypopharynx                       Larynx 
                           8                           18                          232 
                         Lip                    Lower gum                     Mandible 
                           6                            4                            2 
                       Mouth                   Oropharynx    Overlapping lesion of lip 
                          46                           18                          140 
                      Palate                      Pharynx Posterior wall of oropharynx 
                           2                            2                            2 
             Retromolar area                 Supraglottis                       Tongue 
                           2                            2                          260 
                      Tonsil                    Upper Gum    Ventral surface of tongue 
                          92                            2                            2 
> oscc = c("Mouth",
           "Cheek mucosa",
           "Lip",
           "Hard palate",
           "Floor of mouth",
           "gum","Upper Gum","Lower gum",
           "Anterior floor of mouth",
           "Border of tongue","Base of tongue","Tongue",
           "Ventral surface of tongue") #取口腔里發(fā)生的腫瘤
> tumor_site <- tumor_site[tumor_site$anatomic_neoplasm_subdivision %in% oscc,]
> dim(tumor_site) #我們?nèi)〉哪[瘤發(fā)生在口腔里的樣品有528個
[1]528   3

我們需要知道上面82對樣品里有多少是發(fā)生在口腔里的腫瘤樣品：

> tumor_site_unique <- unique(tumor_site$case_submitter_id) #取“腫瘤位置”矩陣?yán)锍霈F(xiàn)過的樣品
> tumor_site_unique = as.data.frame(tumor_site_unique)
> colnames(tumor_site_unique) = "case_submitter_id"
#把“腫瘤位置”里的樣品匹配到“配對腫瘤/正嘲喙”樣品里
> merge_info <- paired_tissue[paired_tissue$case_submitter_id %in% tumor_site_unique$case_submitter_id,]
> dim(merge_info) #這里有36對樣品
[1] 72   4

注意！Ｋ惶俊抱慌！這里不要以為就處理完臨床樣品了逊桦，因?yàn)槟阋粫阂雅R床樣品和甲基化信號矩陣做交集的，所以你還得看你的id和甲基化信號矩陣的列名是不是一樣的格式Ｒ纸Ｇ烤！

讀取甲基化信號矩陣：

> methy_data <- data.table::fread("HumanMethylation450",data.table = F)
|--------------------------------------------------|
|==================================================|
> methy_data[1:4,1:4]
      sample TCGA-CN-6022-01 TCGA-CQ-5327-01 TCGA-CV-5435-11
1 cg13332474          0.0220          0.0264          0.0307
2 cg00651829          0.0209          0.4514          0.0214
3 cg17027195          0.0443          0.2414          0.0330
4 cg09868354          0.0490          0.0546          0.0469

你會發(fā)現(xiàn)寺渗，這個列名是sample_submitter_id的前15位字符匿情，我們上面的sample_submitter_id共有16位字符，所以還要處理一下我們的臨床信息：

> merge_info$sample_submitter_id = substr(merge_info$sample_submitter_id,1,15)
> head(merge_info)
      case_id case_submitter_id sample_submitter_id
3   e7d1f0dd-eec0-4670-a2ba-00ffe38c6382      TCGA-CV-6939     TCGA-CV-6939-11
8   ebfb8c1b-471f-4acb-a1fe-287e02cadba7      TCGA-H7-A6C4     TCGA-H7-A6C4-11
43  91728cd6-ed91-49dc-9279-faedbe211a9d      TCGA-CV-6951     TCGA-CV-6951-11
73  dec3c21b-3f59-472f-8573-27b0e830aa92      TCGA-HD-A6I0     TCGA-HD-A6I0-11
91  93a265b7-6c23-4e5f-b797-c117793744bf      TCGA-CV-6936     TCGA-CV-6936-11
113 f9e82d62-fafc-45b0-97cf-d1a4a3c884a1      TCGA-WA-A7GZ     TCGA-WA-A7GZ-11
            sample_type
3   Solid Tissue Normal
8   Solid Tissue Normal
43  Solid Tissue Normal
73  Solid Tissue Normal
91  Solid Tissue Normal
113 Solid Tissue Normal
#最后千萬別忘記保存
> write.table(merge_info$sample_submitter_id, file = "paired_samples_name_from_clinical.txt",quote = F, row.names = F, col.names = F)

（2）甲基化矩陣預(yù)處理

從580個甲基化樣品里提取成對樣品的甲基化信號
在R里操作：

#讀取甲基化矩陣信殊，這個矩陣很大炬称，需要fread讀取
> methy_data <- data.table::fread("HumanMethylation450",data.table = F)
> methy_data[1:4,1:4]
#將上面36對樣品的腫瘤和正常分別提取出來
> nt_paired = merge_info[merge_info$sample_type=="Solid Tissue Normal",]
> tt_paired = merge_info[merge_info$sample_type=="Primary Tumor",]
#與甲基化信號矩陣進(jìn)行匹配
> methy_nt = methy_data[,colnames(methy_data) %in% nt_paired$sample_submitter_id] #甲基化矩陣?yán)镏挥?9個正常樣品和臨床信息配對
> methy_tt = methy_data[,colnames(methy_data) %in% tt_paired$sample_submitter_id] #甲基化矩陣?yán)镉?6個腫瘤樣品，和臨床信息相符
#匹配后合并腫瘤與正常樣品
> methy_combine = cbind(methy_nt,methy_tt)
> table(colnames(methy_combine)) #這里并不是每一個腫瘤樣品都匹配一個正常對照涡拘，所以我們要把多余的腫瘤樣品去掉
#你會發(fā)現(xiàn)最后7個樣品并不是配對的转砖，我們需要把它們刪掉
TCGA-CV-5436-01 TCGA-CV-5436-11 TCGA-CV-5439-01 TCGA-CV-5439-11 TCGA-CV-5442-01 TCGA-CV-5442-11 
              1               1               1               1               1               1 
TCGA-CV-5970-01 TCGA-CV-5970-11 TCGA-CV-5971-01 TCGA-CV-5971-11 TCGA-CV-5973-01 TCGA-CV-5973-11 
              1               1               1               1               1               1 
TCGA-CV-5976-01 TCGA-CV-5976-11 TCGA-CV-5977-01 TCGA-CV-5977-11 TCGA-CV-5979-01 TCGA-CV-5979-11 
              1               1               1               1               1               1 
TCGA-CV-6003-01 TCGA-CV-6003-11 TCGA-CV-6433-01 TCGA-CV-6433-11 TCGA-CV-6436-01 TCGA-CV-6436-11 
              1               1               1               1               1               1 
TCGA-CV-6441-01 TCGA-CV-6441-11 TCGA-CV-6933-01 TCGA-CV-6933-11 TCGA-CV-6934-01 TCGA-CV-6934-11 
              1               1               1               1               1               1 
TCGA-CV-6936-01 TCGA-CV-6936-11 TCGA-CV-6939-01 TCGA-CV-6939-11 TCGA-CV-6943-01 TCGA-CV-6943-11 
              1               1               1               1               1               1 
TCGA-CV-6951-01 TCGA-CV-6951-11 TCGA-CV-6952-01 TCGA-CV-6952-11 TCGA-CV-6953-01 TCGA-CV-6953-11 
              1               1               1               1               1               1 
TCGA-CV-6954-01 TCGA-CV-6954-11 TCGA-CV-6956-01 TCGA-CV-6956-11 TCGA-CV-6959-01 TCGA-CV-6959-11 
              1               1               1               1               1               1 
TCGA-CV-6961-01 TCGA-CV-6961-11 TCGA-CV-7103-01 TCGA-CV-7103-11 TCGA-CV-7235-01 TCGA-CV-7235-11 
              1               1               1               1               1               1 
TCGA-CV-7238-01 TCGA-CV-7238-11 TCGA-CV-7255-01 TCGA-CV-7255-11 TCGA-CV-7406-01 TCGA-CV-7438-01 
              1               1               1               1               1               1 
TCGA-H7-A6C4-01 TCGA-HD-8635-01 TCGA-HD-A6HZ-01 TCGA-HD-A6I0-01 TCGA-WA-A7GZ-01 
              1               1               1               1               1 

#把只出現(xiàn)過一次列名的列去掉
> deleteuniquecolumn <- function(x) {# x為輸入的數(shù)據(jù)框
  m = substr(colnames(x),1,12)
  stand.col <- m 
  count <- table(stand.col) #table函數(shù)可以得到每個上述列每個數(shù)所出現(xiàn)的頻數(shù)
  if (all(count < 2)) stop("no repeated records") 
  else {
    ind <- sapply(stand.col, function(t) ifelse(count[as.character(t)] > 1, TRUE, FALSE))
  }
  return(x[,ind])
}
> methy_paired = deleteuniquecolumn(methy_combine)
> table(colnames(methy_paired))
> rownames(methy_paired) = methy_data$sample
#保存過濾完樣品的甲基化矩陣
> write.csv(methy_paired,file="OSCC_29paired_methydata.csv")

到這里，我們把甲基化矩陣過濾完成了鲸伴，留下29對口腔腫瘤/正常樣品對，這個結(jié)果與文獻(xiàn)里的不同晋控，我認(rèn)為無所謂汞窗，主要是走一下流程。有可能是我和文獻(xiàn)里過濾樣品的方法不同也有關(guān)系赡译。

過濾樣品后仲吏，甲基化矩陣長這樣：

（3）過濾甲基化矩陣，并整理為ChAMP對象

#讀取上面保存的29對樣品甲基化信號矩陣
> methy = data.table::fread("OSCC_29paired_methydata.csv",data.table = F)
> head(methy)
#載入需要的包
> library(ChAMP)
> library(dplyr)
> library(tibble)
#把myNorm里的列名按照pd里的sample_submitter_id順序排
> methy_sort <- methy[,c("V1","TCGA-CV-5436-01","TCGA-CV-5436-11","TCGA-CV-5439-01",
                       "TCGA-CV-5439-11","TCGA-CV-5442-01","TCGA-CV-5442-11","TCGA-CV-5970-01",
                      "TCGA-CV-5970-11","TCGA-CV-5971-01","TCGA-CV-5971-11","TCGA-CV-5973-01",
                       "TCGA-CV-5973-11","TCGA-CV-5976-01","TCGA-CV-5976-11","TCGA-CV-5977-01",
                       "TCGA-CV-5977-11","TCGA-CV-5979-01","TCGA-CV-5979-11","TCGA-CV-6003-01",
                       "TCGA-CV-6003-11","TCGA-CV-6433-01","TCGA-CV-6433-11","TCGA-CV-6436-01",
                       "TCGA-CV-6436-11","TCGA-CV-6441-01","TCGA-CV-6441-11","TCGA-CV-6933-01",
                       "TCGA-CV-6933-11","TCGA-CV-6934-01","TCGA-CV-6934-11","TCGA-CV-6936-01",
                       "TCGA-CV-6936-11","TCGA-CV-6939-01","TCGA-CV-6939-11","TCGA-CV-6943-01",
                       "TCGA-CV-6943-11","TCGA-CV-6951-01","TCGA-CV-6951-11","TCGA-CV-6952-01",
                       "TCGA-CV-6952-11","TCGA-CV-6953-01","TCGA-CV-6953-11","TCGA-CV-6954-01",
                       "TCGA-CV-6954-11","TCGA-CV-6956-01","TCGA-CV-6956-11","TCGA-CV-6959-01",
                       "TCGA-CV-6959-11","TCGA-CV-6961-01","TCGA-CV-6961-11","TCGA-CV-7103-01",
                       "TCGA-CV-7103-11","TCGA-CV-7235-01","TCGA-CV-7235-11","TCGA-CV-7238-01",
                       "TCGA-CV-7238-11","TCGA-CV-7255-01","TCGA-CV-7255-11")]

> a = column_to_rownames(methy_sort,"V1")
> beta_value = as.matrix(a)
# beta信號值矩陣?yán)锩娌荒苡蠳A值
> beta=impute.knn(beta_value) 
> sum(is.na(beta)) #這里是檢查矩陣?yán)锸欠襁€有NA
#[1] 0
> beta=beta$data
> beta=beta+0.00001
#準(zhǔn)備pd表型文件（實(shí)際上就是樣品的信息）
> pd_1 <- as.data.frame(colnames(beta))
> pd_info <- merge_info[merge_info$sample_submitter_id %in% pd_1$`colnames(beta)`,]
> colnames(pd_1) = "sample_submitter_id"
> pd <- merge(pd_1,pd_info,by = "sample_submitter_id",all.x = TRUE)
#使用ChAMP過濾
> myLoad=champ.filter(beta = beta ,pd = pd) #這一步已經(jīng)自動完成了過濾
> dim(myLoad$beta) #beta就是指的是beta值蝌焚，我們需要的甲基化信號
#[1] 412481     58
#保存這個ChAMP對象
> save(myLoad,file = 'OSCC_29paired_methydata_ChAMPfiltered.Rdata')

了解一下上面的champ.filter這一步都過濾了些什么：

1.過濾掉detection p-value大于0.01的探針
2.過濾掉在至少5%樣本中bead count小于3的探針
3.過濾掉非GpC位點(diǎn)的探針
4.過濾掉所有SNP相關(guān)的探針
5.過濾掉multi-hit探針裹唆，即映射到多個位置的
6.過濾掉X和Y染色體上的探針
參考文章：甲基化芯片入門學(xué)習(xí)-ChAMP包（二）

看一下這個對象長啥樣：