使用DSS包多種方式檢驗差異甲基化訊號區域

  • 2020 年 3 月 27 日
  • 筆記

因為一個有趣的原因，我們來一起學一個R包吧

一個背景

哺乳動物基因組CpG位點通常集中在稱為CpG島（CpG island，CGI）的區域中，並且已知人基因啟動子~60％含有CpG島。CpG島上下游不超過2000個鹼基對（2kb）的基因組區域稱為CpG「島岸」（shores），其中CpG shelves指位於CpG shores 上下游2kb以內的區域，open sea指CpG islands、CpG shores和CpG shelves之外的其他區域。這4種情況形成了CpG resort。CpG位點的密度從island到open sea遞減。

3個技術

主要是 甲基化測序的 WGBS和RRBS，還有 晶片：

全基因組DNA甲基化測序（Whole Genome Bisulfite Sequencing，WGBS）是 DNA 甲基化研究的金標準，它通過 Bisulfite 處理和全基因組 DNA 測序結合的方式，對整個基因組上的甲基化情況進行分析，具有單鹼基解析度，可精確評估單個 C 鹼基的甲基化水平，構建全基因組精細甲基化圖譜。數據量非常大。

簡化甲基化測序 (Reduced representation bisulfite sequencing, RRBS)是一種準確、高效、經濟的DNA甲基化研究方法，通過酶切 (Msp I) 富集啟動子及CpG島區域，並進行Bisulfite測序，同時實現DNA甲基化狀態檢測的高解析度和測序數據的高利用率。作為一種高性價比的甲基化研究方法，簡化甲基化測序在大規模臨床樣本的研究中具有廣泛的應用前景。

Illumina的Infinium BeadChip晶片，包括HumanMethyation450（450K）和MethylationEPIC（850K）。Infinium晶片存在染料偏差、不同探針化學和位置效應的問題，已知這些問題會影響結果，必須在數據處理過程中進行校正。Infinium 450K探針交叉反應和模糊比對到人類基因組中的多個位置影響了485,000個探測器中的約140,000個探針（29％），將可用探針的數量減少到約345,000個。這個問題在新發布850K仍然存在，其包括> 90％的450K探針。

有文章比較這3個技術：Empirical comparison of reduced representation bisulfite sequencing and Infinium BeadChip reproducibility and coverage of DNA methylation in humans

主要是分析 差異甲基化區域（DMRs）與 DMR 相關差異表達基因

數據介紹

這裡選擇 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE52140

提供每個樣本的訊號值矩陣下載

下載並且了解數據:

查看壓縮包內容，如下：

head GSM1084240_A3R_d6_250.cpgs.txt  CHR POS A3R_d6_250  chr1 10525 21/23  chr1 10542 21/23  chr1 10609 1/23  chr1 10617 0/24  chr1 10620 0/24  chr1 10631 0/24  chr1 10633 0/24  chr1 10636 0/24  chr1 10638 0/24

DSS包要求輸入文件數據:每一行代表一個CpG site, 格式如下:

  • 第一列為染色體
  • 第二列為位置
  • 第三列為total reads
  • 第四列為甲基化的reads

所以我們下載的數據需要進行拆分，然後導入到R裡面才能被DSS包使用。

DSS包介紹

主要是把上面項目的數據文件下載，然後導入到R裡面，是有DSS包進行分析。

DSS (Dispersion Shrinkage for Sequencing data)，為基於高通量測序數據的差異分析而設計的Bioconductor包。主要應用於BS-seq（亞硫酸氫鹽測序）中計算不同組別間差異甲基化位點（DML）和差異甲基化區域（DMR）即Call DML or DMR。

DSS包的使用主要包括：

  • 輸入文件的準備
  • 利用DMLtest函數檢驗所有的位點
  • 利用callDML函數挑選統計學顯著的位點
  • 利用callDMR函數Call DMR
  • 利用showOneDMR函數對DMRs可視化

首先我們導入上面GSE52140數據集的文件：

library(data.table)  library(stringr)  library(tidyverse)  allDat <- lapply(list.files('GSE52140_RAW/',pattern='*cpgs.txt.gz'),function(f){    # f="GSM1251242_H2R_d0.cpgs.txt.gz";    print(f);    tmp=fread(file.path('GSE52140_RAW/',f))    chr=as.character(tmp$CHR)    pos=as.character(tmp$POS)    newTmp=separate(tmp,col =3,into = c("methy", "unmethy"), sep = "/")    newTmp$all=as.numeric(newTmp$methy)+as.numeric(newTmp$unmethy)    newTmp=as.data.frame(newTmp[,c(1,2,5,3)])    colnames(newTmp)=c('chr', 'pos' ,'N' ,'X')    return(newTmp)  })    ## 值得注意的是每個樣本的位點數量不一致哦  do.call(rbind,lapply(allDat,dim))  tmp=do.call(cbind,lapply(allDat,head))    sn=gsub('.cpgs.txt.gz','',list.files('GSE52140_RAW/',pattern='*cpgs.txt.gz'))  sn=gsub('GSM.*?_','',sn)  sn

也就是說把這17個文件讀入了，樣本名字是：

> sn   [1] "A0R_d0_rep1"  "A3R_d0_rep1"  "A3R_d6_250"   "A3R_d6_1000"   [5] "A3R_d13_250"  "A3R_d13_1000" "H0R_d0_rep1"  "H3R_d0_rep1"   [9] "H3R_d6_250"   "H3R_d13_250"  "A0R_d0_rep2"  "A3R_d0_rep2"  [13] "H0R_d0_rep2"  "H3R_d0_rep2"  "A1R_d0"       "A2R_d0"  [17] "H2R_d0"

這個時候，這個變數有點大，可能會考驗你的電腦哦。

然後我們用下面的3個例子來說明這個DSS包的用法，需要掌握上面樣本的命名：

# lung cancer cell lines A549 (A) and HTB56 (H)  # normal cell lines (0R)  # a highly metastatic phenotype (3R)  #  5-Azacytidine treatment at low concentrations (250 nM & 1000 nM)  # for 6 days, additional 7 days in regular medium

單樣本VS單樣本

程式碼如下，重要就是構建對象和做統計檢驗

library(DSS)  require(bsseq)  if(T){    BSobj <- makeBSseqData(allDat[1:2],                           c("A0R", "A3R") )[1:1000,]    BSobj     save(BSobj,file = 'single-BSobj.Rdata')    # There is no biological replicates in at least one condition.    dmlTest <- DMLtest(BSobj, group1=c("A0R"), group2=c("A3R"),smoothing=TRUE)    head(dmlTest)  }

多樣本的組VS另一個組

程式碼如下，重要就是構建對象和做統計檢驗，這裡比較"A0R_d0"和"A3R_d0"組別的2個樣本：

if(T){    sn[c(1,11,2,12)]    BSobj <- makeBSseqData(allDat[c(1,11,2,12)],                           sn[c(1,11,2,12)] )[1:1000,]     BSobj    save(BSobj,file = 'group-BSobj.Rdata')    dmlTest <- DMLtest(BSobj, group1=c("A0R_d0_rep1","A0R_d0_rep2"),                       group2=c("A3R_d0_rep1","A3R_d0_rep2"),smoothing=F)    head(dmlTest)  }

多種比較方式

程式碼如下，重要仍然是構建對象和做統計檢驗，但是需要構建屬性矩陣，而且增加了DMLfit步驟。

  sn    sn[grep('rep',sn)]    cellline=substring(sn[grep('rep',sn)],1,1)    type=substring(sn[grep('rep',sn)],2,2)    design=data.frame(cellline=cellline,type=type)    design

得到的屬性矩陣如下：

  sn    sn[grep('rep',sn)]    cellline=substring(sn[grep('rep',sn)],1,1)    type=substring(sn[grep('rep',sn)],2,2)    design=data.frame(cellline=cellline,type=type)    design      # 構建對象特別耗時；    BSobj <- makeBSseqData(allDat[c(grep('rep',sn))],                           sn[grep('rep',sn)])    BSobj    save(BSobj,file = 'multi-BSobj.Rdata')    load(file = 'multi-BSobj.Rdata')    DMLfit=DMLfit.multiFactor(BSobj,design,                              formula = ~cellline+type+cellline:type)      colnames(DMLfit$X)    # 這裡可以使用 『coef』, 『term』, or 『Contrast』我們僅僅是演示 coef    dmlTest=DMLtest.multiFactor(DMLfit,coef=2)    head(dmlTest)

值得注意的是DMLtest.multiFactor結果不需要callDML，只需要callDMR即可!

結果介紹

不管是哪種比較，最後都得到dmlTest變數走後面的流程，包括確定顯著差異甲基化區域及基因，以及可視化展現，程式碼如下：

# 3.Call DML by using callDML function. The results DMLs are sorted by the significance.  dmls <- callDML(dmlTest, p.threshold=0.001)  head(dmls)  ##To detect loci with difference greater than 0.1, do:  dmls2 <- callDML(dmlTest, delta=0.1, p.threshold=0.001)  head(dmls2)    # 4.Call DMR by using callDML function  ##Regions with many statistically significant CpG sites are identified as DMRs.  dmrs <- callDMR(dmlTest, p.threshold=0.01)  head(dmrs)  ##To detect regions with difference greater than 0.1, do:  dmrs2 <- callDMR(dmlTest, delta=0.1, p.threshold=0.05)  head(dmrs2)    # 5.The DMRs can be visualized using showOneDMR function  showOneDMR(dmrs[1,], BSobj)

很明顯，參數都是可以調整的，統計學顯著性的閾值自己把握。

作業

分析文章 Enduring epigenetic landmarks define the cancer microenvironment ，拿到 患者間差異甲基化區域（DMRs）與 DMR 相關差異表達基因（DE-DMRs）

VirMach 便宜 VPS

QNews

QNews