docs/bdgdiff.md

   1 # bdgdiff
   2
   3 ## Overview
   4 The `bdgdiff` command is part of the MACS3 suite of tools and is used
   5 to call differential peaks from four bedGraph tracks of scores,
   6 including treatment and control track for each condition.
   7
   8
   9 ## Detailed Description
  10
  11 The `bdgdiff` command takes four input bedGraph files (two treatment
  12 and two control files) and produces three output files with
  13 differential peaks called. Users should provide paired four bedgraph
  14 files: for each condition, a treatment pileup signal track in bedGraph
  15 format, and a control lambda track in bedGraph format. This
  16 differential calling can only handle one replicate per condition, so
  17 if you have multiple replicates, you should either combine the
  18 replicates when `callpeak`, or choose other tool that can consider
  19 within-group variation (such as DESeq2 or edgeR). The method we use to
  20 define the differential peaks is based on multiple likelihood tests,
  21 based on the poisson distribution. Suppose that we have two conditions
  22 A and B, the unique binding sites in condition A over condition B
  23 should be *more likely* to be a binding event in treatment A over
  24 treatment B, and also *more likely* to be a real binding site in
  25 condition A while comparing treatment A over control A; the unique
  26 binding sites in condition B is defined in the same way; the common
  27 peaks of both condition should be *more likely* to be a real binding
  28 sites in condition A while comparing treatment A and control A, and in
  29 condition B while comparing treatment B over control B, and also the
  30 likelihood test while comparing treatment A and treatment B can't
  31 decide which condition is stronger.
  32
  33 The likelihood function we used while comparing two conditions: ChIP
  34 (enrichment) or control (chromatin bias) is:
  35
  36 $$ln(LR) = x*(ln(x)-ln(y)) + y - x$$
  37
  38 Here $LR$ is the likelihood ratio, x is the signal (fragment pileup)
  39 we observed in condition 1, and y is the signal in condition
  40 2. And $ln$ is the natural logarithm.
  41
  42 Note: All regions on the same chromosome in the bedGraph file should
  43 be continuous so only bedGraph files from MACS3 are acceptable.
  44
  45 ## Command Line Options
  46
  47 The command line options for `bdgdiff` are defined in `/MACS3/Commands/bdgdiff_cmd.py` and `/bin/macs3` files. Here is a brief overview of these options:
  48
  49 - `--t1`: MACS pileup bedGraph for condition 1. Incompatible with
  50   callpeak --SPMR output. REQUIRED
  51 - `--t2`: MACS pileup bedGraph for condition 2. Incompatible with
  52   callpeak --SPMR output. REQUIRED
  53 - `--c1`: MACS control lambda bedGraph for condition 1. Incompatible
  54   with callpeak --SPMR output. REQUIRED
  55 - `--c2`: MACS control lambda bedGraph for condition 2. Incompatible
  56   with callpeak --SPMR output. REQUIRED
  57 - `-C` or `--cutoff`: log10LR cutoff. Regions with signals lower than
  58   the cutoff will not be considered as enriched regions. DEFAULT: 3
  59   (likelihood ratio=1000)
  60 - `-l` or `--min-len`: Minimum length of the differential region. Try
  61   a bigger value to remove small regions. DEFAULT: 200
  62 - `-g` or `--max-gap`: Maximum gap to merge nearby differential
  63   regions. Consider a wider gap for broad marks. The maximum gap
  64   should be smaller than the minimum length (-g). DEFAULT: 100
  65 - `--d1` or `--depth1`: Sequencing depth (# of non-redundant reads in
  66   million) for condition 1. It will be used together with --d2. See
  67   the description for --d2 below for how to assign them. Default: 1
  68 - `--d2` or `--depth2`: Sequencing depth (# of non-redundant reads in
  69   million) for condition 2. It will be used together with --d1. DEPTH1
  70   and DEPTH2 will be used to calculate the scaling factor for each
  71   sample, to down-scale the larger sample to the level of the smaller
  72   one. For example, while comparing 10 million condition 1 and 20
  73   million condition 2, use --d1 10 --d2 20, then the pileup value in
  74   bedGraph for condition 2 will be divided by 2. Default: 1
  75 - `--verbose`: Set the verbose level of runtime messages. 0: only show
  76   critical messages, 1: show additional warning messages, 2: show
  77   process information, 3: show debug messages. DEFAULT: 2
  78 - `--outdir`: If specified, all output files will be written to that
  79   directory. Default: the current working directory
  80 - `--o-prefix`: Output file prefix. Actual files will be named as
  81   PREFIX_cond1.bed, PREFIX_cond2.bed, and PREFIX_common.bed. Mutually
  82   exclusive with -o/--ofile.
  83 - `-o` or `--ofile`: Output filenames. Must give three arguments in
  84   order: 1. file for unique regions in condition 1; 2. file for unique
  85   regions in condition 2; 3. file for common regions in both
  86   conditions. Note: mutually exclusive with --o-prefix.
  87
  88
  89 ## Example Usage
  90
  91 Here is an example of how to use the `bdgdiff` command:
  92
  93 ```bash
  94 macs3 bdgdiff -t1 treatment1.bedGraph -c1 control1.bedGraph -t2 treatment2.bedGraph -c2 control2.bedGraph --depth1 1.0 --depth2 1.0 -o output.bedGraph --minlen 500 --maxgap 1000 --cutoff 1.0
  95 ```
  96
  97 In this example, the program will call differential peaks from the two
  98 pairs of treatment and control bedGraph files and write the result to
  99 `output.bedGraph`. The depth of the first and second condition is set
 100 to 1.0, the minimum length of differential peaks is set to 500, the
 101 maximum gap between differential peaks is set to 1000, and the cutoff
 102 for log10LR to call differential peaks is set to 1.0 (or likelihood
 103 ratio 10).
 104
 105 ## Step-by-step Instruction for calling differential peaks
 106
 107 In this chatper, we will describe how to use MACS3 to identify
 108 differential regions by comparing pileup tracks of two conditions,
 109 starting from the alignment files. Two modules will be involved:
 110 `callpeak` and `bdgdiff` ( `predictd` is optional ). We will use human
 111 ChIP-seq data as example, and filenames of raw data are:
 112 cond1_ChIP.bam and cond1_Control.bam for condition 1; cond2_ChIP.bam
 113 and cond2_Control.bam for condition 2.
 114
 115 ### Step 1: Generate pileup tracks using callpeak module
 116
 117 Purpose of this step is to use `callpeak` with -B option to generate
 118 bedGraph files for both conditions. There are several things to
 119 remember: 1. `--SPMR` is not compatible with `bdgdiff`, so avoid using
 120 it; 2. prepare a pen to write down the number of non-redundant reads
 121 of both conditions -- you will find such information in runtime
 122 message or xls output from `callpeak`; 3. keep using the same
 123 `--extsize` for both conditions ( you can get it from `predictd`
 124 module).
 125
 126 To get a uniform extension size for running `callpeak`, run `predictd`:
 127
 128 ```
 129  $ macs3 predictd -i cond1_ChIP.bam
 130
 131  $ macs3 predictd -i cond2_ChIP.bam
 132 ```
 133
 134 An easy solution is to use the average of two 'fragment size'
 135 predicted in `callpeak`, however any reasonable value will work. For
 136 example, you can use `200` for most ChIP-seq datasets for
 137 transcription factors, or ''147'' for most histone modification
 138 ChIP-seq. The only requirement is that you have to keep using the same
 139 extsize for the following commands:
 140
 141 ```
 142  $ macs3 callpeak -B -t cond1_ChIP.bam -c cond1_Control.bam -n cond1 --nomodel --extsize 120
 143
 144  $ macs3 callpeak -B -t cond2_ChIP.bam -c cond2_Control.bam -n cond2 --nomodel --extsize 120
 145 ```
 146
 147 Pay attention to runtime message, or extract the "tags after filtering in treatment" and "tags after filtering in control" lines from xls to see the effective sequencing depths for both conditions. In our previous command lines, '--to-large' is not used, so the effective sequencing depth is the smaller number of treatment and control. For example:
 148
 149 ```
 150  $ egrep "tags after filtering in treatment|tags after filtering in control" cond1_peaks.xls
 151  # tags after filtering in treatment: 19291269
 152  # tags after filtering in control: 12914669
 153
 154  $ egrep "tags after filtering in treatment|tags after filtering in control" cond2_peaks.xls
 155  # tags after filtering in treatment: 19962431
 156  # tags after filtering in control: 14444786
 157 ```
 158
 159 Then actual effective depths of condition 1 and 2 are: 12914669
 160 and 14444786. Keep record of these two numbers and we will use them
 161 later. After successfully running '''callpeak''', you will have
 162 ''cond1_treat_pileup.bdg'', ''cond1_control_lambda.bdg'',
 163 ''cond2_treat_pileup.bdg'', and ''cond2_control_lambda.bdg'' in the
 164 working directory.
 165
 166 ### Step 2: Call differential regions
 167
 168 The purpose of this step is to do a three ways comparisons to find out
 169 where in the genome has differential enrichment between two
 170 conditions. A basic requirement is that this region should be at least
 171 enriched in either condition. A log10 likelihood ratio cutoff (C) will
 172 be applied in this step. Three types of differential regions will be
 173 reported: 1. those having more enrichment in condition 1 over
 174 condition 2 ( cond1_ChIP > cond1_Control and cond1_ChIP > cond2_ChIP
 175 ); 2. those having more enrichment in condition 2 over condition 1 (
 176 cond2_ChIP > cond2_Control and cond2_ChIP > cond1_ChIP ); those having
 177 similar enrichment in both conditions ( cond1_ChIP > cond1_Control and
 178 cond2_ChIP > cond2_Control and cond1_ChIP ≈ cond1_ChIP ).
 179
 180 Run this:
 181
 182 ```
 183  $ macs3 bdgdiff --t1 cond1_treat_pileup.bdg --c1 cond1_control_lambda.bdg --t2 cond2_treat_pileup.bdg\
 184    --c2 cond2_control_lambda.bdg --d1 12914669 --d2 14444786 -g 60 -l 120 --o-prefix diff_c1_vs_c2
 185 ```
 186
 187 You will get the following three files in working directory:
 188
 189  1. `diff_c1_vs_c2_c3.0_cond1.bed` This file stores regions that are
 190  highly enriched in condition 1 comparing to condition 2. The last
 191  column in the file represent the log10 likelihood ratio to show how
 192  likely the observed signal in condition 1 in this region is from
 193  condition 1 comparing to condition 2. The higher the value, the
 194  greater the difference.
 195
 196  2. `diff_c1_vs_c2_c3.0_cond2.bed` This file stores regions that are
 197  highly enriched in condition 2 comparing to condition 1. The last
 198  column in the file represent the log10 likelihood ratio to show how
 199  likely the observed signal in condition 2 in this region is from
 200  condition 2 comparing to condition 1. Higher the value, bigger the
 201  difference.
 202
 203  3. `diff_c1_vs_c2_c3.0_common.bed` This file stores regions that are
 204  highly enriched in both condition 1 and condition 2, and the
 205  difference between condition 1 and condition 2 is not
 206  significant. The last column in the file represent the difference
 207  between condition 1 and condition 2 in log10 likelihood ratios.