Introduction
DMRfinder 是一款用于WGBS的C位点提取、鉴定以及DMR分析的软件，具体过程是：Bismark 软件比对后，用 DMRfinder 提取CpG位点，然后识别差异甲基化区域（call DMR）并与DSS包的call DMR结果做比较。
Quick start
需要给的参数：
-
/path/to/genome：参考基因组绝对路径 -
C1.fastq：对照组的reads1 -
C2.fastq：对照组的reads2 -
E1.fastq：实验组的reads1 -
E2.fastq：实验组的reads2
目的：寻找对照组和实验组之间差异甲基化区域
Step 1: Alignment 比对
调用 Bismark 软件进行比对
$ bismark_genome_preparation /path/to/genome
$ bismark /path/to/genome C1.fastq
$ bismark /path/to/genome C2.fastq
$ bismark /path/to/genome E1.fastq
$ bismark /path/to/genome E2.fastq
Step 2: Extract methylation counts from alignment files
从比对结果文件中提取C位点：
$ samtools view -h C1_bismark_bt2.bam | python extract_CpG_data.py -i - -o C1.cov -v
$ samtools view -h C2_bismark_bt2.bam | python extract_CpG_data.py -i - -o C2.cov -v
$ samtools view -h E1_bismark_bt2.bam | python extract_CpG_data.py -i - -o E1.cov -v
$ samtools view -h E2_bismark_bt2.bam | python extract_CpG_data.py -i - -o E2.cov -v
Step 3: Cluster CpG sites into regions
将CpG位点合并成甲基化区域：
$ python combine_CpG_sites.py -v -o combined.csv C1.cov C2.cov E1.cov E2.cov
Step 4: Test regions for differential methylation
得到差异甲基化区域：
$ Rscript findDMRs.r -i combined.csv -o results.csv -v -n Control,Exptl C1,C2 E1,E2
DMRfinder软件依赖的软件
- DMRfinder
- Bismark
- Bowtie2
- DSS包 这个包需要R和Bioconductor
- Python
- Samtools 可以使 Bismark 比对后的结果文件变成BAM
1、Alignment 比对
使用Bismark软件对甲基化数据进行比对分析。Bismark将测序的结果和参考基因组都进行了C到T和G到A（反向互补）的转化，将转化后的测序结果和基因组分别进行两两比对，将所有的比对结果输出。
Bismark软件会展示三种序列环境下（CpG/CHG/CHH）的比对结果
比对之后，一些研究人员习惯去除PCR duplicates，我们这里不推荐用Bismark的去重脚本deduplicate_bismark完成去重，因为Bismark的去重脚本只是简单粗暴的留下第一个比对上的文件直接删除比对上同样位置的其他文件，没有考虑reads的序列信息和甲基化信息。
（Bismark 的 version 0.16.0版本有些小bug，我们推荐使用 version 0.15.0版本）
2、Extracting methylation counts
DMRfinder 软件的extract_CpG_data.py脚本会将Bismark比对的结果文件转变成一个展示每个CpG site的methylated/unmethylated数量的表格，这个表格按照染色体顺序和位置进行了排序（是按照SAM文件的头文件中参考基因组顺序排列的）位置坐标是从1开始的（1-based）。
该脚本同时还合并了甲基化位点数量，因为它对每个CpG site统计了C和G的总数。
Usage: python extract_CpG_data.py [options] -i <input> -o <output>
python extract_CpG_data.py \
-i Bismark 比对结果文件SAM
-o 输出文件，计算了CpGs序列环境下甲基化和未甲基化的个数，并进行了排序和位点合并
-m 显示CpG序列环境下甲基化位点的最小覆盖度
-s 在输出文件的第三列显示+/-链信息
-n bed文件列出发生甲基化的区域
-b 内存参数
-e 输出文件按照染色体进行排序
输入文件：-i
-i <input> SAM alignment file produced by Bismark (must have
a header, 'XM' methylation strings, and 'XG'
genome designations; can use '-' for stdin)
Bismark的结果文件SAM文件不管是压缩没压缩的格式都可以作为输入文件，如果是BAM文件作为输入文件那必须保证安装了Samtools，否则就转换成SAM：
$ samtools view -h <BAM> | python extract_CpG_data.py -i - -o <output> -v
输出文件：-o
-o <output> Output file listing counts of methylated and
unmethylated CpGs, merged and sorted
输出文件每一行都代表一个CpG位点，每一行有六列中间用tab键分割：
- chrom 染色体名称
- chromStart 胞嘧啶在CpG中的坐标（1-based）
- chromEnd chromStart + 1
- percent 在此位点甲基化水平
- methylated 覆盖到该位点的所有reads中支持甲基化的reads数
- unmethylated 覆盖到该位点的所有reads中不支持甲基化的reads数
-m
-m <int> Minimum coverage (methylation counts) to report a
CpG site (def. 1)
-m参数用来指定每个CpG site至少出现的甲基化数量，低于指定参数将不会输出到结果文件。默认是1，意思是说所有的CpG site中至少出现一个甲基化C位点才会在输出文件中显示
-s
-s Report strand in third column of output
加上这个参数，输出文件的第三列将用链的信息代替chromEnd，因为输出结果是经过合并的，所以链信息总是 +
-n
-n <file> BED file listing regions for which to collect
linked methylation data
bed输入文件应该列出每个感兴趣的基因区域的染色体名称，开始和结束坐标（1-based），唯一的region 名称，用tab键分开，比如：
chrZ 100 200 regionA
输出文件<file>_linked.txt中将会列出bed文件里reads上每一个region的甲基化数据， 在该region至少有一个CpG。用0表示未甲基化，1 表示甲基化，-表示没有数据，比如：
Region: regionA, chrZ:100-200
Sites: 102, 109, 123, 140, 147, 168
00100- read1
001100 read2
--0000 read3
这表明，在regionA有6个CpG sites，对应于3条reads上的6列数字0或1或-。reads1位置123处显示的是1表示该CpG sites发生了甲基化，而位置102, 109, 140, and 147, 显示的是数字0表示CpG sites未甲基化。最后一个CpG sites在位置168上，显示的是 - 表示reads1上没有覆盖此位置。
-b
-b Memory-saving option (with coordinate-sorted SAM)
加上此参数，脚本将会一次处理一条染色体，而不是将所有甲基化数据存入内存。这需要我们预先对SAM文件进行排序（比如使用：samtools sort）。对于大型输入文件我们推荐加上此选项。搭配下面介绍的-e参数，以及combine_CpG_sites.py脚本里面的-b和-e参数。
-e
-e <file> Output file listing ordered chromosomes
选上此参数，输出文件将会按照染色体排好序，同时搭配参数-b时候后按照染色体排好序的输出文件会提供给combine_CpG_sites.py脚本。前提条件是每一个SAM文件是按照相同的顺序排序的。
注意：这个脚本是专门用来合并 CpG甲基化位点个数的。想要得到全面的甲基化信息，比如：CHG/CHH甲基化数量或者M-bias plots，可以使用Bismark 软件中bismark_methylation_extractor脚本。使用该脚本想要输出类似我们脚本extract_CpG_data.py输出的合并结果，必须运行coverage2cytosine脚本（带有--merge_CpG参数）
3、Clustering CpG sites into regions
DMRfinder软件中的combine_CpG_sites.py脚本在单个 CpG 位点获取一个或者多个样本的甲基化数据，并单个输出。按照下面的三步操走，会将相邻的CpG 位点聚集成一个region，并对每个样本区域内的甲基化数据求和：
Usage: python combine_CpG_sites.py [options] -o <output> [<input>]+
[<input>]+ One or more files, each listing methylation counts
for a particular sample
-o <output> Output file listing genomic regions and combined
methylation counts for each sample
Options:
To consider a particular CpG:
-r <int> Min. number of counts at a position (def. 3)
-s <int> Min. number of samples with -r counts (def. 1)
To analyze a region of CpGs:
-d <int> Max. distance between CpG sites (def. 100)
-c <int> Min. number of CpGs in a region (def. 3)
-x <int> Max. length of a region (def. 500)
To report a particular result:
-m <int> Min. total counts in a region (def. 20)
Other:
-f Report methylation fraction for each sample
-y <file> Report clusters of valid CpG sites;
if <file> exists, use these clusters
-b Analyze one chromosome at a time (memory-saving)
-e <file> File listing ordered chromosome names (comma-
separated; used only with -b option)
[<input>]+
[<input>]+ One or more files, each listing methylation counts
for a particular sample
这些文件没有表头，且每行必须有六个用tab键分开的字段，第一，第二，第五和第六列分别是染色体名称，位置，甲基化个数和未甲基化个数（不使用第三和第四列）。这种格式由DMRfinder软件 extract_CpG_data.py脚本生成的（bismark软件的coverage2cytosine脚本加上--merge_CpG参数也会生成相同的格式）
-o
-o <output> Output file listing genomic regions and combined
methylation counts for each sample
输出文件的表头包括：chr，start，end，CpG，以及输入文件中的两列（不加后缀的样品名后跟 -N，样品名称后面跟 -X）。第四列给出了基因区域的位置（start和end给出的就是在该region中CpG sites第一个和最后一个位置信息，从1开始计数的）同时包含了CpG sites的数量。其余列给出了每个样品的甲基化数据，-N列给出了样品在该区域CpG sites的总个数，-X 列给出了其中发生甲基化的个数。
下面详细讲解如何通过这三步合并CpG sites位点为region的，在附录A中给出了示例说明。
第一步，根据 -r 和 -s 两个参数来确定满足覆盖度标准的CpG sites：
To consider a particular CpG:
-r <int> Min. number of counts at a position (def. 3)
-s <int> Min. number of samples with -r counts (def. 1)
在给定的样本中， CpG sites数量达不到-r指定的最小阈值都会被忽略，注意这里说的CpG sites数量包括甲基化的和未甲基化的，甲基化分数并没有考虑进去（the methylation fraction is not taken into account），默认值为3。当设置-r参数为1时意味着每个site都会被分析。
一旦所有样品的单个CpG sites被确定，-s其实就没有多大用了，默认值为1，意味着只要有一个样本满足-r的条件就行，所以就相当于没有过滤任何sites。
第二步，通过三个参数-d、-c和-x将单个的CpG sites聚类成一个region：
To analyze a region of CpGs:
-d <int> Max. distance between CpG sites (def. 100)
-c <int> Min. number of CpGs in a region (def. 3)
-x <int> Max. length of a region (def. 500)
在-d参数所指定距离内的所有CpG sites都会被聚类到一起。例如，默认距离100bp，CpGs 在位置坐标为100/120/200和300处的都会被聚到一类，而在401位置处的不会被聚到这一类中，因为它距离最近的CpG sites超过了100bp。
-c参数指定了每个region至少含有的CpG 个数，否则该region将不会被展示到结果报告中，默认值是3。如果设置为-c 1意味着所有的region都会被reported。
-x参数设置了每个region最大的长度（默认是500bp）。其实该参数并没有严格的执行，超过指定长度的region会按照-c阈值分成子类，如果region通过这种方式无法分割就会原样输出。注意，分割的每个子类必须满足下面的-m参数指定的条件。
第三步，通过 -m 参数实现最后的过滤，该参数指定了每个样本region中总计数的最小值：
To report a particular result:
-m <int> Min. total counts in a region (def. 20)
对于每个region，每个样本中有效的CpG sites不管有没有发生甲基化都会被计算进去，也就是说，total counts in a region = methylated + unmethylated 。任何一个样本只要total counts in a region小于-m参数给定的阈值，在输出文件中-N和-X列都会显示为NA，否则的话，-N列会显示total counts in a region，-X列会显示the methylated sum。
注意，所谓有效的CpG sites就是要对所有的样本都是有效的，即使那些没有满足-r参数的样本。
-f
Other:
-f Report methylation fraction for each sample
使用此参数，每个样品的输出文件将仅包含一列：甲基化百分数methylation fraction in each region，该值是由该region发生甲基化的CpG sites总数除以该region总CpG sites个数得到的，而不是这些CpG sites甲基化分数的平均值。
也就是说，methylation fraction in each region = the sum of methylated counts at the CpG sites within a region / the total sum of counts
-y
-y <file> Report clusters of valid CpG sites;
if <file> exists, use these clusters
如果给指定文件<file>不存在，该参数会创建一个文件来展示聚类情况，每个聚类是一行（每行展示reference name和CpG sites，用逗号分割）
如果指定文件<file>存在，该参数将用文件中定义好的聚类去计算输入文件中methylation counts总数，除了-m参数以外的其他聚类参数(-r, -s, -d, -c, -x)都会被忽略。
对于处理大量数据或者大量不同样本数据来说，这个参数带来了很大的便利。它可以展示样本的子类的聚类，然后用定义好的聚类去所有样本中提取计数结果，这会节约很多时间消耗和资源消耗。
-b
-b Analyze one chromosome at a time (memory-saving)
这个参数可以让脚本一次分析输入文件中的一条染色体，而不是一次将输入文件中的所有数据都加载到内存中，对于大型的输入文件首选此选项。
这里注意一下，染色体顺序一定要与输入文件中的染色体顺序保持一致。
默认情况下，脚本会尝试从输入文件构建正确的染色体顺序，但在少数情况下（比如当多个样本染色体数目不完整时）脚本无法产生正确的顺序并认为输入文件没有排好序，而且一旦遇到这种情况将会特别浪费资源。为了避免的出现这种情况，可以使用下面将要介绍的-e参数。
-e
-e <file> File listing ordered chromosome names (comma-
separated; used only with -b option)
这个参数可以让我们自己给定一个染色体顺序文件，该文件中各染色体用逗号分隔开，或者写成一个list每行是一个染色体。我们输入文件提供的甲基化数量一定要与这里的染色体顺序一致。 Note that one can use the file generated via the -e option of extract_CpG_data.py.
4、Testing regions for differential methylation
DMRfinder 软件中 findDMRs.r 脚本对样品组进行成对儿的比较以发现差异甲基化区域。底层的统计学原理是基于β二项分布模型（beta-binomial hierarchical modeling）和Wald 检验（Wald test），其中 Wald 检验在DSS包中已经测试过了。附录B中提供了findDMRs.r用法的说明性示例。
Usage: Rscript findDMRs.r [options] -i <input> -o <output> \
<groupList1> <groupList2> [...]
-i <input> File listing genomic regions and methylation counts
-o <output> Output file listing methylation results
<groupList> Comma-separated list of sample names (at least two
such lists must be provided)
Options:
-n <str> Comma-separated list of group names
-k <str> Column names of <input> to copy to <output> (comma-
separated; def. "chr, start, end, CpG")
-s <str> Column names of DSS output to include in <output>
(comma-separated; def. "mu, diff, pval")
-c <int> Min. number of CpGs in a region (def. 3)
-d <float> Min. methylation difference between sample groups
([0-1]; def. 0.10)
-p <float> Max. p-value ([0-1]; def. 0.05)
-q <float> Max. q-value ([0-1]; def. 1)
-up Report only regions hypermethylated in later group
-down Report only regions hypomethylated in later group
-t <int> Report regions with at least <int> comparisons
that are significant (def. 1)
-i
-i <input> File listing genomic regions and methylation counts
-i 指定输入文件，必须是由tab键分割并包含：chr，start，end，和CpG。另外对于每一个分析样本必须包含两列：-N后面的样本名称和-X后面的样本名称，这是由combine_CpG_sites.py脚本产生的格式。
-o
-o <output> Output file listing methylation results
输出文件就是成对儿比较后得到的差异甲基化区域，通过特定的参数（-c / -d /-p / -q / -up / -down / -t 下面会介绍）筛选出差异甲基化区域，每个差异甲基化区域包括同样的四列：chr，start，end，和 CpG。除此之外还包括了DSS包的统计结果（methylation levels, differences, and p-values; see -k, -s below）。
<groupList>
<groupList> Comma-separated list of sample names (at least two
such lists must be provided)
### 以逗号分隔的样本名称列表(必须提供至少两个这样的列表)
这是每组要鉴定差异甲基化的样品列表。样品名一定要与其在输入文件中的表头相匹配（not including the -N and -X）。因为是两两比较找差异甲基化，所以必须至少提供两个list。当提供了两个以上的组别时，该脚本就会两两成对儿进行比较。
一个样本不能出现多次，就算是在不同组别中出现的也不行。
-n
-n <str> Comma-separated list of group names
与<groupList>相对应，各组名之间用逗号分割，如果没有提供，该参数默认情况下组名会使用list中的一对儿样品名中间用-连接。
-k
-k <str> Column names of <input> to copy to <output> (comma-
separated; def. "chr, start, end, CpG")
选择输出文件想要包含输入文件的字段，默认是输入文件中的四列（chr, start, end, CpG），使用-k参数指定想要输出的列会加到这四列的后面。比如，如果输入文件包括基因组注释信息，那-k后加上注释列的列标题输出文件中就会包含基因注释列。
-s
-s <str> Column names of DSS output to include in <output>
(comma-separated; def. "mu, diff, pval")
加上-s参数，输出文件中将包含DSS包的分析结果。默认值是：mu(methylation fraction)，diff（methylation difference），pval（p-value），用-s指定DSS输出文件中想要展示在这里的列，将会添加到默认的这三列的后面。DSS (v2.14.0)可选择的值有：phi，diff.se，stat，和fdr。当-q参数被指定，fdr会自动的被选择。
每组中甲基化分数（methylation fraction）以这种方式：group:mu展示。对于组与组之间的比较结果与给定组中的mu值会有稍微的不同，因为这里展示的值是输出文件中的平均值，report中每组其他值（value）也是这样的情况
每组之间的差异甲基化和p-value在report中以这种方式进行展示：group1->group2:diff和group1->group2:pval。两组之间的差异都是group2相对于group1而言的，所以当value小于0说明group2甲基化程度低。由于report中甲基化水平偶尔会出现前后不一致的情况（前面说过这种情况），所以当group1->group2:diff为负数时并不能完全说明group2:mu就一定小于group1:mu（group2的甲基化程度小于group1）
-c
-c <int> Min. number of CpGs in a region (def. 3)
任何一个不包括 CpG sites 区域都会被排除，这个参数等价于 -c in combine_CpG_sites.py
-d
-d <float> Min. methylation difference between sample groups
([0-1]; def. 0.10)
-d参数来指定最小的差异甲基化显著值，默认是0.1，意味着两组之间至少存在10%的差异甲基化（或正数或负数）。当-d参数设定为0时意味着任何一个区域都会被显示，除非NA
-p
-p <float> Max. p-value ([0-1]; def. 0.05)
该参数用来设定最大的差异甲基化p-value值，默认是0.05，说明每组之间的p-value必须小于0.05。当-p参数指定为1时意味着所有region都会被展示，除非显示的是NA。
-q
-q <float> Max. q-value ([0-1]; def. 1)
-q参数用来设置最大的q-value值（校正后的p-value）。默认是1，意味着任何一个region都不会被过滤掉。当使用了此参数，两组之间比较后的输出文件中会自动添加上fdr列。
-up / -down
-up Report only regions hypermethylated in later group
-down Report only regions hypomethylated in later group
这个参数用来限制哪些region被展示。命令行给出的顺序决定了哪个组别时later。
-t
-t <int> Report regions with at least <int> comparisons
that are significant (def. 1)
这个参数用来指定至少有几组比较组合必须满足参数（-d / -p / -q / -up / -down）的阈值才会显示该region，默认是1，意味着只需有1个比较组合被这些参数认定为显著就会显示该region，其他组合中即使不满足这些参数的阈值，该region也会在report中显示。
特别的，当-t参数指定为0时，意味着-d / -p / -q / -up / -down这些参数将会被忽略，无论DSS结果是怎样（包括NA）所有的region都会被显示。
DMRfinder软件所有脚本的公共参数
The following options work for all scripts in DMRfinder:
-h/--help Print the usage message
--version Print the version and copyright
-v Run the script in verbose mode (printing various
updates and statistics to stdout/stderr)
