通过期望最大化调整读取。基于macs(基于模型的芯片序列数据分析)

AREM的Python项目详细描述


arem 1.0.1自述文件,基于macs 1.4.0rc2
时间戳:<;2011-03-01 18:21:42 jake biesinger>;


*简介

/>以及其他DNA结合蛋白。芯片序列数据分析的一个关键步骤是将高通量测序的短读映射到参考基因组,并识别富含短读的峰区域。尽管已经提出了几种芯片序列分析的方法,但大多数现有方法只考虑可以唯一放置在参考基因组中的读操作,因此对于检测重复序列中的lo-
峰的能力较低。在这里,我们介绍了一种利用所有读操作进行芯片序列数据分析的概率ap-
方法,提供了一个真正的基因组范围的结合模式视图。使用与k富集区和空基因组背景相对应的
混合模型对读取进行建模。我们使用最大似然估计富集区的位置,并实现期望最大化(e-m)al-
gorithm,称为arem,以更新读到不同基因组位置的每个
的对齐概率。

.ics.uci.edu/arem


随着测序技术的改进,染色质
免疫沉淀和高通量测序(chip seq)
正逐渐流行于全基因组蛋白质-dna相互作用的研究。为了解决芯片序列分析方法的不足,我们提出了一种基于模型的芯片序列分析(macs)算法,用于识别转录因子结合位点。macs捕获基因组复杂度的影响来评估富含
的芯片区域的意义,并通过结合测序标签
位置和方向的信息来提高
结合位点的空间分辨率。macs可以很容易地单独用于chip seq数据
,或者随着特定性的增加与对照样本一起使用。


原始macs包可在以下位置获得:http://liulab.dfci.harvard.edu/macs/

*install

>请检查发行版中的"install"文件。

*用法

用法:arem<;-t文件>;[-n名称]-g基因组化][选项]

示例:arem-t chip.bam-c control.bam-f bam-g h-n test-w--调用子峰值



arem--根据基于模型的芯片排序(mac)分析,按期望最大化对齐读取


选项:
--版本显示程序的版本号并退出
-h,--help sh显示此帮助消息并退出。
-t文件,--treatment=t文件
chip seq处理文件。必修的。当选择elandmultipet
时,必须提供两个由
逗号分隔的文件,例如
s_1_u eland_u multi.txt,s_1_u eland_u multi.txt
-c cfile,--control=cfile
控制文件。当选择elandmultipet时,
必须提供两个用逗号分隔的文件,例如
s庘2庘u eland庘u multi.txt,s庘2庘u eland庘u multi.txt
-n name,--name=name实验名称,将用于生成输出
文件名。默认值:"na"
-f格式,--format=format
标记文件格式,"auto"、"bed"或"eland"或
"elandmulti"或"elandmultipet"或"elandexport"或
"sam"或"bam"或"bowtie"。默认的自动选项
将允许macs决定文件的格式。如果选择el
和/elandmulti/elandmultipet/elandexport/sam/bam/bowti
e,请检查00readme文件中的定义。默认值:"auto"
--petdist=petdist成对结束标记之间的最佳距离。仅在格式为"elandmultipet"时可用。默认值:200
-g g size,--gsize=gsize
有效基因组大小。它可以是1.0e+9或1000000000,
或快捷键:人类的"hs"(2.7e9)、
(1.87e9)的"mm"、C.elegans的"ce"(9e7)和
果蝇的"dm"(1.2e8),默认值:hs
-s tsize,--tsize=tsize
标签大小。这将覆盖自动检测的标签
大小。默认值:25
--bw=bw带宽。此值仅在构建
移位模型时使用。默认值:300
-p pValue,--pValue=pValue
p峰值检测的值截止。默认值:1e-5
-m mfold,--mfold=mfold
选择mfold范围内高-
背景置信度富集比的区域以构建模型。区域必须低于上限,高于下限。默认值:10,30
--nolambda如果为真,macs将使用固定的背景lambda作为每个峰值区域的本地lambda。通常,macs
计算一个动态局部lambda,以反映由于潜在染色质结构而产生的局部
偏差。
--slocal=small local基对中较小的附近区域,以计算动态lambda。这用于捕捉峰值区域附近的偏差。如果没有控件
数据,则无效。默认值:1000
--llocal=largelocal基对中的较大邻近区域,以计算动态lambda。这用于捕获环绕
偏移。默认值:10000。
--关闭自动是否关闭自动配对模型进程。如果没有
设置,当macs无法构建成对模型时,它将
使用nomodel设置,'--shiftsize'参数
移动和扩展每个标记。默认值:false
--nomodel是否构建移位模型。如果为真,
macs将不会构建模型。默认情况下,它意味着
移位大小=100,尝试将shiftSize设置为更改
。默认值:false
--shiftSize=shiftSize
bp中的任意移位大小。当nomodel为true时,
macs将使用此值作为片段大小的1/2。
默认值:100
--keep dup=keepduplicates
它控制macs对重复标记的行为位置——相同的协调
和相同的链。默认的"auto"选项使
macs使用1e-5作为
p值截止值,基于二次正态分布,在完全相同的位置计算最大标记,
并且"all"选项保留每个标记。
如果给定一个整数,则最多将此数目的标记
保留在同一位置。默认值:auto
--设置为small时,将较大的数据集缩小到较小的
数据集,默认情况下,较小的数据集将
缩小到较大的数据集。默认值:false
-w,-wig是否将每个wigeextend bps处的扩展片段pileup保存到一个wiggle文件中。当--single-
profile打开时,整个基因组只有一个文件被保存。警告:此过程耗费时间/空间!!
-b,--bdg是否将每个bp的扩展片段pileup保存到bedGraph文件中。打开时,-w,
--空格和--调用子峰值将被忽略。当
--启用"单个配置文件"时,只保存整个
基因组的一个文件。警告:此过程耗费时间/空间!!
-s,--single profile设置后,将保存一个wiggle文件用于
处理和输入。默认值:false
--space=space保存摇摆文件的重新排序,默认情况下,
macs将每10 bps保存一次原始标记计数。仅与"--wig"选项一起使用。
--调用子峰值如果设置,macs将通过系统调用调用mali salmon的peaksplitter
soft。如果找不到peaksplitter
,将显示下载
和安装peaksplitter包的说明。-w option
需要打开,-b应该关闭才能工作。
默认值:false
--verbose=verbose set verbose level。0:仅显示关键消息,1:
显示其他警告消息,2:显示进程
信息,3:显示调试消息。默认值:2
--诊断是否生成诊断报告。这要花费9倍的时间。有关
详细信息,请查看00readme文件。默认值:FuxMI= FEMIN用于诊断,min倍富集考虑。
缺省值:0
FEmax = FEMAX用于诊断,最大倍富集考虑。
缺省:最大折叠富集
Fe步骤=FESTEP。对于诊断,折叠浓缩步骤。默认值:20
--没有e m不通过e-m迭代对齐多个读取。多个-
读取概率将基于质量分数或
一致(如果--没有quals)默认值:false
--em converge diff=min_change
在停止e-m步骤之前,迭代之间的最小熵变化。默认值:1e-05
--em min score=min_score
最小富集分数。低于此阈值的窗口
在校准器中看起来都是一样的。默认值:1.5
--em max score=max_score
最大富集分数。高于此阈值的窗口
在校准器中看起来都一样,默认值:否
最大值
--e m show graphs生成e-m的诊断图(需要
matplotlib)。默认值:false
--质量刻度=质量刻度
初始对齐概率由读取
质量和不匹配确定。每个可能的对齐都是
当没有
不匹配时,从所有基上的乘积中分配一个概率
为1-p(readerror_base),或当调用基
不匹配时,分配一个概率统一初始化。读取质量刻度必须
是["auto"、"sanger+33"、"illumina+64"]之一。
默认值:auto
--随机多读通过为每个读随机选择一个对齐将所有多读转换为唯一读。默认值:false
--没有多个丢弃所有具有多个对齐方式的读取
--没有贪婪的调用方使用arem default peak caller而不是贪婪的调用方。这通常会导致较宽、较低浓度的
峰值,特别是在多次读取时。默认值:false
--没有map quals在每个更新步骤中不使用映射概率作为prior;只使用相对丰富。默认值:
false
--prior snp=prior_prob_snp

基因组的任何碱基处发生snp的先验概率。默认值:0.001
--write-read probs写出所有最终读取,包括它们作为bed文件对齐的概率。默认值:false

**参数:

**-t/--treatment filename

这是macs唯一需要的参数。如果格式是
elandmultipet,则用户必须提供两个由
逗号分隔的处理文件,例如s_1_u eland_u multi.txt、s_1_u eland_u multi.txt。

**-c/--control

bed格式或--format选项指定的任何eland输出格式的控件或模拟数据文件。请按照与-t/--处理相同的方向进行操作。

**-n/--name

macs将使用此字符串名
创建输出文件,如"name_peaks.xls"、"name_negative_peaks.xls"、"name_peaks.bed"、"name_summits.bed"、"name_model.r"等。因此,请避免这些文件名与您现有的
文件之间的任何冲突。


**-f/--格式格式


标记文件的格式,可以是"eland"、"bed"、"elandmulti"、"elandexport"、"elandmultipet"(用于对端标记)、"sam"、"bam"或
"bowtie"。默认值为"auto",允许mac自动决定格式。请仅在组合不同格式的文件时使用"auto"。


bed格式在"http://genome.ucsc.edu/faq/faq format\format1"中定义。


>如果格式是eland,则文件必须是eland结果输出文件,
每行只能表示一个标记,字段为:

1。序列名(如果格式不是fasta,则从文件名和行号派生)
2。顺序
3.匹配类型:
nm-未找到匹配。
qc-未完成匹配:qc失败(基本上ns太多)。
rm-未完成匹配:重复屏蔽(如果指定repeatfile.txt,则可能会看到)。
u0-找到的最佳匹配是唯一的完全匹配。
u1-找到的最佳匹配是唯一的1-错误匹配。
u2-找到的最佳匹配是唯一的2错误匹配。
r0-找到多个精确匹配。
r1-找到多个1-错误匹配,没有精确匹配。
r2-找到多个2-错误匹配,没有精确匹配或1-错误匹配。
4。找到的完全匹配项数。
5。找到1个错误匹配项的数目。
6。找到的2个错误匹配数。
只有找到唯一的最佳匹配时,才能看到其余字段(即字段3中的匹配代码以"u"开头)。
7。找到匹配的基因组文件。
8.匹配位置(文件中的基从1开始编号)。
9。匹配方向(F=正向,R=反向)。
10.如何解释read中的n个字符:(""=不适用,"d"=删除,"i"=插入)。
只有在唯一的不精确匹配(即匹配代码是u1或u2)的情况下才能看到其余字段。
11。第一替换错误的位置和类型(例如12a:基12是a,而不是读取的任何内容)。
12。第一个替换错误的位置和类型,如上所述。

如果格式是elandmulti,则文件必须是eland output file from
multiple match mode,每行只能表示一个标记,其中
字段为:
<Br/>1。序列名
2。顺序
3.NM,QC,Rm(如上所述)或以下各项:
4.x:y:z其中x、y和z是找到的精确、单错误和双错误匹配数
5。空白,如果没有找到匹配项,或者如果找到太多匹配项,或者以下内容:
bac廑u plus廑vector.fa:163022R171028F2,e廑coli.fa:3909847R1这表示
与bac廑u plus廑vector有两个匹配项。fa:one in reverse
direction starting at position 160322 with one error,one in
从位置170128开始向前,有两个
错误。如果数据来自对端测序,也有一个与大肠杆菌fa.

您可以将格式
分隔为eland multiple match pair end tags,
然后--treat(和--control,如果需要)参数必须是两个用逗号分隔的文件名。每个文件必须采用上述eland multiple match
格式。例如,

macs14——格式elandmultipet-t s_1_u eland_u multi.txt,s_2_u eland_u multi.txt…


如果使用elandmultipet,则可能需要修改--petdist参数。

如果格式为bam/sam,请检查
(http://samtools.sourceforge.net/samtools.shtml)中的定义。pair end mapping
结果可以保存在单个BAM文件中,如果是,MACS将自动保留left mate(5'end)标记。

如果格式为bowtie,则需要提供后缀为".map"的ASCII bowtie输出文件。请注意,您需要确保
在蝴蝶结输出中,您只为一个
读取保留一个位置。如果您想在
(http://bowtie bio.sourceforge.net/manual.shtml)


查看蝴蝶结手册以获取详细信息,以下是我从上述网页复制的ASCII字符的蝴蝶结输出定义:

1。对齐

2的读的名称。在对齐中读取的方向,-对于反向补码,
+否则

>3。发生对齐的引用序列的名称,或者如果未提供名称,则序号id


4。向前参考链的0基偏移,其中最左边的对齐字符出现


5。读取顺序(如果方向-

6,则反向补全。ascii编码的读取质量(如果方向为-,则相反)。
编码的质量值在phred刻度上,编码的ascii偏移量为33(ascii char!).

7.与此对齐中对齐的引用字符相同的读取与对齐的其他实例数。这不是读取与相同数量的不匹配项对齐的其他位置的数目。此列中的数字
通常不是该数字的良好代理(例如,
此列中的数字可能为"0",而具有相同不匹配数的其他对齐数可能较大)。此列以前称为"保留"。


>8。逗号分隔的不匹配描述符列表。如果对齐中没有不匹配项,则此字段为空。单个描述符的格式偏移量为:reference base>;read base。
偏移量表示为从高质量(5')
读取结束处开始的0基偏移量。


注意:


1)对于床格式,
macs需要第6列钢绞线信息。请注意bed格式的坐标是
基于零且半开放的
(http://genome.ucsc.edu/faq/faq tracks\tracks1)。


计算中包含3个以上的错误。如果原始eland文件中包含一个
标记的多次命中,请删除冗余
以保持该排序标记的最佳命中。

3)对于具有多个重复的实验,建议将多个芯片序列处理文件连接到一个文件中。要执行此操作,请在Unix/Mac或Cygwin(对于Windows OS)下键入:


$cat replicate1.bed replicate2.bed replicate3.bed>;all廑replicates.bed

<4)eland导出格式支持有时可能无法在您的
数据集上运行,因为人们可能会错误地标记第11和第12列。macs
使用第11列作为序列名,序列名应该是染色体名。

**--petdist=petdist

仅在格式为
"elandmultipe"时可用。默认值为200bps。当macs读取5'标记和3'标记的映射位置
时,它将使用
此最佳距离参数来确定它们的最佳配对。一个简单的评分系统使用如下,

score=abs(abs(p5-p3)-200)+e5+e5

,其中p5是5'标签位置之一,e5是5'标签映射位置的
不匹配/错误。p3和e3用于
3'标签。得分最低的配对被认为是最佳配对。该对的5'标记位置保存在模型构建和
峰值调用中。

**-g/--gsize

请根据您的需要分配此参数!

是可映射的基因组大小或有效的基因组大小,即
定义为可测序的基因组大小。由于染色体上的
重复特征,实际可绘制的基因组大小
将小于原始大小,约占基因组大小的90%或70%。对于ucsc human hg18
程序集,建议使用默认的hs--2.7e9。有效基因组大小的所有预编译参数如下:

-g hs=-g 2.7e9
-g mm=-g 1.87e9
-g ce=-g 9e7
-g dm=-g1.2e8

**-s/--tsize

如果不指定,macs将尝试使用输入处理文件的前10个序列来确定标记大小。指定它将覆盖自动确定的标签
大小。

**--bw


您可以将此参数设置为湿实验中预期的声波片段大小
。先前对峰值
检测过程的副作用已经消除。所以这个参数只影响模型的建立。默认值为1e-5。

**-m/--mfold


此参数用于选择mfold范围内
高置信度富集比的区域,建立
模型。区域必须低于上限,高于褶皱富集的
下限。默认值:10,30表示使用所有区域
不太低(>;10)也不太高(<;30)来构建成对峰值
模型。如果macs找不到超过100个区域来构建模型,它将使用--shiftsize参数继续峰值检测。


check related*--off auto*和*--shiftsize*获取详细信息。

**--nolambda

拉姆达这意味着macs将不考虑峰值
候选区域的局部偏压。

**--slocal,--llocal


这两个参数控制在峰值区域周围检查哪两个级别的区域,以计算最大lambda为
局部lambda。默认情况下,对于小的局部
区域(--slocal),macs考虑1000bp;对于大的局部区域(--llocal)
,macs考虑10000bps,这些区域捕获来自长距离效应的偏差,如开放的染色质域。您可以根据您的
项目调整这些选项。请记住,如果区域设置得太小,输入数据中的尖峰可能会杀死有效峰值。

**--关闭自动

如果未设置,当
macs无法构建成对模型时,它将使用nomodel设置,
使用"--shiftsize"参数移动和扩展每个标记。如果设置,
如果配对的峰值模型失败,macs将终止。

**--nomodel


启用此选项时,macs将绕过生成移位模型。

**--shiftsize

例如,如果您的
转录因子的结合区大小为200bp,并且您希望绕过macs构建的
模型,则此参数可以设置为100。此选项只有在设置了--nomodel或macs无法构建
成对峰值模型时才有效。

**--keep dup

默认的
"auto"选项使macs使用1e-5作为pvalue cutoff在完全相同的
位置基于二次正态分布计算最大标记;
,"all"选项保留每个标记。如果给定一个整数,则在
大多数标记数将保持在同一位置。默认值:
auto

**--to small

每个染色体的rmat。gzip的wiggle文件将存储在名为name+"u macs_wiggle/treat"的子目录中,用于处理数据
,名为+"u macs_wiggle/control"的子目录中,用于控制数据。--单个profile
选项可以组合起来为整个基因组生成单个wig文件。

**-b/--bdg

bedGraph文件通常比wiggle文件小得多。但是,由于理论上1bp分辨率数据将被保存,因此该过程将比-w选项稍长一点。bedGraph文件将被gzip压缩并存储在子目录
named name+'''u macs'/bedGraph/treat'中进行处理,并
name+''u macs'/bedGraph/control'中进行控制数据。--单个剖面
选项可以组合起来为
整个基因组生成一个单一的bedgraph文件。

**-s/--单个剖面(以前是单个wig)

y
染色体。gzip的wiggle文件将存储在子目录
named experiment\u name+''u macs\u wiggle'+''u macs\u wiggle/treat/'
+experiment\u name+'treat\u afterfiting\u all.wig.gz'或
'treat\u afterfiting\u all.bdg.gz'中,用于处理数据,以及
experiment\u name+''u macs\u wiggle'+''u macs\u wiggle/control/'
+对于控制数据,实验名+安装完所有.wig.gz的"控制"或
"安装完所有.bdg.gz的"控制"。

**--space=space

您可以使用"--wig"选项来更改它。

如果找不到peaksplitter,将显示用于下载和安装peaksplitter包的指令。
峰值分配器可以细化macs峰,并将宽峰分割为
较小的子峰。有关详细信息,请检查以下URL:

http://www.ebi.ac.uk/bertone/software/peaksplitter_cpp_usage.txt


请注意,如果-b/--bdg处于打开状态,则此选项不起作用。

**--verbose

如果在运行Mac时不想看到任何消息,请将其设置为0。但关键信息永远不会被隐藏。如果您希望
看到丰富的信息,例如每个
染色体有多少峰被调用,您可以将其设置为3或大于3。

**--diag

是选项。这个报告可以帮助你得到一个关于序列饱和的假设。这一功能仅处于β阶段。

**-fe min,--fe max&;-fe step

例如,"--fe min 0--fe max 40--fe step 10"将
让macs选择要考虑的以下折叠富集范围:
[0,10)、[10,20)、[20,30)和[30,40"。


*输出文件


1。name_peaks.xls是一个表格文件,包含有关
称为peaks的信息。您可以在excel中打开它,并使用excel
函数进行排序/筛选。信息包括:染色体名称,
峰的起始位置,
峰的结束位置,峰区的长度,与峰区起始位置相关的峰顶位置,
峰区的标记数,-10*log10(pValue)为峰区(如pValue
=1e-10,则该值应为100),该
区域的褶皱富集度与局部lambda、fdr在
百分比中的随机泊松分布相对应。xls中的坐标是基于1的,这与bed
格式不同。

2。name_peaks.bed是包含峰值
位置的床格式文件。您可以将其加载到ucsc基因组浏览器或affymetrix igb
软件中。name_summits.bed采用bed格式,其中包含每个山峰的峰顶
位置。文件中的第5列是碎片堆积的顶点。如果您想在
绑定站点找到基序,建议使用此文件。

4。name_negative_peaks.xls是一个表格文件,包含有关负峰值的
信息。通过交换芯片序列和控制通道来调用负峰值。name_model.r是一个r脚本,可用于根据数据生成关于模型的pdf
图像。通过以下方式将其加载到r:

$r--vanilla<;name_model.r


,然后将在当前的
目录中生成pdf文件name_model.pdf。注意,绘制此图需要r。

>6。name_treat/control_afterfiting.wig.gz name_macs_wiggle
目录中的文件是wiggle格式的文件,可以导入ucsc
genome browser/gmod/affy igb。.bdg.gz文件是bedgraph
格式,也可以导入到ucsc基因组浏览器中,或者转换成更小的bigwig文件。name_diag.xls是诊断报告。第一列用于不同的
富集范围;第二列为该fc
范围的峰数;第三列后为
采样90%、80%、70%后覆盖的峰百分比……以及总标签的20%。

>8。name_peaks.subpeaks.bed是一个不采用bed
格式的文本文件。此文件由peaksplitter
(<;http://www.ebi.ac.uk/bertone/software/peaksplitter-cpp-usage.txt>;)
在设置--call subpeaks选项时生成。


*其他有用链接

cistrome用于芯片/序列分析的web服务器:http://cistrome.org/ap/


标签:

  • 文件
  • name
  • 标记
  • 模型
  • 基因组
  • 格式
  • 序列
  • 芯片
  • macs
  • 欢迎加入QQ群-->: 979659372 Python中文网_新手群

    推荐PyPI第三方库


    热门话题
    Java swing,将数组列表添加到组件   gwt开发模式运行时的java服务器端测试用例   网格窗格中的java居中按钮和文本(JavaFX)   java heroku Getting at=error code=H13 desc=“连接关闭而无响应”method=GET path=“/api/test”service=28795ms status=503   java在批处理中使用JDBC preparedStatement   按钮文本的java Selenium css选择器   Spring boot未在maven项目中从soap wsdl生成java类。这是我的pom。xml   java如何向POST restful服务发送json数据   servlet线程的java数据库连接?   java ViewPager事件在从用户注销时给我错误   在Java8中检查列表是否同时为空和非空   java相对路径:无法为访问配置文件选择正确的相对路径   用playframework存储二进制密码哈希的java   java第三方JAR创建日志文件   安卓在java中使用常量1的按位移位的优势   Java使用for循环平均从键盘读取的分数集   java我想在同一活动中停止后重新启动相机源   java为什么if-else语句不能按预期工作?输入问题