Bismark甲基化提取器

描述

以下是控制Bismark甲基化提取器所有选项的简要说明。该脚本读取由Bismark亚硫酸氢盐映射器生成的亚硫酸氢盐读取比对结果文件（BAM/CRAM/SAM格式），并提取单个胞嘧啶的甲基化信息。此信息位于甲基化调用字段中，该字段可包含以下字符：

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
~~~   X   表示CHG上下文的甲基化C                      ~~~
~~~   x   表示CHG上下文的非甲基化C                    ~~~
~~~   H   表示CHH上下文的甲基化C                      ~~~
~~~   h   表示CHH上下文的非甲基化C                    ~~~
~~~   Z   表示CpG上下文的甲基化C                      ~~~
~~~   z   表示CpG上下文的非甲基化C                    ~~~
~~~   U   表示未知上下文（CN或CHN）的甲基化C          ~~~
~~~   u   表示未知上下文（CN或CHN）的非甲基化C        ~~~
~~~   .   表示任何不涉及胞嘧啶的碱基                  ~~~
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

甲基化提取器会输出CpG、CHG和CHH上下文中胞嘧啶的结果文件（这种区分实际上在Bismark本身中就已完成）。由于分析的每个C的甲基化信息可能会生成轻松达到数千万甚至数亿行的文件，文件大小可能会变得非常大，处理起来也更加困难。C甲基化信息还会将胞嘧啶甲基化调用拆分为给定亚硫酸氢盐读取比对的四种可能链之一：

OT：原始顶链
CTOT：与原始顶链互补
OB：原始底链
CTOB：与原始底链互补

因此，默认情况下，每个输入文件会生成12个单独的输出文件（除非指定了--comprehensive，见下文）。输出文件可以导入到基因组查看器（如SeqMonk）中，如果需要，可以重新组合成一个数据组（实际上，除非亚硫酸氢盐读取是在保留方向性的情况下生成的，否则按链特异性方式查看数据没有任何意义）。可以选择跳过链特异性输出文件，在这种情况下，只会生成三个用于CpG、CHG或CHH上下文的输出文件。对于链特异性和综合输出，也可以选择将两个非CpG上下文（CHG和CHH）合并为一个单一的非CpG上下文。

输出文件采用以下格式（制表符分隔）：

<序列ID>    <链>    <染色体>    <位置>    <甲基化调用>

使用方法

bismark_methylation_extractor [选项] <文件名>

参数

<文件名>：以空格分隔的Bismark结果文件列表，格式为SAM，从中提取读取中每个胞嘧啶的甲基化信息。

选项

-s/–single-end：输入文件是由单端读取数据生成的Bismark结果文件。如果既未设置-s也未设置-p，则将自动确定实验类型。

-p/–paired-end：输入文件是由双端读取数据生成的Bismark结果文件。如果既未设置-s也未设置-p，则将自动确定实验类型。

–no_overlap：对于双端读取，理论上读取1和读取2可能会重叠。此选项可避免对重叠的甲基化调用进行两次计分（在该过程中仅使用读取1的甲基化调用，因为历史上读取1的碱基调用质量高于读取2）。虽然此选项可消除测序片段中心甲基化调用过多的偏差，但实际上可能会去除相当比例的数据。对于双端数据，此选项默认开启，但可使用--include_overlap禁用。默认值：开启。

–include_overlap：对于双端数据，无论甲基化调用是否重叠，都将提取所有甲基化调用。默认值：关闭。

–ignore <整数>：在处理甲基化调用字符串时，忽略读取1（或单端比对文件）5’端的前<整数>个碱基对。这可以去除例如每个读取开头的限制酶位点或任何其他偏差来源（如PBAT-Seq数据）。

–ignore_r2 <整数>：仅在处理双端测序结果时，忽略读取2 5’端的前<整数>个碱基对。由于亚硫酸氢盐测序实验中读取2的前几个碱基由于用未甲基化的胞嘧啶对超声处理的片段进行末端修复而表现出严重的非甲基化偏差（见M偏差图），建议在开始下游分析之前去除读取2的前几个碱基对。有关更多详细信息，请参阅Bismark用户指南中关于M偏差图的部分。

–ignore_3prime <整数>：在处理甲基化调用字符串时，忽略读取1（或单端比对文件）3’端的最后<整数>个碱基对。这可以去除读取末端的不必要偏差。

–ignore_3prime_r2 <整数>：仅在处理双端测序结果时，忽略读取2 3’端的最后<整数>个碱基对。这可以去除读取末端的不必要偏差。

–comprehensive：指定此选项将把所有四种可能的链特异性甲基化信息合并到依赖于上下文的输出文件中。默认上下文为：

CpG上下文
CHG上下文
CHH上下文

–merge_non_CpG：在--comprehensive模式下，这将生成两个输出文件；在默认模式下，将生成八个链特异性输出文件，用于以下两种情况的C：

CpG上下文
非CpG上下文

–no_header：抑制所有输出文件中的Bismark版本标题行，以便更方便地进行批量处理。

-o/–output_dir DIR：允许指定不同的输出目录（绝对或相对路径）。如果未明确指定，则输出将写入当前目录。

–samtools_path：Samtools安装的路径，例如/home/user/samtools/。如果Samtools已在PATH中，则无需明确指定。

–gzip：甲基化提取器文件（如CpG_OT_...、CpG_OB_...等）将以GZIP压缩形式写出，以节省磁盘空间。此选项也会传递给全基因组胞嘧啶报告。BedGraph和覆盖度文件默认将以.gz格式写出。

–version：显示版本信息。

-h/–help：显示此帮助文件并退出。

–mbias_only：甲基化提取器将读取整个文件，但仅输出M偏差表和图以及报告（可选），然后退出。默认值：关闭。

–mbias_off：甲基化提取器将像往常一样处理整个文件，但不写出任何M偏差报告。仅建议那些故意想保留早期版本M偏差报告的用户使用。默认值：关闭。

–parallel <整数>：也可以是--multicore <整数>。设置用于甲基化提取过程的核心数。如果系统资源充足，这是加快提取过程的可行选项（我们观察到使用多达10个核心时，速度几乎呈线性增加）。请注意，提取BAM文件并写出.gz输出流的典型过程实际上每个--parallel <整数>指定的值会使用约3个核心（1个用于甲基化提取器本身，1个用于Samtools流，1个用于GZIP流），因此--parallel 10可能会使用系统资源的约30个核心。此选项对bismark2bedGraph或全基因组胞嘧啶报告过程没有影响。

–yacht：此选项（Yet Another Context Hunting Tool）会写出关于甲基化调用所属读取的额外信息，其输出旨在输入到NOMe_filtering脚本中。此选项会写出一个单一的any_C_context文件，其中连续包含读取的所有甲基化调用。其预期用途是单细胞NOMe-Seq数据，因此此选项仅在单端模式下有效（双端读取通常会受到嵌合体问题的影响…）。--yacht将在标准甲基化调用文件中添加三列：

<读取起始位置>
<读取结束位置>
<读取方向>

对于正向读取（+方向），起始位置是最左侧的位置；对于反向读取（-方向），起始位置是最右侧的位置。

BedGraph特定选项

–bedGraph：完成甲基化提取后，将甲基化输出写入排序后的bedGraph文件，该文件报告给定胞嘧啶的位置及其甲基化状态（以百分比表示，详见下文）。甲基化提取器输出会临时拆分为每个染色体一个的临时文件（写入当前目录或使用-o/--output指定的文件夹）；然后使用这些临时文件进行排序，排序后将其删除。默认情况下，仅对CpG上下文中的胞嘧啶进行排序。可以使用--CX_context选项报告所有胞嘧啶，而不考虑序列上下文（这将花费更长时间！）。排序过程中临时文件的默认文件夹是输出目录。bedGraph转换步骤由外部模块bismark2bedGraph执行；此脚本需要与bismark_methylation_extractor位于同一文件夹中。
–zero_based：写出一个额外的覆盖度文件（以.zero.cov结尾），该文件使用0基基因组起始坐标和1基基因组结束坐标（0基，半开区间），类似于bedGraph文件中使用的坐标，而不是始终使用1基坐标。默认值：关闭。
–cutoff [阈值]：在报告核苷酸的甲基化百分比之前，任何甲基化状态（甲基化或非甲基化）必须被观察到的最小次数。默认值：1。
–remove_spaces：将序列ID字段中的空格替换为下划线，以便进行排序。
–CX/–CX_context：排序后的bedGraph输出文件包含实验中覆盖的每个单个胞嘧啶的信息，而不考虑其序列上下文。这适用于正向和反向链。请注意，此选项可能会生成大型临时文件和输出文件，并且排序可能需要很长时间（长达数小时）。默认值：关闭（即默认仅考虑CpG上下文）。
–buffer_size <字符串>：允许在排序甲基化信息时指定主内存排序缓冲区。可以通过附加%指定物理内存的百分比（例如--buffer_size 50%），或者指定1024字节的倍数，例如K表示乘以1024，M表示乘以1048576，依此类推（例如--buffer_size 20G）。有关排序的更多信息，请在命令行中输入info sort。默认值为2G。
–scaffolds/–gazillion：处理包含数十万个甚至数十万个支架/重叠群/染色体的未完成基因组的用户经常会遇到将读取预排序到单个染色体文件时的错误。这些错误是由操作系统在任何时候可以写入的文件句柄数量限制（通常为1024；要查看Linux上的此限制，请输入ulimit -a）引起的。为了绕过打开文件句柄的限制，--scaffolds选项不会将甲基化调用预排序到单个染色体文件中。相反，所有输入文件将临时合并为一个文件（除非只有一个文件），然后使用Unix排序命令按染色体和位置对该文件进行排序。请注意，此选项可能需要很长时间才能完成，具体取决于输入文件的大小以及为该过程分配的内存（见--buffer_size）。不过，它似乎是可行的。
–ample_memory：使用此选项将不使用UNIX的sort命令对染色体位置进行排序，而是使用两个数组对甲基化和非甲基化调用进行排序。这可能会加快非常大文件的排序过程，但代价是更大的内存占用（两个长度与人类最大染色体1（约2.5亿bp）相同的数组将消耗约16GB的RAM）。由于创建和遍历这些数组的开销，对于小文件（数百万个比对），它实际上可能会更慢，我们目前正在测试何时适合使用此选项。请注意，--ample_memory与--scaffolds/--gazillion选项不兼容（因为它需要预先排序的文件）。
–ucsc：写出一个额外的bedGraph文件（以_UCSC.bedGraph.gz结尾），该文件与UCSC基因组浏览器兼容。如果比对是针对Ensembl版本的基因组进行的，此步骤将在染色体名称前加上chr，并将线粒体染色体从MT重命名为chrM。此外，还会写出一个以.chromosome_sizes.txt结尾的制表符分隔文件，以启用诸如bedGraphToBigWig之类的工具。

全基因组胞嘧啶甲基化报告特定选项

–cytosine_report：完成转换为bedGraph后，--cytosine_report选项将生成全基因组所有胞嘧啶的甲基化报告。默认情况下，输出使用1基染色体坐标（可选0基起始坐标），并且仅报告CpG上下文（可选所有胞嘧啶上下文）。输出考虑正向和反向链上的所有C，并报告它们的位置、链、三核苷酸内容和甲基化状态（如果未覆盖，计数为0）。胞嘧啶报告转换步骤由外部模块coverage2cytosine执行；此脚本需要与bismark_methylation_extractor位于同一文件夹中。
–CX/–CX_context：输出文件包含基因组中每个单个胞嘧啶的信息，而不考虑其上下文。这适用于正向和反向链。请注意，对于人类或小鼠等哺乳动物基因组，这将生成超过11亿行的输出文件。默认值：关闭（即默认仅考虑CpG上下文）。
–zero_based：使用0基基因组坐标而不是1基坐标。默认值：关闭。
–genome_folder <路径>：输入要用于提取序列的基因组文件夹（仅全路径）。接受的格式是扩展名为.fa或.fasta的FastA文件。必须指定基因组文件夹路径。
–split_by_chromosome：将输出写入每个染色体的单独文件中，而不是一个单一的输出文件。文件名将包含输入文件名和染色体编号。

输出

bismark_methylation_extractor默认输出格式

<序列ID>  <甲基化状态*>  <染色体>  <起始位置（=结束位置）>  <甲基化调用>

甲基化的胞嘧啶具有+方向。
非甲基化的胞嘧啶具有-方向。

bismark_methylation_extractor使用--yacht（可选）时的输出格式

<序列ID>  <甲基化状态*>  <染色体>  <起始位置（=结束位置）>  <甲基化调用>  <读取起始位置>  <读取结束位置>  <读取方向>

甲基化的胞嘧啶具有+方向。
非甲基化的胞嘧啶具有-方向。

BedGraph输出（可选）格式（制表符分隔；0基起始坐标，1基结束坐标）

track type=bedGraph (标题行)
<染色体>  <起始位置>  <结束位置>  <甲基化百分比>

覆盖度输出格式（制表符分隔，1基基因组坐标；使用--zero_based可获得0基半开区间坐标）

<染色体>  <起始位置>  <结束位置>  <甲基化百分比>  <甲基化计数>  <非甲基化计数>

全基因组胞嘧啶甲基化输出文件格式（制表符分隔）

<染色体>  <位置>  <链>  <甲基化计数>  <非甲基化计数>  <C上下文>  <三核苷酸上下文>