3．1短读的准备和导入

为了描述样本，我们必须提供一个至少有3列的data.frame对象:

• Sample_name:样本的唯一名称，
• file_name: MeDIP seq对齐文件的路径，一般为bam文件
• 分组:将样本分配到分组，如“处理组”和“对照组”。

进一步描述样品的列是可选的，可以在下游分析中考虑:

• 性别:如果提供，则考虑性染色体的数量。
• 输入文件:如果输入库的排序(富集之前)可用，则可以通过指定相应的列在addCNV函数中使用此数据。
• 在统计测试或创建结果表时，可以使用其他分组安排或临床数据。

为了演示QSEA分析步骤，我们使用MEDIPSData包中。该包包含3个NSCLC肿瘤样本和邻近正常组织的染色体20、21和22的对齐MeDIP seq数据。

data(samplesNSCLC, package="MEDIPSData")针织::kable(samples_NSCLC)

路径=系统。file("extdata"， package="MEDIPSData") $file_name=paste0(path，"/"，samples_NSCLC$file_name)

现在，必须加载QSEA包。

库(qsea)

为了访问参考基因组的信息，我们加载预先安装的(见章节)安装) hg19库:

库(BSgenome.Hsapiens.UCSC.hg19)

富集实验的所有相关信息，包括样本信息、基因组读覆盖、CpG密度等参数都存储在一个“qseaSet”对象中。要创建这样一个对象，调用函数" createQseaSet "，它接受以下参数:

• BSgenome:由BSgenome定义的参考基因组名称。(必需)
• Window_size:定义计算短读覆盖的基因组分辨率。(默认:250基地)
• 空空的。select:只处理指定染色体上的数据。(默认值:NULL=所有染色体)
• 区域:如果指定，则只处理指定区域内的数据。

qseaSet = createQseaSet (sampleTable = samples_NSCLC BSgenome = " BSgenome.Hsapiens.UCSC。hg19”,空空的。select=paste0("chr"， 20:22)， window_size=500) ## ====创建qsea set ==== ##对chr20, chr21, chr22的限制性分析##划分BSgenome.Hsapiens.UCSC的选定染色体。500年hg19 windows nt qseaSet # # qseaSet  ## ======================================= ## 6样本:# # 1 t 1 n, 2 t, 2 n, n 3 t, 3 # # 324919地区3 BSgenome.Hsapiens.UCSC.hg19的染色体

现在我们读取对齐文件并计算每个窗口的MeDIP覆盖率。对于这一步，我们必须指定以下参数:

• qs: QSEA集合，在上面的步骤中准备
• uniquePos:如果设置了，起始和结束位置完全相同的片段被认为是PCR重复片段，并用一个有代表性的片段替换
• minMapQual:读取时考虑的最小映射质量。通常，映射质量为0反映了大多数对齐工具中的模糊对齐。
• Paired:如果设置了，则读取被认为是成对结束。在这种情况下，会导入精确的片段大小。
• fragment_size:对于单端测序数据(paired=FALSE)，该参数提供片段大小。

qseaSet=addCoverage(qseaSet, uniquePos=TRUE, paired=TRUE)

# #正常化

QSEA模型可以考虑拷贝数变异(CNV)，并考虑它们对MeDIP读密度的影响。CNV既可以从外部计算并导入QSEA，也可以直接在QSEA内估算。QSEA内部使用生物导体封装HMMcopy
从输入库的测序数据或MeDIP库(参数“MeDIP”=TRUE)估计CNV。在后一种情况下，只考虑不包含CpG二核苷酸的片段。

外部计算的cnv可以通过addCNV函数通过指定" cnv "参数导入:

• cnv:外部计算的cnv作为GRange对象，包含一个具有log2 cnv值的元数据列(相对于正常的Copy Number)。对每个样本进行分析。列的名称必须与示例表中指定的示例名称相匹配。

如果没有提供“cnv”对象，则通过以下参数控制内部cnv分析:

• file_name:指定样本表中用于CNV分析的包含测序比对文件的列。
• window_size:用于CNV分析的窗口大小。
• 配对，“fragment_size”:参见上面的addCoverage步骤
• mu:各个状态下拷贝数的建议中值，参见HMMcopy::HMMsegment
• Normal_idx:样本，假定没有拷贝数变化。默认情况下，在样本表中选择分组为“正常”或“控制”的样本，
• MeDIP:通过CpG甲基化富集测序数据估计CNV(如MeDIP seq)
• plot_dir:可选地在指定目录下生成染色体的CNV图

qseaSet=addCNV(qseaSet, file_name="file_name"，window_size=2e6, paired=TRUE, parallel=FALSE, MeDIP=TRUE)

注意，CNV估计将受益于分析整个基因组。对于这一步，通常不建议将估计限制在单个染色体上。

标度库因子QSEA通过使用m值的修剪平均值方法估计每个样本的有效库大小来解释测序深度和库组成的差异(TMM, Robinson and Oshlack 2010)．或者，QSEA可以使用addLibraryFactors中的" factors "参数导入用户提供的标准化因子。

qseaSet = addLibraryFactors (qseaSet)

估算转化为绝对甲基化值的模型参数

由于MeDIP seq读取覆盖率依赖于片段的CpG密度，我们估计每个基因组窗口每个片段的平均CpG密度。

qseaSet=addPatternDensity(qseaSet， "CG"， name="CpG")

从没有cpg的区域，我们可以估计背景读取的覆盖偏移。

qseaSet = addOffset(qseaSet, enrichment pattern = "CpG") ##选择低CpG密度的窗口用于背景读取估计## 1.511%的窗口具有最多0.01个片段的富集模式密度##，并用于背景读取估计

为了估算相对富集，我们需要建立富集效率模型，富集效率取决于CpG密度。这个模型的参数可以使用addEnrichmentParameters()函数来估计。对于这项任务，我们需要为窗口子集提供已知值，例如通过验证实验或其他研究获得的值。然后QSEA对这些值拟合一个s型函数，以平滑和外推所提供的值。富集参数和估计的背景读数偏移应用于转化为绝对甲基化值。

对于示例数据集，我们可以使用TCGA肺腺癌的甲基化数据(LUAD, Collisson EA 2014)还有鳞状细胞癌(LUSC, Hammerman et al. 2012)研究。为此目的，DNA甲基化数据集已从https://gdc.cancer.gov/，并在500个基窗口中求平均值，并过滤所有样本中高度甲基化的窗口。对象的“tcga_lung”对象中包含该数据的相关染色体MEDIPSData包中。

data(tcga_luad_lusc_450kmeth, package="MEDIPSData") wd=findOverlaps(tcga_luad_lusc_450kmeth, getRegions(qseaSet)， select="first") signal=as.matrix(mcols(tcga_luad_lusc_450kmeth)[，rep(1:2,3)]) qseaSet= addmentparameters (qseaSet, mentpattern ="CpG"， windowIdx=wd, signal=信号)

如果分析的数据集无法获得此类信息，则可以使用粗略估计(“盲校准”)对富集模型进行校准。例如，在成人组织样本中，我们可以假设CpG密度低的区域有高甲基化，而CpG密度高的区域平均甲基化水平呈线性下降。然而，对于不同类型的样本，如细胞系或胚胎细胞，这些假设将是不合适的，并导致错误的结果。

wd=which(getRegions(qseaSet)$CpG_density>1 & getRegions(qseaSet)$CpG_density<15) signal=(15-getRegions(qseaSet)$CpG_density[wd])*.55/15+。25 qseaSet_blind= addmentparameters (qseaSet, mentpattern ="CpG"， windowIdx=wd, signal=signal)

质量控制为了检查MeDIP富集的质量和在前面步骤中所做的模型假设，我们可以检查QSEA估计的模型参数，这些参数描述了信噪比和富集效率。首先，我们检查估计的背景读取的比例:

getOffset(qseaSet, scale="fraction") ## 1T 1N 2T 2N 3T 3N ## 0.13407727 0.10738857 0.11573308 0.09726356 0.10060349 0.10387696

这表明，9%至13%的碎片是由于背景而不是富集。较高的值(例如> 0.9)表明富集效率不足。

为了检查样本CpG密度依赖的富集曲线，由addEnrichmentParameters()估计，QSEA提供了plotemmatrix()函数:

plotEPmatrix(qseaSet) ## 1t# # 1n# # 2t# # 3t# # 3N

在提供的“信号”中观察到的平均富集用黑色表示，用绿色表示拟合的s形函数。剖面平坦的样品可能表明富集效率较低，或与校准数据的一致性较差。

# #探索性分析

为了分析数据的主要特征，并可视化样本之间的关系，QSEA提供了一套用于探索性数据分析的工具。

plotCNV()函数提供了关于CNV的基因组范围的概述。删除用绿色表示，放大用绿色表示，正常拷贝数用蓝色表示。每个样本都表示为一行，这是基于CNV轮廓的相似性排序。

plotCNV (qseaSet)

该图显示，样本2T中chr 20的一半被扩增，而样本1T中chr21的大部分被删除。其他样本，包括3T，拷贝号正常。没有信息的区域(例如，在参考文献中被掩盖的区域)用灰色表示。

为了探索样本之间的关系，主成分图可以是直观的。默认情况下，主成分分析(PCA)是基于日志转换的归一化计数值计算的。为了使用到绝对甲基化值的转换，“norm_method”参数被设置为“beta”。为了获得不同注释级别上的效果的印象，可以将PCA限制在感兴趣的区域上，作为GRanges对象提供。为了演示这个特性，我们导入一个包含不同UCSC注释轨迹的注释对象。为了获得CpG岛启动子区域，我们将CpG岛与转录起始位点(TSS)交叉。

data(annotation, package="MEDIPSData") CGIprom=intersect(ROIs[["CpG岛"]]，ROIs[["TSS"]]，忽略。strand=TRUE) pca_cgi=getPCA(qseaSet, norm_method="beta"， ROIs=CGIprom) col=rep(c("red"， "green")， 3) plotPCA(pca_cgi, bg=col, main="基于CpG岛启动子的PCA图")

在该图中，我们可以看出肿瘤样本和正常样本之间的DNA甲基化是不同的，肿瘤样本是一个异质组，而正常样本紧密地聚集在一起。

QSEA中的差异甲基化分析基于广义线性模型(GLMs)和似然比检验，类似于检测恭顺表达基因的测试，在中实现DESeq2或刨边机．简单地说，我们用负二项分布的均值和色散参数来模拟读计数。对于每个窗口，我们用对数链接函数拟合GLM。对于显著性检验，我们拟合了一个简化模型，并将模型的似然比(LR)与卡方分布进行比较。这些步骤在以下函数中实现:

设计=模型。matrix(~group, getSampleTable(qseaSet)) qseaGLM=fitNBglm(qseaSet, design, norm_method="beta") qseaGLM=addContrast(qseaSet,qseaGLM, coef=2, name="TvN")

首先，使用样本表的group列定义完整模型的模型矩阵(“设计”)。该模型包含两个系数:截距和一个虚拟变量，区分“正常”和“肿瘤”。“fitNBglm”同时拟合了全模型的色散参数和模型系数。在下一步中，“addContrast”拟合了一个模型，减少了第二个成分(因此，只保留了截取)。此外，进行似然比检验。所有结果都存储在一个QSEA GLM对象中。注意，一个GLM对象可以包含多个对比和测试结果(对于更复杂的实验设计)。

注释，探索和导出结果

现在，我们可以创建一个结果表，其中包含原始计数(_reads)估计%甲基化值(_beta)，以及所有依次甲基化窗口估计值的95%可信区间[_betaLB， _betaUB]。“isSignificant”选择低于定义的显著性阈值的窗口。

library(GenomicRanges) sig=isSignificant(qseaGLM, fdr_th=.01) ##选择区域从对比tvn# #选中1006/169189窗口result=makeTable(qseaSet, glm=qseaGLM, groupMeans=getSampleGroups(qseaSet)， keep=sig, annotation=ROIs, norm_method="beta") ##添加测试结果从TvN ##获得原始值的6个样本在1006个窗口##导出beta值…##添加注释##创建表…knitr:: kable(头(结果))

makeTable函数的实现非常通用，它提供了以下选项来指定表的区域和内容:

• Norm_methods:归一化方法名称的向量。有效的名称包括**计数:原始(非规范化)读计数。** nrpkm:每千碱基(基因组窗口)每百万映射读取的CNV归一化读取。** β:转化为绝对甲基化。
• 样本:作为单独列包含在表中的样本名称。
• groupMeans:样本名称的命名列表，汇总为组平均值。
• Glm:可选参数，用于在结果表中包含测试统计信息。
• keep:由isSignificant()获取的窗口id向量
• roi:定义感兴趣区域的GRange对象，结果表将被限制在该区域。
• annotation: GRange对象的命名列表，包含基因组注释。

为了评估基因组注释中不同甲基化区域的富集程度，使用了regionsStats()函数

sigList=list(gain=isSignificant(qseaGLM, fdr_th=.1,direction="gain")， loss=isSignificant(qseaGLM, fdr_th=.1,direction="loss")) ##从对比度tvn# # selected 725/169189窗口中选择区域##从对比度tvn# # selected 2842/169189窗口中选择区域roi_stats=regionStats(qseaSet, subsets=sigList, ROIs=ROIs) ##获得与total qs重叠的数量##获得与gain重叠的数量##获得与loss重叠的数量knitr::kable(roi_stats)

roi_stats_rel=roi_stats[，-1]/roi_stats[，1] x=barplot(t(roi_stats_rel)*100,ylab=" roi的百分比[%]"，names。arg=rep(""， length(ROIs)+1)， beside=TRUE, legend=TRUE, las=2, args.legend=list(x="topleft")， main="特征丰富肿瘤vs正常DNA甲基化")text(x=x[2，]，y=-.15,labels=rownames(roi_stats_rel)， xpd=TRUE, srt=30, cex=1, adj=c(1,0)))

与预期的肿瘤样本一样，低甲基化区域在全基因组范围内更常见，而CpG岛则富含高甲基化。

如果我们对某个特定的基因组区域感兴趣，可以用类似基因组浏览器的方式来描述它:(找到一个更有趣的例子，添加注释……)

plotCoverage(qseaSet, samples=getSampleNames(qseaSet)， chr="chr20"， start=38076001, end=38090000, norm_method="beta"， col=rep(c("red"， "green")， 3)， yoffset=1,space=。1, reorder="clust"， regions=ROIs["TFBS"]，regions_offset=。5, cex =。7) ##选择指定区域##选择指定roi ##在28个窗口中获取6个样本的原始值##导出beta值…##创建表…# # TFBS

#并行化

运行时的很大一部分用于处理对齐文件。方法可以使这些步骤并行化BiocParallel包:

library("BiocParallel") register(MulticoreParam(workers=3)) qseaSet=addCoverage(qseaSet, uniquePos=TRUE, paired=TRUE, parallel=TRUE)

“parallel”参数也可用于addCNV函数。注意，对于文件的并行扫描，读取速度可能是一个限制因素。