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

打开脚本

这个例子展示了如何执行转录因子的全基因组分析拟南芥(泰利·克雷斯)模式生物。

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

介绍

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

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

数据集

本例探讨了Wang生成的成对端ChIP-Seq数据出版社。使用Illumina®平台。该数据集已提交给Gene Expression Omnibus数据库，登录号为GSM424618。未映射的对端读取可以从NCBI FTP站点。

这个例子假设你:

(1)下载包含未映射短读的数据，并将其转换为FASTQ格式文件NCBI SRA工具包。

(2)利用BWA[2]、Bowtie或SSAHA2([1]的作者使用的映射器)将短读段映射到Thale Cress参考基因组，生成SAM格式文件;

(3)先按参考名称排序，再按基因组位置排序。

对于本例的已发布版本，使用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 with properties: SequenceDictionary: {5x1 cell}参考:[14637324x1文件索引属性]签名:[14637324x1文件索引属性]开始:[14637324x1文件索引属性]MappingQuality: [14637324x1文件索引属性]标志:[14637324x1文件索引属性]MatePosition: [14637324x1文件索引属性]质量:[14637324x1文件索引属性]序列:[14637324x1文件索引属性]Header: [14637324x1文件索引属性]NSeqs: 146373324名称:"

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

getSummary (bm)

BioMap摘要:名称:" 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 with properties: SequenceDictionary: 'Chr1' Reference: [3151847x1文件索引属性]Signature: [3151847x1文件索引属性]Start: [3151847x1文件索引属性]MappingQuality: [3151847x1文件索引属性]Flag: [3151847x1文件索引属性]MatePosition: [3151847x1文件索引属性]Quality: [3151847x1文件索引属性]Sequence: [3151847x1文件索引属性]Header: [3151847x1文件索引属性]NSeqs: 3151847 Name: "

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

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分辨率下探索覆盖信号(也称为堆积配置文件)。探索数据中观测到的一个位于4598837的大峰。

P1 = 4598837-1000;P2 = 4598837+1000;figure plot(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 = getSubset(bm1,i，“inmemory”假)getCompactAlignment (bmx, 4599029, 4599145)

bmx = BioMap with properties: SequenceDictionary: 'Chr1' Reference: [62x1文件索引属性]Signature: [62x1文件索引属性]Start: [62x1文件索引属性]MappingQuality: [62x1文件索引属性]Flag: [62x1文件索引属性]MatePosition: [62x1文件索引属性]Quality: [62x1文件索引属性]Sequence: [62x1文件索引属性]Header: [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质量，这表明制图者很可能在参考基因组中发现了多个匹配，因此分配了较低的制图质量。

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

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

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

i = find(bitget(getFlag(bm1)，2));bm1_filtered = getSubset(bm1,i)

bm1_filtered = BioMap with properties: SequenceDictionary: 'Chr1' Reference: [3040724x1文件索引属性]Signature: [3040724x1文件索引属性]Start: [3040724x1文件索引属性]MappingQuality: [3040724x1文件索引属性]Flag: [3040724x1文件索引属性]MatePosition: [3040724x1文件索引属性]Quality: [3040724x1文件索引属性]Sequence: [3040724x1文件索引属性]Header: [3040724x1文件索引属性]NSeqs: 3040724 Name: "

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

i = find(getMappingQuality(bm1_filtered)==60);bm1_filtered = getSubset(bm1_filtered,i)

bm1_filtered = BioMap with properties: SequenceDictionary: 'Chr1' Reference: [2313252x1文件索引属性]Signature: [2313252x1文件索引属性]Start: [2313252x1文件索引属性]MappingQuality: [2313252x1文件索引属性]Flag: [2313252x1文件索引属性]MatePosition: [2313252x1文件索引属性]Quality: [2313252x1文件索引属性]Sequence: [2313252x1文件索引属性]Header: [2313252x1文件索引属性]NSeqs: 2313252名称:"

使用这两种方法再次可视化过滤后的数据集，整个染色体的粗分辨率为1000 bp，小区域的精细分辨率为20,000 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;figure plot(p1:p2,getBaseCoverage(bm1_filtered,p1,p2)) xlim([p1,p2])%设置x轴限制fixGenomicPositionLabels%格式化勾号标签并添加数据指针包含(“基本立场”) ylabel (“深度”)标题(筛选后1号染色体的覆盖率）

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

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

利用对端reads重建测序片段。首先，获取每对中正向和反向读取的索引。的第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 aaccctaaacctctgaatccttaatccctaaccctaatccctaabbbbbbbbcbcb ?2?BBBBB@7;BBC?7=7?BCC4*)3 SRR054715.sra。6849385 83 229 60 40M ccattttcttgttgttccttcttttcttagcattgga 06BBBB= bbbbbbbbbbbbbbb6 @@@9<* 9bbba @>BBBBB SRR054715.sra。6992346 99 20 60 40M aaccctaaacctctgaatccttaccctaatccctaaa =B?BCB=2;BBBBB=B8BBCCBBBBBBBC66BBB=BC8BB SRR054715.sra。60 239 147 6992346 40米GTGGTTTTCTTTCCTTCACTTAGCTATGGATGGTTTATCT BBCBB6B ? B0B8B < @ SRR054715.sra .BBBBBBBB = BBBBB6BBBBB; * 6。8438570 163 47 60 40M ctaaatccctaaatcttaatcctacatccatgaatccc BC=BBBBCBB?==BBB;BB;?BBB8BCB??B-BB<*4?+< bb6bb66> ?77年bbcbc@4abb-bbbccbb SRR054715.sra。1676744 83 283 60 40m ttgttattattagatacaagcttttgattagagcattaggattacatatatgccb6bbb93 < bbbb >> @b ?<@ bbbbbbbbbbb SRR054715.sra。5658991 163 103 60 40m caaacccgaaaccgttttctctggttgaaactcattgtgt 7? bbbbbb;= bbbb ? 8b; b -; bbb - b <49< 6b8-bb ?+?B SRR054715.sra。5658991 83 311 60 40m <- bbbbbbbb ?= bbbba < abbbbbbb ? 79bbb ?BB SRR054715.sra。4625439 163 143 60 40m atataatgataattattagtagttttagattgcaattgcttatt 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;

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

n = numel(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)“片段大小分布\n %d配对末端片段映射到1号染色体”, n))包含(片段大小的) ylabel (“数”）

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

bm1_fragments = BioMap“序列”seq,“签名”雪茄,“开始”, J)

bm1_fragments = BioMap with properties: SequenceDictionary: {0x1 cell} Reference: {0x1 cell} Signature: {1156626x1 cell} Start: [1156626x1 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”，“马克斯”）;图plot(bin,cov_reads,bin,cov_fragments) xlim([x1,x2])%设置x轴限制fixGenomicPositionLabels%格式化勾号标签并添加数据指针包含(“基本立场”) ylabel (“深度”)标题(的报道比较)传说(“短读”，“碎片”）

以bp分辨率执行相同的比较。在这个数据集中，Wang等。[1]研究了一个基本的螺旋-环-螺旋(bHLH)转录因子。bHLH蛋白质通常结合在一个称为an的一致序列上E-box(CANNTG主题)。使用fastaread要加载参考染色体，搜索E-box在3'和5'方向上的基序，然后将基序位置叠加在覆盖信号上。这个例子适用于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;motif = [mp1 mp2];图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(sort(motifs));图hist(motif_sep(motif_sep<500)，50) title(“相邻E-box motif之间的距离(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峰都接近于1E-box主题[1]。当结合位点的预期密度很高时，这加强了以尽可能高的分辨率(bp分辨率)进行峰检测的重要性，就像在这种情况下一样E-box主题。这个例子还说明，对于这种类型的分析，应该考虑对端测序而不是单端测序[1]。

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

参考文献

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

[2]李，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岁,没有。(2009年8月15日):2078-79。

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

[5]霍夫曼、布拉德·G和史蒂文·J·M·琼斯。“利用染色质免疫沉淀结合流式细胞测序的dna -蛋白质相互作用的全基因组鉴定。”内分泌学杂志201年,没有。(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。