【原】3種方法注釋你的甲基化探針

健明 2021-07-14

展開全文

關于表達芯片的公共數據庫挖掘我這邊以及差不多把改寫的推文在2年前就寫完了，但表達芯片畢竟只占芯片市場的半壁江山,，還有大量的非表達芯片,，比如大名鼎鼎的甲基化芯片。關于甲基化,，我們公眾號教程非常少,，主要是因為我本人在短暫的6年生物信息學工作經驗中并沒有實際負責過這樣的項目，而我們公眾號90%教程都是我寫的,，極少的投稿里面,，只有 850K甲基化芯片數據的分析讓我印象深刻，所以我的學徒也是跟著這個教程在學習甲基化芯片數據處理,。
但是學徒反饋說這個教程沒有提到芯片的探針如何注釋到基因,，所以就安排他做了這個補充，下面讓我們看學徒的表演,！

甲基化芯片背景

甲基化芯片原理：https://www.jianshu.com/p/c4f758e0399d

芯片主要分為EPIC和450k兩種,，EIPC也就是850k，兩種探針的都是以cg開頭的數字編號,，所謂注釋也就是提取這些探針的所對應的信息,，例如,，探針序列的CpG位置信息，對應的基因信息,，染色體上的位置信息,，等等。很多包在安裝的時候都會自動下載這些注釋信息,，并包裝在一起,，如果我們想要自己注釋這些探針，就要考慮如何獲取獨立的注釋信息,。而所需要注釋數據的,，大部分都來自于兩個數據庫，GEO和TCGA,。

下面介紹三種提取注釋信息的方法

# 方法一：從UCSC Xena下載

直接從UCSC Xena相應的癌癥甲基化數據庫里下載對應的文件,。可以看到是來自GPL16304平臺的芯片,，其實和下面要介紹的從GEO下載注釋信息是一樣的,，不過TCGA的探針數可能會少于45w，大約39w,，因為提前過濾了一些低質量的探針,。

# 方法二：從GEO下載對應平臺的注釋文件

在GEO的官網platform下搜索Illumina HumanMethylation450，可以看到450k的芯片主要來自三個平臺,，探針數也是不一樣的,，TCGA中下載時一般都會標明來自那個平臺，從GEO中下載數據都會得知平臺的信息,。直接進入對應平臺的介紹就可以了,。

看一下最常見的GPL13534平臺的內容

可以看到有1295個GSE數據集來自這個平臺，可以利用的數據相當多,，這里給出了一部分數據的概覽

表格中展示了部分信息,，直接下載然后就可以提取我們需要的注釋信息了,，485577個探針一個不差,，可能是因為我網速的問題，只有下載CSV這個的時候速度比較快,，其他速度都非常感人

可以看到信息非常全面了,，但實際上我們并用不到這么多，有下面這些就夠了

ids<-a[,c("IlmnID",
          "Infinium_Design_Type",
          "CHR",
          "UCSC_RefGene_Name",
          "UCSC_RefGene_Group",
          "UCSC_CpG_Islands_Name",
          "Relation_to_UCSC_CpG_Island")]

#可以看到,，包含了探針的名稱,，探針類型（可用來去批次），染色體編號,，對應的基因名,，在染色體上的位置信息,，以及與CpG島的位置信息，這樣就可以注釋差異的探針了

# 方法三：從ChAMP包中提取

這個方法嚴格來說其實是從ChAMP依賴的兩個注釋包中提取的,，但是我又懶又笨,，懶得看原始的包里數據藏在哪里了，ChAMP包在做甲基化分析的時候也很方便,，而其中champ.filter函數直接就提取好了

myimport <- champ.import(directory=system.file("extdata",package="ChAMPdata"))
myImport=myimport#包里的演示代碼有個小細節(jié)錯了,，沒有區(qū)分大小寫，無傷大雅的
myfilter <- champ.filter(beta=myImport$beta,pd=myImport$pd,detP=myImport$detP,beadcount=myImport$beadcount)
View(hm450.manifest.hg19)
View(probe.features)#兩個信息是一樣的,，包已經把我們需要的信息自動提取了