主要内容

探索DNA甲基化谱的全基因组差异

打开脚本

本例展示了如何使用基因组测序对人类的DNA甲基化进行全基因组分析。

注意:为了增强性能，MathWorks建议您在64位平台上运行此示例，因为内存占用接近2 GB。在32位平台上，如果你收到“内存不足”错误时，请尝试增加操作系统的虚拟内存(或交换空间)或尝试设置3 gbswitch(仅适用于32位Windows®XP)。这些技术将在本文中描述文档．

简介

DNA甲基化是一种表观遗传修饰，在正常和疾病过程中调节基因表达和基因组组织的维持。DNA甲基化可以定义细胞的不同状态，并且在细胞复制过程中是可遗传的。异常的DNA甲基化模式与癌症和肿瘤抑制基因有关。

在本例中，您将探索两种人类癌细胞的DNA甲基化谱:亲本HCT116结肠癌细胞和DICERex5细胞。DICERex5细胞来源于截断DICER1等位基因后的HCT116细胞。Serre等人在[1]中提出利用MBD2蛋白作为甲基CpG结合域来研究DNA甲基化谱，随后使用高通量测序(HTseq)。这种技术通常被称为MBD-Seq。两个样本的两个重复的短读已提交给NCBI的SRA存档由[1]的作者。还有其他技术可用来与HTseq结合来询问CpG位点的DNA甲基化状态，例如MeDIP-seq或限制性内切酶的使用。您还可以按照本例中提供的方法分析这种类型的数据集。

数据集

您可以从NCBI获得四个测序实验的未映射单端reads。使用Illumina®的Genome Analyzer II进行短读。平均插入尺寸为120 bp，短读取长度为36 bp。

这个例子假设你:

(1)下载文件SRR030222.sra，SRR030223.sra，SRR030224.sra而且SRR030225.sra分别包含来自DICERex5样本的两个重复和来自HCT116样本的两个重复的未映射短读NCBI SRA运行选择器并使用NCBI SRA工具包将它们转换为fastq格式的文件。

(2)使用Bowtie[2]算法将短读序列映射到参考人类基因组(NCBI Build 37.5)，生成sam格式的文件。只报告唯一映射的读。

(3)使用SAMtools[3]将SAM格式的文件压缩到BAM中，先按引用名排序，再按基因组位置排序。

本例还假设您下载了参考人类基因组(GRCh37.p5)。您可以使用bowtie-inspect命令直接从领结索引重建人体参考。或者你可以通过取消注释从NCBI存储库下载引用:

% getgenbank(‘NC_000009’,‘FileFormat’,‘fasta’,‘去整理’,‘hsch9.fasta’);

创建bam格式文件的MATLAB®接口

为了探究HCT116样本的信号覆盖范围，您需要构建一个BioMap．BioMap提供了一个接口，可以直接访问存储在bam格式文件中的映射短读，从而最大限度地减少实际加载到内存中的数据量。使用函数baminfo获取现有引用的列表以及映射到每个引用的短读的实际数量。

信息= baminfo(“SRR030224.bam”，“ScanDictionary”,真正的);流(“-35年代% % s \ n”，“参考”，“读取数”）;为i = 1:数字(info.ScannedDictionary) fprintf(“-35年代% % d \ n”,信息。我ScannedDictionary {},.．.info.ScannedDictionaryCount (i));结束
参考数读gi|224589800|ref b| NC_000001.10 b| 205065 gi|224589811|ref| NC_000003.11| 73986 gi|224589816| |NC_000004.11| 96898 gi|224589818|ref|NC_000006.11| 120890 gi|224589819|ref| NC_000008.10| 111229 gi|224589821|ref| NC_000010.10| 104466 gi|224589802|ref|NC_000011.9| NC_000012.11| 87091gi|224589804|ref|NC_000013.10| 53638 gi|224589805|ref| NC_000015.9| 64049 gi|224589807|ref|NC_000016.9| 146868 gi|224589808|ref| NC_000018.9| 60344 gi| 2245898010 |NC_000019.9| 166420 gi|224589810|ref| NC_000020.10| 148950 gi|224589812|ref| NC_000021.8| 310048 gi|224589822|ref|NC_000023.10| 32437 gi|224589822| NC_000023.10| 32437 gi|224589822|ref|NC_000023.10| 32437 gi|224589822| NC_000023.10| NC_001807.4| 1015 Unmapped 6805842

在这个例子中，你将重点分析9号染色体。创建一个BioMap两个HCT116样本重复。

bm_hct116_1 = BioMap(“SRR030224.bam”，“SelectRef”，“gi | 224589821 | ref | NC_000009.11 |”) bm_hct116_2 = BioMap(“SRR030225.bam”，“SelectRef”，“gi | 224589821 | ref | NC_000009.11 |”）
bm_hct116_1 = BioMap与属性:SequenceDictionary: 'gi|224589821|ref|NC_000009.11|'引用:[106189x1文件索引属性]签名:[106189x1文件索引属性]开始:[106189x1文件索引属性]MappingQuality: [106189x1文件索引属性]标志:[106189x1文件索引属性]MatePosition: [106189x1文件索引属性]质量:[106189x1文件索引属性]序列:[106189x1文件索引属性]头:[106189x1文件索引属性]NSeqs: 106189名称:" bm_hct116_2 = BioMap with properties: SequenceDictionary: 'gi|224589821|ref|NC_000009.11|'引用:[107586x1文件索引属性]签名:[107586x1文件索引属性]开始:[107586x1文件索引属性]MappingQuality: [107586x1文件索引属性]标志:[107586x1文件索引属性]MatePosition: [107586x1文件索引属性]质量:[107586x1文件索引属性]序列:[107586x1文件索引属性]头:[107586x1文件索引属性]NSeqs:107586名称:

类提供的分箱算法getBaseCoverage方法，您可以绘制用于初始检查的两个副本的覆盖率。作为参考，您还可以将人类染色体9的表意符号添加到图中chromosomeplot函数。

图ha = gca;持有在N = 141213431;9号染色体长度的%[cov,bin] = getBaseCoverage(bm_hct116_1,1,n，“binWidth”, 100);H1 = plot(bin,cov，“b”）;%表示bm_hct116_1的分箱覆盖率[cov,bin] = getBaseCoverage(bm_hct116_2,1,n，“binWidth”, 100);H2 = plot(bin,cov，‘g’）;%表示bm_hct116_2的分组覆盖率chromosomeplot (“hs_cytoBand.txt”9“AddToPlot”,哈)%沿x轴绘制表意符号轴(ha，[1 n 0 100])%放大y轴fixGenomicPositionLabels (ha)%格式标记标签和添加数据游标传奇((h1 h2),“HCT116-1”，“HCT116-2”，“位置”，“东北”) ylabel (“报道”)标题(“覆盖两个重复的HCT116样本”)图= gcf;fig. position = max(图位置，[0 0 900 0]);%调整大小窗口

因为短读表示DNA的甲基化区域，在对齐覆盖率和DNA甲基化之间存在相关性。观察染色体端粒附近DNA甲基化增加;众所周知，DNA甲基化与端粒维持染色体完整性的作用之间存在关联。在覆盖图中，你还可以看到染色体着丝粒上有一个很长的间隙。这是由于着丝粒中存在重复序列，这阻止了我们将短读序列对齐到该区域的唯一位置。在这个例子中使用的数据集中，只有大约30%的短读被唯一地映射到参考基因组。

CpG岛与DNA甲基化的相关性

DNA甲基化通常发生在CpG二核苷酸中。DNA甲基化模式的改变可导致转录沉默，特别是在基因启动子CpG岛。但是，我们也知道，DNA甲基化可以阻止CTCF结合，并可以沉默miRNA转录等相关功能。一般来说，映射的读取应该最好对齐到CpG丰富的区域。

从参考文件中加载人类9号染色体hs37.fasta．方法从Bowtie索引中恢复了引用bowtie-inspect命令;因此hs37.fasta包含了人类所有的染色体。要只加载9号染色体，可以使用选项中值对BLOCKREAD与fastaread函数。

Chr9 = fastaread(“hs37.fasta”，“blockread”9);chr9。头
ans = 'gi|224589821|ref|NC_000009.11|智人9号染色体，GRCh37主参考装配'

使用cpgisland函数来查找CpG集群。使用CpG岛[4]的标准定义，200个或更多bp岛，CpGobserved/CpGexpected比例为60%或更高，可在9号染色体上发现1682个GpG岛。

cpgi = cpgisland(chr9.Sequence)
cpgi = struct with fields:开始:[10783 29188 73049 73686 113309 114488 116877 117469 117987…停止:[11319 29409 73624 73893 114336 114809 117105 117985 118203…］

使用getCounts计算CpG岛内各基地排列比例的方法。对于样本HCT116的第一次复制，该比例接近45%。

aligned_bases_in_CpG_islands = getCounts(bm_hct116_1,cpgi.Starts,cpgi.Stops，“方法”，“和”) aligned_bases_total = getCounts(bm_hct116_1,1,n，“方法”，“和”) ratio = aligned_bases_in_CpG_islands ./ aligned_bases_total
aligned_bases_in_CpG_islands = 1724363 aligned_bases_total = 3822804 ratio = 0.4511

您可以探索两个样本重复的高分辨率覆盖图，并观察信号如何与CpG岛相关。例如，勘探23,820,000 - 23,830,000 bp之间的区域。这是人类基因ELAVL2的5'区域。

R1 = 23820001;%设置区域限制R2 = 23830000;fhELAVL2 =图;%保留数字句柄，以便以后使用持有在bm_hct116_1的% plot高分辨率覆盖率h1 = plot(r1:r2,getBaseCoverage(bm_hct116_1,r1,r2，“binWidth”1),“b”）;bm_hct116_2的% plot高分辨率覆盖率h2 = plot(r1:r2,getBaseCoverage(bm_hct116_2,r1,r2，“binWidth”1),‘g’）;%标记[r1 r2]区域内的CpG岛屿为i = 1:数字(cpgi.Starts)如果cpgi.Starts(i)>r1 && cpgi.Stops(i)% CpG岛在[r1 r2]里面?Px = [cpgi.]启动([i i]) cpgi。停止([我])];补丁的% x坐标Py = [0 max(ylim) max(ylim) 0];补丁的% y坐标HP = patch(px,py，“r”，“FaceAlpha”、1。“EdgeColor”，“r”，“标签”，“cpgi”）;结束结束轴([r1 r2 0 20])%放大y轴fixGenomicPositionLabels (gca)%格式标记标签和添加数据游标图例([h1 h2 hp]，“HCT116-1”，“HCT116-2”，“CpG岛”) ylabel (“报道”)包含(“9号染色体位置”)标题(“覆盖两个重复的HCT116样本”）

计数数据的统计建模

要找到包含比偶然预期的更多映射读取的区域，可以遵循与Serre等人描述的方法类似的方法。非重叠连续100 bp窗口的计数数被统计建模。

首先，使用getCounts方法计算从每个窗口开始的映射读的数量。在本例中，您使用了一种只考虑每个映射读的起始位置的分箱方法，遵循Serre等人的方法。但是，您也可以使用重叠而且方法中的名称-值对getCounts计算更准确的统计数据。例如，为了获得考虑碱基对分辨率的每个窗口的最大覆盖，设置重叠到1和方法来马克斯．

n = numel(chr9.Sequence);染色体长度百分比W = 1:100:n;%窗口为100 bpcount = getCounts(bm_hct116_1,w,w+99，“独立”,真的,“重叠”，“开始”）;count = getCounts(bm_hct116_2,w,w+99，“独立”,真的,“重叠”，“开始”）;

首先，试着对计数进行建模，假设所有有计数的窗口都具有生物学意义，因此来自相同的分布。利用负双项分布拟合计数数据的模型。

NBP = nbinfit(counts_1);

在经验数据的直方图上绘制拟合模型。

图保存在强调= histc(counts_1,0:100);计算经验分布栏(0:100 emphist. /笔(emphist),“c”，“分组”）%柱状图情节(0:100 nbinpdf(0:100平衡(1),所以(2)),“b”，“线宽”2);%图拟合模型轴([0 50 0 .001])图例(经验分布的，“负二项式拟合”) ylabel (“频率”)包含(“计数”)标题(计数频率，100bp窗口(HCT116-1)）

较差的拟合表明，观测到的分布可能是由于两个模型的混合，一个代表背景，一个代表甲基化DNA窗口中的计数数据。

一个更现实的场景是假设具有少量映射读取的窗口主要是背景(或空模型)。Serre等人假设包含四个或更多读取的100 bp窗口不太可能偶然生成。为了估计一个较好的零模型近似值，您可以将经验分布的左侧主体拟合为截断的负二项分布。函数拟合截断分布大中型企业函数。首先，需要定义一个匿名函数，该函数定义的右截断版本nbinpdf．

Rtnbinpdf = @(x,p1,p2,t) nbinpdf(x,p1,p2) ./ nbincdf(t-1,p1,p2);

使用另一个匿名函数定义拟合函数。

Rtnbinfit = @(x2,t) mle(x2, t)“pdf”@ (x3, p1, p2) rtnbinpdf (x3, p1, p2, t),“开始”nbinfit (x2),“低”[0 0]);

在拟合真实数据之前，让我们用来自已知分布的一些抽样数据来评估拟合过程。

NBP = [0.5 0.2];已知系数%X = nbinrnd(nbp(1)，nbp(2)，10000,1);%随机样本Trun = 6;设置截断阈值Nbphat1 = nbinfit(x);拟合非截断模型到所有数据Nbphat2 = nbinfit(x(x将未截断的模型拟合为截断的数据(错误)Nbphat3 = rtnbinfit(x(x将截断模型拟合为截断数据图保存在强调= histc(x,0:100);计算经验分布栏(0:100 emphist. /笔(emphist),“c”，“分组”）%柱状图H1 = plot(0:100,nbinpdf(0:100,nbphat1(1)，nbphat1(2))，“这”，“线宽”2);H2 = plot(0:100,nbinpdf(0:100,nbphat2(1)，nbphat2(2)))，“r”，“线宽”2);H3 = plot(0:100,nbinpdf(0:100,nbphat3(1)，nbphat3(2))，‘g’，“线宽”2);轴([0 25 0 .2])图例([h1 h2 h3]，“负二项式拟合到所有数据”，.．.“负数二项式拟合截断数据”，.．.“截断的负数二项式拟合截断的数据”) ylabel (“频率”)包含(“计数”）

识别显著甲基化区域

对于HCT116样本的两个重复，使用右截断的负二项分布拟合观察到的零模型rtnbinfit之前定义的匿名函数。

Trun = 4;设置截断阈值(如[1])Pn1 = rtnbinfit(counts_1(counts_1%适合HCT116-1计算Pn2 = rtnbinfit(counts_2(counts_2%适合HCT116-2计算

计算空分布的每个窗口的p值。

Pval1 = 1 - nbincdf(counts_1,pn1(1)，pn1(2));Pval2 = 1 - nbincdf(counts_2,pn2(1)，pn2(2));

方法计算错误发现率mafdr函数。使用名称-值对BHFDR使用线性上升(LSU)程序([6])来计算FDR调整p值。将FDR设置为< 0.01允许您识别显著甲基化的100 bp窗口。

Fdr1 = mafdr(pval1，“bhfdr”,真正的);Fdr2 = mafdr(pval2，“bhfdr”,真正的);W1 = fdr1<.01;%逻辑向量表示HCT116-1的显著窗口W2 = fdr2<.01;%逻辑向量表示HCT116-2的重要窗口W12 = w1 & w2;%逻辑向量，表示两个复制中都有重要的窗口Number_of_sig_windows_HCT116_2 = sum(w2) number_sig_windows_hct116 = sum(w12)
Number_of_sig_windows_HCT116_1 = 1662 number_sig_windows_hct116_2 = 1674 number_sig_windows_hct116 = 1346

总的来说，你们在HCT116样本的两个重复中鉴定出1662个和1674个不重叠的100 bp窗口，这表明存在DNA甲基化的显著证据。在两个复制中都有1346个显著窗口。

例如，再次观察ELAVL2人类基因的启动子区域，您可以观察到在两个样本复制中，多个100 bp的窗口被标记为显著性。

图(fhELAVL2)还原先前绘制的图形情节(w (w1) + 50, counts_1 (w1),“废话”，“HandleVisibility”，“关闭”）%绘制HCT116-1的显著窗口情节(w (w2) + 50, counts_2 (w2),“gs”，“HandleVisibility”，“关闭”）%绘制HCT116-2的显著窗口轴([r1 r2 0 100])“在HCT116样本的两个重复中都有显著的100个bp窗口”）

寻找具有显著甲基化启动子区域的基因

DNA甲基化参与了基因表达的调节。例如，众所周知，高甲基化与几种肿瘤抑制基因的失活有关。你可以在这个数据集中研究基因启动子区域的甲基化。

首先，从Ensembl下载一个标签分离值(TSV)表，其中所有蛋白质编码基因到一个文本文件，ensemblmart_genes_hum37.txt．对于本例，我们使用的是Ensembl发行版64。通过运用的BioMart服务中，您可以选择具有以下属性的表:染色体名称、基因生物型、基因名称、基因起始/结束和链方向。

使用提供的helper函数ensemblmart2gff将下载的TSV文件转换为GFF格式的文件。然后使用GFFAnnotation将文件加载到MATLAB中，并创建一个只存在于9号染色体上的基因的子集。结果在Ensembl数据库中有800个注释的蛋白质编码基因。

GFFfilename = ensemblmart2gff(“ensemblmart_genes_hum37.txt”）;a = GFFAnnotation(GFFfilename) a9 = get子集(a，“参考”，“9”) numGenes = a9。NumEntries
a = GFFAnnotation with properties: FieldNames: {1x9 cell} NumEntries: 21184 a9 = GFFAnnotation with properties: FieldNames: {1x9 cell} NumEntries: 800 numGenes = 800

找到每个基因的启动子区域。在本例中，我们将近端启动子视为-500/100上游区域。

下游= 500;上游= 100;geneDir = strcmp(a9。链,“+”）;%逻辑向量，表示向前方向的股%计算基因正向启动子的起始位置promoterStart(geneDir) = a9.Start(geneDir) -下游;%计算基因正向启动子的末端位置promoterStop(geneDir) = a9.Start(geneDir) + upstream;%反向计算基因启动子的起始位置promoterStart(~geneDir) = a9.Stop(~geneDir) - upstream;%反向计算基因启动子的末端位置promoterStop(~geneDir) = a9.Stop(~geneDir) +下游;

使用一个数据集作为启动子信息的容器，稍后我们可以添加新的列来存储基因计数和p值。

Promoters = dataset({a9。特性,“基因”});启动子。Strand = char(a9.Strand);启动子。Start = promoterStart';启动子。Stop = promoterStop';

通过观察在指定的启动子区域至少有一个碱基对重叠的短读序列的数量，找到在启动子区域具有显著DNA甲基化的基因。

启动子。Counts_1 = getCounts(bm_hct116_1,promoters.Start,promoters.Stop，.．.“重叠”, 1“独立”,真正的);启动子。Counts_2 = getCounts(bm_hct116_2,promoters.Start,promoters.Stop，.．.“重叠”, 1“独立”,真正的);

为每个样本复制拟合一个空分布并计算p值:

Trun = 5;设置截断阈值pn1 = rtnbinfit(promoters.Counts_1(promoters.Counts_1%符合HCT116-1启动子计数pn2 = rtnbinfit(promoters.Counts_2(promoters.Counts_2%符合HCT116-2启动子计数启动子。pval_1 = 1 - nbincdf(promoters.Counts_1,pn1(1)，pn1(2));HCT116-1中每个启动子的p值%启动子。pval_2 = 1 - nbincdf(promoters.Counts_2,pn2(1)，pn2(2));HCT116-2中每个启动子的p值%Number_of_sig_promoters = sum(promoters.pval_1<。01 & promoters.pval_2<.01) ratio_of_sig_methyllated_promoter = number_of_sig_promoter ./numGenes .01 & promoter .pval_2<.01
Number_of_sig_promoters = 74 ratio_of_sig_methyllated_promoters = 0.0925

观察到9号染色体上只有74个基因(800个基因中)存在显著的DNA甲基化区域(两个重复的pval均<0.01)。显示具有最显著甲基化启动子区域的30个基因的报告。

[~，order] = sort(promoters.pval_1.*promoters.pval_2);促销员(订单(1:30)，[1 2 3 4 5 7 6 8])
ans =基因链开始停止Counts_1{‘DMRT3} + 976464 977064 223{‘CNTFR} - 34590021 34590621 219{‘GABBR2} - 101471379 101471979 404{‘CACNA1B} + 101471979 140772341 454{‘BARX1} - 96717554 140772341 264{‘FAM78A} - 134151834 134152434 497{‘FOXB2} + 134152434 79634671 163{‘TLE4} + 79634671 82186788 157{‘ASTN2} - 82186788 120177848 141{‘FOXE1} + 120177848 100615636 149{‘MPDZ} - 13279489 100615636 129{‘PTPRD} - 10612623 10613223 145{‘PALM2-AKAP2} + 10613223112542689 134{‘FAM69B} + 112542689 139607122 112{‘WNK2} + 139607122 95947298 108{‘IGFBPL1} - 95947298 38424944 110{‘AKAP2} + 38424944 112542869 107{‘C9orf4} - 112542869 111930071 102{‘COL5A1} + 111930071 137533720 84{‘LHX3} - 137533720 139097455 74{‘OLFM1} + 139097455 137967368 75{‘NPR2} + 137967368 35792251 68{‘DBC1} - 122131645 35792251 61{‘SOHLH1} - 138591274 138591874 56{‘PIP5K1B} + 138591874 71320675 59{‘PRDM12} + 71320675 133540081 53{‘ELAVL2} -23826235 23826835 50 {'ZFP37'} - 115818939 115819539 59 {'RP11-35N6.1'} + 103790491 103791091 60 {'DMRT2'} + 103790491 103791091 60 {'DMRT2'} + 1049854 1050454 pval_1 Counts_2 pval_2 6.6613e-16 253 5.5511e-16 6.6613e-16 226 5.5511e-16 6.6613e-16 408 5.5511e-16 6.6613e-16 286 5.5511e-16 6.6613e-16 499 5.5511e-16 1.4e-13 165 6.0352e-13 3.5649e-13 151 4.7347e-12 4.3566e-12 163 8.0969e-13 1.2447e-12 133 6.7598e-11 2.3068e -11 127 5.0276e-11 4.1911e-10144 1.3295e-11 7.897e-10 125 2.2131e-10 5.7523e-10 114 1.1364e-09 9.2538e- 09 2.0467e-09 96 1.6795e-08 3.6266e-08 97 1.4452e-08 1.6457e -07 69 1.0074e- 08 1.8093e -07 73 5.4629e-07 1.9575e -06 3.4322e-06 63 2.5345e-06 5.6364e-06 62 2.9575e-06 2.0943e-06 47 3.0746e-05 4.7762e-06 46 3.6016e-05

发现显著甲基化的基因间区

Serre et al.[1]报道，在这些数据集中，大约90%的唯一映射reads位于5'基因启动子区域之外。使用与以前类似的方法，您可以找到具有基因间甲基化区域的基因。为了弥补基因长度的变化，您可以使用逐个碱基计算的最大覆盖，而不是原始的映射短读数。另一种由基因长度规范化计数的方法是设置方法到的名值对rpkm在getCounts函数。

Intergenic = dataset({a9。特性,“基因”});基因间。Strand = char(a9.Strand);基因间。Start = a9.Start;基因间。Stop = a9.Stop;基因间。Counts_1 = getCounts(bm_hct116_1,intergenic.Start,intergen . stop，.．.“重叠”，“全部”，“方法”，“马克斯”，“独立”,真正的);基因间。Counts_2 = getCounts(bm_hct116_2,intergenic.Start,intergenic.Stop，.．.“重叠”，“全部”，“方法”，“马克斯”，“独立”,真正的);Trun = 10;设置截断阈值pn1 = rtnbinfit(intergenic.Counts_1(intergenic.Counts_1%符合HCT116-1基因间计数pn2 = rtnbinfit(intergenic.Counts_2(intergenic.Counts_2%符合HCT116-2基因间计数基因间。pval_1 = 1 - nbincdf(intergenic.Counts_1,pn1(1)，pn1(2));HCT116-1基因间区p值%基因间。pval_2 = 1 - nbincdf(intergenic.Counts_2,pn2(1)，pn2(2));HCT116-2基因间区p值%Number_of_sig_genes = sum(intergenic.pval_1<。01 & intergenic.pval_2<.01)基因甲基化率=基因数量。/numGenes [~，order] = sort(intergenic.pval_1.*intergenic.pval_2);间型(顺序(1:30)，[1 2 3 4 5 7 6 8])
Number_of_sig_genes = 62 Ratio_of_sig_methylated_genes = 0.0775 ans =基因链开始停止Counts_1{‘AL772363.1} - 140762377 140787022 106{‘CACNA1B} + 140787022 141019076 106{‘SUSD1} - 141019076 114937688 88{‘C9orf172} + 114937688 139741797 99{‘NR5A1} - 127243516 139741797 86{‘BARX1} - 96713628 96717654 77{‘KCNT1} + 96717654 138684992 58{‘GABBR2} - 138684992 101471479 65{‘FOXB2} + 79634571 79635869 51{‘NDOR1} + 79635869 140113813 54{‘KIAA1045} + 14011381334984679 50{‘ADAMTSL2} + 34984679 136440641 55{‘PAX5} - 136440641 37034476 48{‘OLFM1} + 37034476 138013025 55{‘PBX3} + 138013025 128729656 45{‘FOXE1} + 128729656 100618986 49{‘MPDZ} - 100618986 13279589 51{‘ASTN2} - 13279589 120177348 43{‘ARRDC1} + 120177348 140509812 49{‘IGFBPL1} - 140509812 38424444 45{‘LHX3} - 38424444 139096955{‘爸爸’}+ 44 139096955 119164601 44{‘CNTFR} - 34551430 119164601 41{‘DMRT3} + 976964 991731 40{‘TUSC1} - 25676396 2567885646 {'ELAVL2'} - 23690102 23826335 35 {'SMARCA2'} + 2015342 2193624 36 {'GAS1'} - 89559279 89562104 34 {'GRIN1'} + 140032842 140063207 36 {'TLE4'} + 82186688 82341658 36 pval_1 Counts_2 pval_2 8.6597e-15 98 1.8763e-14 8.6597e-15 98 1.8763e-14 2.2904e-12 112 7.7716e-16 7.4718e-14 96 3.5749e-14 4.268e-12 90 2.5457e-13 7.0112e-11 62 2.569e-09 2.5424e-08 73 9.5469e- 11 2.5525e-08 3.0134e-07 55 2.5525e-08 6.4307e-085.585e-07 49 1.8188e-07 6.4307e-08 42 1.7861e-06 2.6058e-06 43 1.2894e-06 4.1027e-07 36 1.2564e-05 1.9155e-06 36 1.7377e-05 4.8199e-06 37 9.0815e-06 6.5537e-06 37 1.7377e-05 3.4736e -06 2.2358e-05 40 2.4736e-06 2.2358e-05 38 6.5629e-06 2.2358e-05 37 9.0815e-06

例如，探索BARX1基因的甲基化特征，这是前面列表中第六种具有基因间甲基化的重要基因。GTF格式的文件ensemblmart_barx1.gtf包含该基因的结构信息，从Ensembl使用BioMart服务。

使用GTFAnnotation将结构信息加载到MATLAB中。这个基因有两个注释的转录本。

barx1 = GTFAnnotation(“ensemblmart_barx1.gtf”) getTranscriptNames(barx1)
barx1 =带有属性的GTFAnnotation: FieldNames: {1x11 cell} NumEntries: 18 transcripts = 2x1 cell array {'ENST00000253968'} {'ENST00000401724'}

用碱基对分辨率绘制两个HCT116样本复制的DNA甲基化谱。覆盖CpG岛，并沿着图的底部绘制两个转录本的每个外显子。

range = barx1.getRange;R1 =范围(1)-1000;%设置区域限制R2 = range(2)+1000;图保存在bm_hct116_1的% plot高分辨率覆盖率h1 = plot(r1:r2,getBaseCoverage(bm_hct116_1,r1,r2，“binWidth”1),“b”）;bm_hct116_2的% plot高分辨率覆盖率h2 = plot(r1:r2,getBaseCoverage(bm_hct116_2,r1,r2，“binWidth”1),‘g’）;%标记[r1 r2]区域内的CpG岛屿为i = 1:数字(cpgi.Starts)如果cpgi.Starts(i)>r1 && cpgi.Stops(i)% CpG岛在[r1 r2]里面?Px = [cpgi.]启动([i i]) cpgi。停止([我])];补丁的% x坐标Py = [0 max(ylim) max(ylim) 0];补丁的% y坐标HP = patch(px,py，“r”，“FaceAlpha”、1。“EdgeColor”，“r”，“标签”，“cpgi”）;结束结束在轴的底部标记外显子为i = 1: number (transcripts) exons = get子集(barx1，“记录”成绩单我{},“功能”，“外显子”）;为J = 1:外显子。NumEntriespx = [exons.Start([j j]);exons.Stop([j j])]';补丁的% x坐标Py = [0 1 1 0]-i*2-1;补丁的% y坐标Hq = patch(px,py，“b”，“FaceAlpha”、1。“EdgeColor”，“b”，“标签”，“外显子”）;结束结束轴([r1 r2 -数字(副本)*2- 80])%放大y轴fixGenomicPositionLabels (gca)%格式标记标签和添加数据游标ylabel (“报道”)包含(“9号染色体位置”)标题(BARX1基因的高分辨率覆盖([h1 h2 HP hq]，“HCT116-1”，“HCT116-2”，“CpG岛”，“外显子”，“位置”，“西北”）

观察5'启动子区高度甲基化区域(CpG岛最右侧)。回想一下，这个基因的转录发生在反向链上。更有趣的是，观察高度甲基化的区域，这些区域重叠在两个注释转录本的起始处(两个中间的CpG岛)。

甲基化模式的差异分析

Serre等人的研究还分析了另一种细胞系。在截断DICER1等位基因后，新的细胞(DICERex5)来自相同的HCT116结肠癌细胞。文献[5]报道，在少量基因启动子上存在DNA甲基化的局部性改变。在本例中，采用与亲本HCT116细胞相同的方法，您将在DICERex5细胞的两个样本复制中发现显著的100 bp窗口，然后您将搜索两个细胞系之间的统计显著差异。

辅助函数getWindowCounts捕获与前面类似的步骤，以查找具有重要覆盖率的Windows。getWindowCounts返回每个新复制的包含计数、p值和错误发现率的向量。

bm_dicer_1 = BioMap(“SRR030222.bam”，“SelectRef”，“gi | 224589821 | ref | NC_000009.11 |”）;bm_dicer_2 = BioMap(“SRR030223.bam”，“SelectRef”，“gi | 224589821 | ref | NC_000009.11 |”）;[counts_3,pval3,fdr3] = getWindowCounts(bm_dicer_1,4,w,100);[counts_4,pval4,fdr4] = getWindowCounts(bm_dicer_2,4,w,100);W3 = fdr3<.01;%表示DICERex5_1中重要窗口的逻辑向量W4 = fdr4<.01;%表示DICERex5-2中重要窗口的逻辑向量W34 = w3 & w4;%逻辑向量，表示两个复制中都有重要的窗口Number_of_sig_windows_DICERex5_1 = sum(w3) Number_of_sig_windows_DICERex5_2 = sum(w4) number_of_windows_dicerex5 = sum(w34)
Number_of_sig_windows_DICERex5_1 = 908 Number_of_sig_windows_DICERex5_2 = 1041 Number_of_sig_windows_DICERex5 = 759

要执行差异分析，需要使用至少一个样本(HCT116或DICERex5)中显著的100 bp窗口。

Wd = w34 | w12;%逻辑向量，表示差异分析中包含的Windows计数= [counts_1(wd) counts_2(wd) counts_3(wd) counts_4(wd)];Ws = w(wd);%窗口开始为每一行计数

使用函数manorm将数据规范化。的百分位名称-值对允许您在规范化时过滤掉具有大量计数的窗口，因为这些窗口主要是由于工件造成的，例如参考染色体中的重复区域。

Counts_norm = round(manm(计数，“百分比”, 90)。* 100);

使用函数垫执行双样本t检验，以识别来自两个不同细胞系的不同覆盖的Windows。

pval =垫(counts_norm (:, (1 2)), counts_norm (: [3 - 4]),“引导”,真的,.．.“showhist”,真的,“showplot”,真正的);

创建一个包含25个最重要的差异覆盖窗口的报告。在创建报表时使用helper函数findClosestGene确定窗口是基因间的，基因内的，还是在近端启动子区域。

[~，ord] = sort(pval);流(' windows Pos Type p-value HCT116 DICERex5\n\n'）;为I = 1:25 j = ord(I);[~，msg] = findClosestGene(a9，[ws(j) ws(j)+99]);流('%10d %-25s %7.6f%5d%5d %5d%5d\n'，.．.ws (j)、味精、pval (j), counts_norm (j,:));结束
窗口Pos类型p值HCT116 DICERex5 140311701基因间(EXD3) 0.000026 13 13 104 105 139546501基因内0.001827 21 21 91 93 10901基因内0.002671 258 257 434 428 120176801基因间(ASTN2) 0.002733 266 270 155 155 139914801基因间(ABCA2) 0.002980 64 63 26 25 126128501基因间(CRB2) 0.003193 94 93 129 130 71939501基因启动子(FAM189A2) 0.005549 107 101 00 124461001基因间(DAB2IP) 0.005618 77 76 39 37 140086501基因间(TPRN) 0.006520 47 42 123 124 79637201基因内0.007512 52 51 32 31 136470801基因内0.007512 52 51 32 31 140918001基因间(CACNA1B) 0.008115 176 169 71 68 100615901基因间(FOXE1) 0.008346 262 253 123 118 98221901基因间(PTCH1) 0.009934 26 30 104 99 138730601基因间(CAMSAP1) 0.010273 26 21 97 93 89561701基因间(GAS1) 0.010336 77 76 6 12 977401基因间(PAX5) 0.010559 133 127 207 211 139744401基因间(PHPT1) 0.010869 47 46 32 31 126771301基因内0.011458 43 46 97 93 93922501基因内0.011486 34 34 149 161 94187101基因内0.011507 73 80 6 6 136044401基因内0.011567 39 34 110 105 139611201基因间(FAM69B) 0.011567 39 34 110 105 139716201基因间(C9orf86) 0.011832 73 72 136 130

绘制基因FAM189A2启动子区域的DNA甲基化谱，该启动子区域与前一列表中差异最大。覆盖CpG岛和FAM189A2基因。

range = getRange(get子集(a9，“功能”，“FAM189A2”));R1 =范围(1)-1000;R2 = range(2)+1000;图保存在所有重复的高分辨率覆盖率百分比图h1 = plot(r1:r2,getBaseCoverage(bm_hct116_1,r1,r2，“binWidth”1),“b”）;h2 = plot(r1:r2,getBaseCoverage(bm_hct116_2,r1,r2，“binWidth”1),‘g’）;h3 = plot(r1:r2,getBaseCoverage(bm_dicer_1,r1,r2，“binWidth”1),“r”）;h4 = plot(r1:r2,getBaseCoverage(bm_dicer_2,r1,r2，“binWidth”1),“米”）;%标记[r1 r2]区域内的CpG岛屿为i = 1:数字(cpgi.Starts)如果cpgi.Starts(i)>r1 && cpgi.Stops(i)% CpG岛在[r1 r2]里面?Px = [cpgi.]启动([i i]) cpgi。停止([我])];补丁的% x坐标Py = [0 max(ylim) max(ylim) 0];补丁的% y坐标HP = patch(px,py，“r”，“FaceAlpha”、1。“EdgeColor”，“r”，“标签”，“cpgi”）;结束结束%标记在轴的底部Px = range([1 1 2 2]);Py = [0 1 1 0]-2;Hq = patch(px,py，“b”，“FaceAlpha”、1。“EdgeColor”，“b”，“标签”，“基因”）;轴([r1 r1+4000 -4 30])%被放大fixGenomicPositionLabels (gca)%格式标记标签和添加数据游标ylabel (“报道”)包含(“9号染色体位置”)标题(FAM189A2基因启动子区域的DNA甲基化谱)传说([h1 h2 h3 h4 HP hq]，.．.“HCT116-1”，“HCT116-2”，“DICERex5-1”，“DICERex5-2”，“CpG岛”，“FAM189A2基因”，.．.“位置”，“东北”）

观察到CpG岛在两次DICERex5重复中均明显未甲基化。

参考文献

[1] Serre, D.， Lee, b.h.，和Ting a.h.，“mbd分离的基因组测序提供了一个高通量和全面的人类基因组DNA甲基化调查”，核酸研究，38(2):391- 9,2010。

[2] Langmead, B.， Trapnell, C.， Pop, M.，和Salzberg, S.L.，“超快和内存高效的短DNA序列与人类基因组的对齐”，基因组生物学，10(3):R25, 2009。

李海燕，李海燕，等，“序列比对/图谱(SAM)格式与SAM工具”，生物信息学，25(16):2078- 9,2009。

王晓明，刘志刚，“CpG岛在脊椎动物基因组中的应用”，中国生物工程学报，32(2):366 - 366,1987。

王晓明，王晓明，等，“人体肿瘤细胞中DICER蛋白对CpG岛高甲基化的影响”，中华癌症杂志，38(4):457 - 457,2008。

[6] Benjamini, Y.和Hochberg, Y.“控制错误发现率:一种实用而有效的多重测试方法”，皇家统计学会杂志，57(1):289-300,1995。

生物信息学工具箱文档

金宝app

试试MATLAB、Sim金宝appulink和其他产品下载188bet金宝搏

现在就去审判吧