从配对端ChIP-Seq数据中探索蛋白质- dna结合位点

打开脚本

这个例子显示了如何执行转录因子的全基因组分析拟南芥(Thale Cress)模式生物。

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

简介

ChIP-Seq是一种用于识别与特定DNA位点相互作用的转录因子的技术。第一次染色质免疫沉淀使用结合特定蛋白质的抗体来丰富dna -蛋白质复合物。然后，所有产生的片段都使用高通量测序进行处理。测序片段被映射回参考基因组。通过检查过度代表的区域，有可能标记dna -蛋白质相互作用的基因组位置。

在本例中，短读取由配对端Illumina®平台生成。每个片段都是由成功映射的两个短读重建的，这样就可以计算出片段的确切长度。使用来自序列读取的配对端信息最大限度地提高了预测dna -蛋白质结合位点的准确性。

数据集

本例探讨Wang生成的配对端ChIP-Seq数据出版社。[1]使用Illumina®平台。数据集已提交到基因表达Omnibus存储库，登录号为GSM424618。未映射的配对端读可以从NCBI FTP站点．

这个例子假设你:

(1)下载包含未映射短读的数据，并使用NCBI SRA工具包．

(2)通过使用BWA[2]、Bowtie或SSAHA2(这是[1]作者使用的映射器)将短读序列映射到Thale Cress参考基因组，从而生成SAM格式的文件，并且，

(3)先按引用名，再按基因组位置对SAM格式化文件进行排序。

对于此示例的已发布版本，使用BWA映射器[2]映射了8,655,859个对端短读。BWA生成了一个SAM格式的文件(aratha.sam)，有17,311,718条记录(8,655,859 x 2)。随机选择重复命中，只报告了一次命中，但映射质量较低。在加载到MATLAB之前，使用SAMtools[3]将SAM文件排序并转换为BAM格式的文件。

示例的最后一部分还假设您下载了Thale Cress模式生物(包括五条染色体)的参考基因组。取消注释以下代码行，从NCBI存储库下载引用:

% getgenbank(‘NC_003070’,‘FileFormat’,‘fasta’,‘去整理’,‘ach1.fasta’);% getgenbank(‘NC_003071’,‘FileFormat’,‘fasta’,‘去整理’,‘ach2.fasta’);% getgenbank(‘NC_003074’,‘FileFormat’,‘fasta’,‘去整理’,‘ach3.fasta’);% getgenbank(‘NC_003075’,‘FileFormat’,‘fasta’,‘去整理’,‘ach4.fasta’);% getgenbank(‘NC_003076’,‘FileFormat’,‘fasta’,‘去整理’,‘ach5.fasta’);

创建一个MATLAB®接口到BAM格式化文件

为了创建局部对齐并查看覆盖范围，我们需要构造一个BioMap．BioMap有一个接口，提供对存储在BAM格式化文件中的映射短读取的直接访问，从而最大限度地减少实际加载到工作区的数据量。创建一个BioMap访问BAM格式文件中映射的所有短读。

bm =生物地图(“aratha.bam”）

bm = BioMap与属性:SequenceDictionary: {5x1 cell}引用:[14637324x1文件索引属性]签名:[14637324x1文件索引属性]开始:[14637324x1文件索引属性]MappingQuality: [14637324x1文件索引属性]标志:[14637324x1文件索引属性]MatePosition: [14637324x1文件索引属性]质量:[14637324x1文件索引属性]序列:[14637324x1文件索引属性]头:[14637324x1文件索引属性]NSeqs: 14637324名称:"

使用getSummary方法获取现有引用的列表以及映射到每个引用的短读的实际数量。

getSummary (bm)

名称:" Container_Type: '数据是文件索引。'Total_Number_of_Sequences: 14637324 Number_of_References_in_Dictionary: 5 Number_of_Sequences Genomic_Range Chr1 3151847 1 30427671 Chr2 3080417 1000 19698292 Chr3 3062917 94 23459782 Chr4 2218868 1029 18585050 Chr5 3123275 11 26975502

这个例子的剩余部分集中在对五条染色体中的一条进行分析，Chr1．创建一个新的BioMap通过设置第一个染色体的子集来访问映射到第一个染色体上的短读。

bm1 = get子集(bm，“SelectReference”，“Chr1”）

bm1 = BioMap与属性:SequenceDictionary: 'Chr1'引用:[3151847x1文件索引属性]签名:[3151847x1文件索引属性]开始:[3151847x1文件索引属性]MappingQuality: [3151847x1文件索引属性]标志:[3151847x1文件索引属性]MatePosition: [3151847x1文件索引属性]质量:[3151847x1文件索引属性]序列:[3151847x1文件索引属性]头:[3151847x1文件索引属性]NSeqs: 315184747名称:"

通过访问映射短读的开始和停止位置，可以获得基因组范围。

x1 = min(getStart(bm1)) x2 = max(getStop(bm1))

X1 = uint32 1 x2 = uint32 30427671

探索不同分辨率下的覆盖范围金宝搏官方网站

为了探索整个染色体范围的覆盖率，需要一种分箱算法。的getBaseCoverage方法根据有效对齐生成覆盖信号。它还允许您指定一个bin宽度来控制输出信号的大小(或分辨率)。然而，内部计算仍然在碱基对(bp)分辨率下进行。这意味着尽管设置了较大的宾尺寸，但仍然可以观察到覆盖信号中的窄峰。一旦绘制了覆盖信号，您可以在使用工具提示时编程图形的数据光标来显示基因组位置。您可以缩放和平移该图，以确定ChIP-Seq峰值的位置和高度。

[cov,bin] = getBaseCoverage(bm1,x1,x2，“binWidth”, 1000,“binType”，“马克斯”）;图绘图(bin,cov)轴([x1,x2,0,100])%设置轴限制fixGenomicPositionLabels%格式标记标签和添加数据游标包含(“基本立场”) ylabel (“深度”)标题(“1号染色体的覆盖率”）

也可以在bp分辨率(也称为bp分辨率)下探索覆盖信号堆积配置文件)。探索数据中在位置4598837处观测到的一个大峰值。

P1 = 4598837-1000;P2 = 4598837+1000;(p1:p2,getBaseCoverage(bm1,p1,p2)) xlim([p1,p2])%设置x轴限制fixGenomicPositionLabels%格式标记标签和添加数据游标包含(“基本立场”) ylabel (“深度”)标题(“1号染色体的覆盖率”）

用工件识别和过滤区域

观察位置4599029和4599145之间覆盖深度为800+的大峰值。研究这些读数是如何与参考染色体对齐的。您可以检索这些读取的子集，以满足25的覆盖深度，因为这足以理解该区域中正在发生的事情。使用getIndex获取该子集的索引。然后使用getCompactAlignment以显示短读的相应多次对齐。

i = getIndex(bm1,4599029,4599145，“深度”25);bmx = get子集(bm1,i，“inmemory”假)getCompactAlignment (bmx, 4599029, 4599145)

bmx = BioMap与属性:SequenceDictionary: 'Chr1'引用:[62x1文件索引属性]签名:[62x1文件索引属性]开始:[62x1文件索引属性]MappingQuality: [62x1文件索引属性]标志:[62x1文件索引属性]MatePosition: [62x1文件索引属性]质量:[62x1文件索引属性]序列:[62x1文件索引属性]头:[62x1文件索引属性]NSeqs: 62名称:“俺们= 35 x117 char数组AGTT AATCAAATAGAAAGCCCCGAGGGCGCCATATCCTAGGCGC AAACTATGTGATTGAATAAATCCTCCTCTATCTGTTGCGG GAGGACTCCTTCTCCTTCCCCTTTTGG”“AGTGC TCAAATAGAAAGCCCCGAGGGCGCCATATTCTAGGAGCCC GAATAAATCCTCCTCTATCTGTTGCGGGTCGAGGACTCCT CTCCTGCCCCTTTTGG”“AGTTCAA CCCGAGGGCGCCATATTCTAGGAGCCCAAACTATGTGATT TATCTGTTGCGGGTCGAGGACTCCTTCTCCTTCCCCTTCT”“AGTTCAATCAAATAGAAAGC TTCTAGGAGCCCAAACTATGTGATTGAATAAATCCTCCTC AGGACTCCTTCTCCTTCCCCTTTTGG”“AGTT AAGGAGCCCAAAATATGTGATTGAATAAATCCACCTCTATGGACTCCTTCTCCTTCCCCTTTTGG“AGTACAATCAAATAGAAAGCCCCGAGGGCGCCATA TAGGAGCCCAAACTATGTGATTGAATAAATCCTCCTCTAT CCTTCACCTTCCCCTTTTGG“CGTACAATCAAATAGAAAGCCCCGAGGGCGCCATATTC GGAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCT TTCCCCTTTTGG“CGTACAATCAAATAGAAAGCCCCGAGGGCGCCATATTC GGAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCT“CGTACAATCAAATAGAAAGCCCCGAGGGCGCCATATTC GGAGCCCAAGCTATGTGATTGAATAAATCCTCCTCTATCT“CGTACAATCAAATAGAAAGCCCCGAGGGCGCCATATTC GGAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCT”“AGTTCAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTA GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG‘GATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTA GAGCCCAAACTATGTGATTGAATAAATCTTCCTCTATCTG‘GATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTA GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG”“GATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTA GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG”“GATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTA GAGCCCAAATTATGTGATTGAATAAATCCTCCTCTATCTG”ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTGTTG“ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG CACAAACTATGTGATTGAATAAATCCTCCTCTATCTGTTG“ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG CCAAACTATGTGATTGAATAAATCCTCCTCTATCTGTTGC“ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG“ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTCG“ATACAATCAAATAGAAAGCCCCGGGGGCGCCATATTCTAG“ATTGAGTCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG“ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG“ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG ' 'ATACAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAG“CAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAGGAG“CAATCAAATAGAAAGCCCCGAGGGCGCCATATTCTAGGAG“TAGGAGCCCAAACTATGTGATTGAATAAATCCTCCTCTAT“TAGGAGCCCAAACTATGCCATTGAATAAATCCTCCGCTAT“GGAGCCCAAGCTATGTGATTGAATAAATCCTCCTCTATCT“GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG“GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG“GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG“GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG ' ' GAGCCCAAACTATGTGATTGAATAAATCCTCCTCTATCTG '

除了在视觉上确认对齐之外，您还可以探索该区域中所有短读取的映射质量，因为这可能会提示潜在的问题。在这种情况下，只有不到1%的短读具有60的Phred质量，这表明作图者很可能在参考基因组中发现了多个命中，因此分配了较低的作图质量。

figure i = getIndex(bm1,4599029,4599145);嘘(双(getMappingQuality (bm1, i)))标题(' 4599029和4599145之间的读取的映射质量')包含(“Phred质量评分”) ylabel (“读取数”）

该数据集中的大多数大峰值是由于卫星重复区域或由于其接近着丝粒[4]而出现的，并显示出与刚才探索的例子相似的特征。您可以使用相同的程序探索其他具有较大峰值的区域。

为了防止这些有问题的区域，使用了两种技术。首先，假定所提供的数据集使用成对端测序，通过删除未以适当的对对齐的读取，可以减少潜在的对齐器错误或歧义的数量。您可以通过探索国旗字段的SAM格式文件，其中第二个低有效位用于指示短读是否映射为正确的对。

i = find(bitget(getFlag(bm1)，2));bm1_filtered = get子集(bm1,i)

bm1_filtered = BioMap with properties: SequenceDictionary: 'Chr1'引用:[3040724x1文件索引属性]签名:[3040724x1文件索引属性]开始:[3040724x1文件索引属性]MappingQuality: [3040724x1文件索引属性]标志:[3040724x1文件索引属性]MatePosition: [3040724x1文件索引属性]质量:[3040724x1文件索引属性]序列:[3040724x1文件索引属性]头:[3040724x1文件索引属性]NSeqs: 3040724名称:"

第二，只考虑唯一映射的读取。通过查看映射质量，可以检测到同样映射到参考序列的不同区域的读，因为BWA为这种类型的短读分配了较低的映射质量(小于60)。

i = find(getMappingQuality(bm1_filtered)==60);bm1_filtered = get子集(bm1_filtered,i)

bm1_filtered = BioMap with properties: SequenceDictionary: 'Chr1'引用:[2313252x1文件索引属性]签名:[2313252x1文件索引属性]开始:[2313252x1文件索引属性]MappingQuality: [2313252x1文件索引属性]标志:[2313252x1文件索引属性]MatePosition: [2313252x1文件索引属性]质量:[2313252x1文件索引属性]序列:[2313252x1文件索引属性]头:[2313252x1文件索引属性]NSeqs: 2313252名称:"

再次使用这两种方法可视化过滤后的数据集，整个染色体的粗分辨率为1000 bp，小区域的精细分辨率为20000 bp。大多数由人工制品造成的大山峰已经被移除。

[cov,bin] = getBaseCoverage(bm1_filtered,x1,x2，“binWidth”, 1000,“binType”，“马克斯”）;图绘图(bin,cov)轴([x1,x2,0,100])%设置轴限制fixGenomicPositionLabels%格式标记标签和添加数据游标包含(“基本立场”) ylabel (“深度”)标题(“筛选后1号染色体的覆盖率”) p1 = 24275801-10000;P2 = 24275801+10000;图(p1:p2,getBaseCoverage(bm1_filtered,p1,p2)) xlim([p1,p2])%设置x轴限制fixGenomicPositionLabels%格式标记标签和添加数据游标包含(“基本立场”) ylabel (“深度”)标题(“筛选后1号染色体的覆盖率”）

从pair - end Reads中恢复测序片段

在Wang的论文[1]中，假设成对端测序数据有可能提高DNA相关蛋白染色体结合位点识别的准确性，因为片段长度可以准确地推导出来，而当使用单端测序时，必须采用片段长度的统计近似，并将其模糊地用于所有假定的结合位点。

使用配对的末端读取来重建测序片段。首先，获取每对中正向和反向读取的索引。的第5位捕获该信息国旗字段，根据SAM文件格式。

fow_idx = find(~bitget(getFlag(bm1_filtered)，5));rev_idx = find(bitget(getFlag(bm1_filtered)，5));

sam格式的文件使用相同的头字符串来标识配对。通过对标题字符串，可以确定如何读取短内容BioMap都是成对的。要对标题字符串进行配对，只需将它们按升序排列，并使用排序索引(高频而且人力资源)来链接未排序的头字符串。

[~，hf] = sort(getHeader(bm1_filtered,fow_idx));[~，hr] = sort(getHeader(bm1_filtered,rev_idx));Mate_idx = 0 (numel(fow_idx)，1);Mate_idx (hf) = rev_idx(hr);

使用结果fow_idx而且mate_idx检索配对的变量。例如，检索前10个片段的成对末端读取。

为j = 1:10 disp(getInfo(bm1_filtered, fow_idx(j))) disp(getInfo(bm1_filtered, mate_idx(j))))结束

SRR054715.sra。6849385 163 20 60 40M AACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAA BBBBBBBBBBCBCB?2?BBBBB@7;BBC?7=7?BCC4*)3 SRR054715.sra。6849385 83 229 60 40M cctatttcttgttggttttcttttccttcacttagctatgga 06BBBB= bbbbbbbbbbba6 @@@9<*9BBA@>BBBBB SRR054715.sra。6992346 99 20 60 40M AACCCTAAACCTCTGAATCCTTAATCCCTAAATCCCTAAA =B?BCB=2;BBBBB= b8bbccbbbbbb66bbb =BC8BB SRR054715.sra。6992346 147 239 60 40M gtggttttttttttccttcacttagctatggatggtttatct BBCBB6B?B0B8B<'. bbbbbbbb =BBBBB6BBBBB;*6@ SRR054715.sra。8438570 163 47 60 40M CTAAATCCCTAAATCTTTAAATCCTACATCCATGAATCCC BC=BBBBCBB?==BBB;BB;?BBB8BCB??B-BB<*4?+< bb6bb66bbc ?77年bbcbc@4abb-bbbccbb SRR054715.sra。1676744 83 283 60 40m ttgttattattggatacaagctttgctacgatctacatttg ccb6bbb93 < bbbb >> @b ?<@ bbbbbbbbbbb SRR054715.sra。5658991 163 103 60 40m caaacccgaaaccggtttctctggttgaaactcattgtgt 7? bbbbbb;= bbbb ? 8b; b -; bcb-b <49< 6b8-bb ?+?B SRR054715.sra。5658991 83 311 60 40m gatctacatttgggaatgtgagtctcttattgtaacctta 3，< -bbcbbbbbb ?= bbbba < abbbbbbbbb ? 79bbb ?BB SRR054715.sra。4625439 163 143 60 40m atataatgataattttagcgttttttttgcaattgcttatt bbbbb@，*< 8bbb ++ 2b6b;+ 6b8b;8+ 9bb0，' 9b =。= B SRR054715.sra。4625439 83 347 60 40M CTTAGTGTTGGTTTATCTCAAGAATCTTATTAATTGTTTG +BB8B0BBB?BBBB-BBBB22?BBB-BB6BB-BBBBBB?B SRR054715.sra.1007474 163 210 60 40M ATTTGAGGTCAATACAAATCCTATTTCTTGTGGTTTGCTT BBBBBBBB;.>BB6B6',BBBCBB-08BBBBB;CB9630< SRR054715.sra.1007474 83 408 60 40M TATTGTCATTCTTACTCCTTTGTGGAAATGTTTGTTCTAT BBB@AABBBCCCBBBBBBB=BBBCB8BBBBB=B6BCBB77 SRR054715.sra.7345693 99 213 60 40M TGAGGTCAATACAAATCCTATTTCTTGTGGTTTTCTTTCT B>;>BBB9,<6?@@BBBBBBBBBBBBBB7<9BBBBBB6*' SRR054715.sra.7345693 147 393 60 40M TTATTTTTGGACATTTATTGTCATTCTTACTCCTTTGGGG BB-?+?C@>9BBBBBB6.

使用成对的末端索引来构造一个新的BioMap用最少的信息来表示测序片段。首先，计算插入的大小。

J = getStop(bm1_filtered, fow_idx);K = getStart(bm1_filtered, mate_idx);L = k - j - 1;

通过使用由适当数量的跳过雪茄符号分隔的短读原始签名(N）.

n =数字(L);雪茄= cell(n,1);为I = 1:n支雪茄{I} = sprintf(% dN的L (i));结束雪茄= strcat(getSignature(bm1_filtered, fow_idx)，.．.雪茄,.．.getSignature (bm1_filtered mate_idx));

通过连接成对的末端短读的各自序列来重建片段的序列。

seqs = strcat(getSequence(bm1_filtered, fow_idx)，.．.getSequence (bm1_filtered mate_idx));

计算并绘制碎片大小分布。

J = getStart(bm1_filtered,fow_idx);K = getStop(bm1_filtered,mate_idx);L = k - j + 1;图hist(double(L)，100) title(sprintf(“片段大小分布\n %d配对端片段映射到1号染色体”, n))包含(片段大小的) ylabel (“数”）

构建一个新的BioMap表示测序片段。有了这个，您将能够探索覆盖信号以及片段的局部对齐。

bm1_fragments =生物地图(“序列”seq,“签名”雪茄,“开始”, J)

bm1_fragments = BioMap with properties: SequenceDictionary: {0x1 cell} Reference: {0x1 cell} Signature: {1156626x1 cell} Start: [1156626x1 cell uint32] MappingQuality: [0x1 uint8] Flag: [0x1 uint16] MatePosition: [0x1 uint32] Quality: {0x1 cell} Sequence: {1156626x1 cell} Header: {0x1 cell} NSeqs: 1156626 Name: "

使用片段对齐探索覆盖率

将使用重构片段获得的覆盖信号与使用单个配对端读获得的覆盖信号进行比较。请注意，用峰表示的富集结合位点可以更好地从背景信号中区分出来。

cov_reads = getBaseCoverage(bm1_filtered,x1,x2，“binWidth”, 1000,“binType”，“马克斯”）;[cov_fragments,bin] = getBaseCoverage(bm1_fragments,x1,x2，“binWidth”, 1000,“binType”，“马克斯”）;(bin,cov_reads,bin,cov_fragments) xlim([x1,x2])%设置x轴限制fixGenomicPositionLabels%格式标记标签和添加数据游标包含(“基本立场”) ylabel (“深度”)标题(的报道比较)传说(“短读”，“碎片”）

在bp分辨率下执行相同的比较。在这个数据集中，Wang等人。[1]研究了一个基本的螺旋-环-螺旋(bHLH)转录因子。bHLH蛋白质通常结合到一个一致的序列，称为E-box(CANNTG主题)。使用fastaread要加载参考染色体，请搜索E-box在3'和5'方向上的Motif，然后在覆盖信号上覆盖Motif位置。本例工作在1-200,000区域，但是为了更好地描述细节，数字限制被缩小到3000 bp区域。

P1 = 1;P2 = 200000;cov_reads = getBaseCoverage(bm1_filtered,p1,p2);[cov_fragments,bin] = getBaseCoverage(bm1_fragments,p1,p2);Chr1 = fastaread(“ach1.fasta”）;mp1 = regexp(chr1.Sequence(p1:p2)，“CA . . TG”) + 3 + p1;mp2 = regexp(chr1.Sequence(p1:p2)，“GT . .交流”) + 3 + p1;图案= [mp1 mp2];Figure plot(bin,cov_reads,bin,cov_fragments) hold住在情节([1;1;1]*图案,(0;马克斯(ylim);南),“r”xlim([111000 114000])%设置x轴限制fixGenomicPositionLabels%格式标记标签和添加数据游标包含(“基本立场”) ylabel (“深度”)标题(的报道比较)传说(“短读”，“碎片”，“E-box主题”）

观察到，不可能将覆盖信号中的每个峰值与E-box主题。这是因为测序片段的长度与平均基序距离相当，模糊了接近的峰。画出距离的分布E-box主题网站。

Motif_sep = diff(排序(motifs));图hist(motif_sep(motif_sep<500)，50) title(“相邻E-box图案之间的距离(bp)”)包含(“距离(bp)”) ylabel (“计数”）

发现覆盖信号中的显著峰值

使用函数mspeaks用小波去噪对片段对齐的覆盖信号进行峰值检测。使用高度过滤器过滤假定的ChIP峰，以去除未被考虑的结合过程丰富的峰。

Putative_peaks = mspeaks(bin,cov_fragments，“noiseestimator”, 20岁,.．.“heightfilter”10“showplot”,真正的);持有在传奇(“关闭”)情节([1;1;1]*图案(图案> p1 &主题< p2),(0;马克斯(ylim);南),“r”xlim([111000 114000])%设置x轴限制fixGenomicPositionLabels%格式标记标签和添加数据游标传奇(“碎片覆盖”，“小波去噪覆盖率”，“假定的芯片峰值”，“E-box图案”)包含(“基本立场”) ylabel (“深度”)标题(“ChIP-Seq峰值检测”）

使用knnsearch函数来找到与每个假定的峰值最接近的母题。正如预期的那样，大多数富集的ChIP峰接近anE-box主题[1]。这加强了在最佳分辨率(bp分辨率)下进行峰值检测的重要性，当期望结合位点密度高时，就像在的情况下E-box主题。这个例子也说明了对于这种类型的分析，应该考虑对端测序而不是单端测序[1]。

H = knnsearch(motifs'，putative_peaks(:，1));距离= putative_peaks(:，1)-motifs(h(:))';图hist(距离(abs(距离)<200)，50)“距离最近的E-box Motif为每个检测到的峰值”)包含(“距离(bp)”) ylabel (“计数”）

参考文献

[1]王聪茂，徐洁，张大生，Zoe A Wilson，张大兵。“从配对端ChIP-Seq数据中识别体内蛋白质- dna结合位点的有效方法。”BMC生物信息学11日,没有。1(2010): 81。

[２]李，H.和R.德宾。“快速和准确的短读对齐与Burrows-Wheeler变换。”生物信息学25岁,没有。14(2009年7月15日):1754-60。

[3]李，H.， B. Handsaker, A. Wysoker, T. Fennell, J. Ruan, N. Homer, G. Marth, G. Abecasis, R. Durbin, 1000基因组计划数据处理亚组。序列比对/映射格式和SAMtools生物信息学25岁,没有。16(2009.8.15): 2078-79。

[4]Jothi, R.， S. Cuddapah, A. Barski, K. Cui和K. Zhao。“从ChIP-Seq数据中鉴定体内蛋白质- dna结合位点的全基因组。”核酸研究36岁的没有。16(2008年8月1日):5221-31。

［5］霍夫曼，布拉德·G和史蒂文·J·M·琼斯。利用染色质免疫沉淀偶联流式细胞测序对dna -蛋白质相互作用的全基因组鉴定内分泌学杂志201年,没有。1(2009年4月):1 - 13。

[6]Ramsey, Stephen A.， Theo A. Knijnenburg, Kathleen A. Kennedy, Daniel E. Zak, Mark Gilchrist, Elizabeth S. Gold, Carrie D. Johnson等，“全基因组组蛋白乙酰化数据提高了哺乳动物转录因子结合位点的预测。”生物信息学26日,没有。17(2010年9月1日):2071-75。