1.1输入:

这个小插图需要大的输入文件，可以在这里下载:MIRA_Ewing_Vignette_Files.tgz．这个档案包含:

• 本文分析了12例尤文氏肉瘤和12例肌管(肌肉)标本的单核苷酸分辨率DNA甲基化数据(RRBS数据)DNA甲基化异质性决定了尤文氏肉瘤的病谱．
  
• 两组区域-一组尤因特异性DHSs，另一组肌肉特异性DHSs监管元素数据库．所有数据都在hg38参考基因组上。

2.1最初的贯通

首先，让我们加载库并读入数据文件。包中包含了一个示例注释文件:

library(MIRA) library(data.table) #用于函数:fread, setkey, merge library(genome icranges) #用于函数:GRanges, resize library(ggplot2) #用于函数:ylab exampleAnnoDT2 <- fread(system. table) #文件(“extdata”、“exampleAnnoDT2.txt”、包=“米拉”))

接下来，为DNA甲基化样本和区域集构造文件名(通过输入部分)。

# 12 Ewing样本:T1-T12 pathToData <- "/path/to/MIRA_Ewing_Vignette_Files/" ewingFiles <- paste0(pathToData， "EWS_T"， 1:12， ".bed") # 12肌肉相关样本，每种类型3个muscleFiles <- c("Hsmm_"， "Hsmmt_"， "Hsmmfshd_"，"Hsmmtubefshd_") muscleFiles <- paste0(pathToData, muscleFiles, rep(1:3, each=4)， ".bed") RRBSFiles <- c(ewingFiles, muscleFiles) regionSetFiles <- paste0(pathToData, muscleFiles, rep(1:3, each=4)， "sknmc_specific. txt ")床”、“muscle_specific.bed”))

让我们把文件读入R:

BSDTList <- lapply(X=RRBSFiles, FUN=BSreadBiSeq) regionSets <- lapply(X=regionSetFiles, FUN=fread)

在将数据插入MIRA之前，我们需要做一些预处理。的BSreadBiSeq函数从软件Biseqmethcalling读取DNA甲基化调用，将甲基化调用读入data.table对象。我们在区域内阅读从文件中读他们也将如此data.table对象。首先，让我们把样本名称添加到甲基化的列表中data.tables．然后我们把区域集从data.table对象农庄对象，包括染色体、开始和结束信息。在大多数情况下，MIRA不需要链信息，应该省略。

names(BSDTList) <- tools::file_path_sans_ext(basename(RRBSFiles)) regionSets <- lapply(X=regionSets, FUN=function(X) setNames(X, c("chr"， "start"， "end")) regionSets <- lapply(X=regionSets, FUN=function(X) GRanges(seqnames= X $chr, ranges=IRanges(X $start, X $end))

我们还需要选择初始区域大小。为了获得我们感兴趣地点周围的甲基化剖面，我们将区域大小调整到5000 bp。此外，对集合中的所有区域使用相同的区域大小将使最终的甲基化概要更容易解释。与区域集的平均区域大小(两个区域集都在151 bp左右)相比，这是一个较宽的初始猜测，但我们稍后将对此进行调整。

regionSets[[1]] <- resize(regionSets[[1]]， 5000, fix="center") regionSets[[2]] <- resize(regionSets[[2]]， 5000, fix="center") names(regionSets) <- c("Ewing_Specific"， "Muscle_Specific")

现在我们准备在样本的一个子集上运行MIRA——三个尤因样本和三个肌肉样本。首先，我们将甲基化数据聚合到aggregateMethyl．MIRA将每个区域划分为不同的容器，并在每个区域的每个容器中发现甲基化。然后在所有区域聚合匹配的箱子(所有的箱子1在一起，所有的箱子2在一起，等等)。垃圾箱编号(binNum)将影响数据的分辨率或噪声——较高的仓数将允许更高的分辨率，因为每个仓对应的基因组区域较小，但由于每个仓中的甲基化数据较少，可能会增加噪声。为了对称和MIRA评分，建议使用奇数仓号。

subBSDTList <- BSDTList[c(1,5,9,13,17,21)] bigBin <- lapply(X=subBSDTList, FUN= aggregatem乙基，GRList=regionSets, binNum=21, minBaseCovPerBin=0)

然后我们将所有样本合并为一个data.table并根据列表项的名称为示例名称添加一列:

bigBinDT1 <- rbindNamedList(bigBin, newColName = "sampleName")

如果您在自己的设备上运行实际示例，则可以跳过以下步骤。否则，让我们加载打包的甲基化数据(bigBinDT1)，这是通过上一步得到的。

MIRA包数据(bigBinDT1)

接下来我们添加注释数据:每个样本的实验条件或样本类型。这是根据样本类型对样本进行分组所必需的，我们将在绘制甲基化剖面和MIRA评分时进行分组。我们将使用函数data.table包来做此操作，合并基于sampleName列。

setkey(bigBinDT1, sampleName) bigBinDT1 <- merge(bigBinDT1, exampleAnnoDT2, all.x=TRUE)

让我们看看MIRA概要文件，看看是否应该更改区域大小。

plotMIRAProfiles (binnedRegDT = bigBinDT1)

2．2选择合适的区域大小

虽然甲基化剖面看起来很好，但当观察肌肉特定区域集时，对肌肉相关样本来说，降幅有点太窄。如果凹陷过窄，MIRA中的评分函数可能会遇到问题，因为它预计凹陷从外缘到外缘至少有5个箱宽(这个数字包括外缘)。让我们对两个区域集使用较小的区域大小:尤因DHS区域为4000 bp，肌肉DHS区域为1750 bp(选择有点随意)。区域的大小不宜缩小太多;如果区域太小，可能会有更多的噪音，因为每个容器的甲基化数据会更少。此外，区域必须足够宽，以捕捉倾角的外边缘，因为外边缘用于MIRA评分。

要实现我们的新区域大小:

regionSets[[1]] <- resize(regionSets[[1]]， 4000, fix="center") # Ewing regionSets[[2]] <- resize(regionSets[[2]]， 1750, fix="center") #肌肉

2．3全面的米拉分析

现在我们可以用所有样本运行MIRA。让我们把甲基化数据再次汇总aggregateMethyl．

(X=BSDTList, FUN= aggregatem乙基，GRList=regionSets, binNum=21, minBaseCovPerBin=0) bigBinDT2 <- rbindNamedList(bigBin)

如果您在自己的设备上运行实际示例，则可以跳过以下步骤。否则，让我们加载打包的甲基化数据(bigBinDT2)，这是通过上一步得到的。

# MIRA包数据(bigBinDT2)

我们使用与前面相同的过程添加注释数据:

#数据。table函数用于将注释对象中的信息添加到bigBinDT2 setkey(exampleAnnoDT2, sampleName) setkey(bigBinDT2, sampleName) bigBinDT2 <- merge(bigBinDT2, exampleAnnoDT2, all.x=TRUE)

让我们看看MIRA的配置文件:

plotMIRAProfiles (binnedRegDT = bigBinDT2)

接下来，我们规范化归档的甲基化数据(MIRA配置文件)。通过从每个样本/区域集组合中减去最小甲基化bin值进行归一化处理，通常会使中心MIRA剖面得到相似的甲基化值，使MIRA得分更多地与剖面的形状有关，而较少地与平均甲基化有关。我们用data.table语法。

bigBinDT2[， methylProp:= methylProp - min(methylProp) + .05, by=。(featureID sampleName))

normPlot <- plotMIRAProfiles(binnedRegDT=bigBinDT2)

现在我们可以计算所有样本的MIRA分数calcMIRAScore函数。

- calcMIRAScore(bigBinDT2, shoulderShift="auto"， regionSetIDColName="featureID"， sampleIDColName="sampleName") head(sampleScores)

## 1: Ewing_Specific EWS_T1 2.180698 ## 2: Ewing_Specific EWS_T10 2.114649 ## 3: Ewing_Specific EWS_T11 2.064900 ## 4: Ewing_Specific EWS_T12 2.290313 ## 5: Ewing_Specific EWS_T2 1.987881 ## 6: Ewing_Specific EWS_T3 2.230074

最后，让我们重新添加注释数据并查看分数!

setkey(exampleAnnoDT2, sampleName) setkey(sampleScores, sampleName) sampleScores <- merge(sampleScores, exampleAnnoDT2, all.x=TRUE) plotMIRAScores(sampleScores)

##警告:不赞成指定' guide = FALSE '来删除参考线。请##使用' guide = "none"代替。