hmmratac
Description
HMMRATAC (macs3 hmmratac) is a dedicated peak calling algorithm
based on Hidden Markov Model for ATAC-seq data. The basic idea behind
HMMRATAC is to digest ATAC-seq data according to the fragment length
of read pairs into four signal tracks: short fragments,
mono-nucleosomal fragments, di-nucleosomal fragments and
tri-nucleosomal fragments. Then integrate the four tracks using Hidden
Markov Model to consider three hidden states: open region, nucleosomal
region, and background region. The orginal
paper was
published in 2019, and the original software was written in JAVA, by
the then PhD student Evan Tarbell, a mohawk bioinformatician. In MACS3
project, we implemented HMMRATAC idea in Python/Cython and optimize
the whole process using existing MACS functions and hmmlearn.
Here’s an example of how to run the hmmratac command:
$ macs3 hmmratac -i yeast.bam -n yeast
or with the BEDPE format
$ macs3 hmmratac -i yeast.bedpe.gz -f BEDPE -n yeast
Note: you can convert BAMPE to BEDPE by using
$ macs3 filterdup --keep-dup all -f BAMPE -i yeast.bam -o yeast.bedpe
Please use macs3 hmmratac --help to see all the options. Here we
list the essential ones.
Essential Options
-i INPUT_FILE [INPUT_FILE ...] / --input INPUT_FILE [INPUT_FILE ...]
This is the only REQUIRED parameter for hmmratac. Input files
containing the aligment results for ATAC-seq paired end reads. If
multiple files are given as ‘-t A B C’, then they will all be read and
pooled together. The file should be in BAMPE or BEDPE format (aligned
in paired end mode). Files can be gzipped. Note: all files should be
in the same format. REQUIRED.
-f {BAMPE,BEDPE} / --format {BAMPE,BEDPE}
Format of input files, “BAMPE” or “BEDPE”. If there are multiple
files, they should be in the same format – either BAMPE or
BEDPE. Please note that the BEDPE only contains three columns –
chromosome, left position of the whole pair, right position of the
whole pair– and is NOT the same BEDPE format used by BEDTOOLS. To
convert BAMPE to BEDPE, you can use this command macs3 filterdup --keep-dup all -f BAMPE -i input.bam -o output.bedpe. DEFAULT:
“BAMPE”.
--outdir OUTDIR
If specified all output files will be written to that directory. Default: the current working directory
-n NAME/ --name NAME
Name for this experiment, which will be used as a prefix to generate output file names. DEFAULT: “NA”
--modelonly
This option will only generate the HMM model as a JSON file and
quit. This model can then be applied using the --model
option. Default: False
--model
If provided, HMM training will be skipped and a JSON file generated from a previous HMMRATAC run will be used instead of creating new one. Default: NA
-t HMM_TRAINING_REGIONS / --training HMM_TRAINING_REGIONS
Customized training regions can be provided through this option. -t
takes the filename of training regions (previously was BED_file) to
use for training HMM, instead of using foldchange settings to
select. Default: NA
--min-frag-p MIN_FRAG_P
We will exclude the abnormal fragments that can’t be assigned to any of the four signal tracks. After we use EM to find the means and stddevs of the four distributions, we will calculate the likelihood that a given fragment length fit any of the four using normal distribution. The criteria we will use is that if a fragment length has less than MIN_FRAG_P probability to be like either of short, mono, di, or tri-nuc fragment, we will exclude it while generating the four signal tracks for later HMM training and prediction. The value should be between 0 and 1. Larger the value, more abnormal fragments will be allowed. So if you want to include more ‘ideal’ fragments, make this value smaller. Default = 0.001
--cutoff-analysis-only
Only run the cutoff analysis and output a report. After generating
the report, the process will stop. The report will help user decide
the three crucial parameters for -l, -u, and -c. So it’s highly
recommanded to run this first! Please read the report and
instructions in Choices of cutoff values
on how to decide the three crucial parameters.
-u HMM_UPPER / --upper HMM_UPPER
Upper limit on fold change range for choosing training sites. This is
an important parameter for training so please read. The purpose of
this parameter is to EXCLUDE those unusually highly enriched chromatin
regions so we can get training samples in ‘ordinary’ regions
instead. It’s highly recommended to run the --cutoff-analysis-only
first to decide the lower cutoff -l, the upper cutoff -u, and the
pre-scanning cutoff -c. The upper cutoff should be the cutoff in the
cutoff analysis result that can capture some (typically hundreds of)
extremely high enrichment and unusually wide peaks. Default: 20
-l HMM_LOWER / --lower HMM_LOWER
Lower limit on fold change range for choosing training sites. This is
an important parameter for training so please read. The purpose of
this parameter is to ONLY INCLUDE those chromatin regions having
ordinary enrichment so we can get training samples to learn the common
features through HMM. It’s highly recommended to run the
--cutoff-analysis-only first to decide the lower cutoff -l, the
upper cutoff -u, and the pre-scanning cutoff -c. The lower cutoff
should be the cutoff in the cutoff analysis result that can capture
moderate number ( about 10k ) of peaks with normal width ( average
length 500-1000bps long). Default: 10
-c PRESCAN_CUTOFF / --prescan-cutoff PRESCAN_CUTOFF
The fold change cutoff for prescanning candidate regions in the whole
dataset. Then we will use HMM to predict/decode states on these
candidate regions. The higher the prescan cutoff, the fewer regions
will be considered. Must be > 1. This is an important parameter for
decoding so please read. The purpose of this parameter is to EXCLUDE
those chromatin regions having noises/random enrichment so we can have
a large number of possible regions to predict the HMM states. It’s
highly recommended to run the --cutoff-analysis-only first to decide
the lower cutoff -l, the upper cutoff -u, and the pre-scanning
cutoff -c. The pre-scanning cutoff should be the cutoff close to the
BOTTOM of the cutoff analysis result that can capture a large number
of possible peaks with normal length (average length 500-1000bps). In
most cases, please do not pick a cutoff too low that captures almost
all the background noises from the data. Default: 1.2
Choices of cutoff values
Before you proceed, it’s highly recommended to run with
--cutoff-analysis-only for the initial attempt. When this option is
activated, hmmratac will use EM to estimate the best parameters for
fragment sizes of short fragments, mono-, di-, and tri-nucleosomes,
pileup fragments, convert the pileup values into fold-change, and
analyze each possible cutoff. This analysis includes the number of
peaks that can be called, their average peak length, and their total
length. After the report is generated, you can review its contents and
decide on the optimal -l, -u, and -c.
The report consists of four columns:
Score: the possible fold change cutoff value.
npeaks: the number of peaks.
lpeaks: the total length of all peaks.
avelpeak: the average length of peaks.
While there’s no universal rule, here are a few suggestions:
The lower cutoff should be the cutoff in the report that captures a moderate number (about 10k) of peaks with a normal width (average length 500-1000bps long).
The upper cutoff should capture some (typically hundreds of) extremely high enrichment and unusually wide peaks in the report. The aim here is to exclude abnormal enrichment caused by artifacts such as repetitive regions.
The pre-scanning cutoff should be the cutoff close to the BOTTOM of the report that can capture a large number of potential peaks with a normal length (average length 500-1000bps). However, it’s recommended not to use the lowest cutoff value in the report as this may include too much noise from the genome.
Tune the HMM model
It’s highly recommended to check the runtime message of the HMM model after training. An example is like this:
#4 Train Hidden Markov Model with Multivariate Gaussian Emission
# Extract signals in training regions with bin size of 10
# Use Baum-Welch algorithm to train the HMM
# HMM converged: True
# Write HMM parameters into JSON: test_model.json
# The Hidden Markov Model for signals of binsize of 10 basepairs:
# open state index: state2
# nucleosomal state index: state1
# background state index: state0
# Starting probabilities of states:
# bg nuc open
# 0.7994 0.1312 0.06942
# HMM Transition probabilities:
# bg nuc open
# bg-> 0.9842 0.01202 0.003759
# nuc-> 0.03093 0.9562 0.01287
# open-> 0.007891 0.01038 0.9817
# HMM Emissions (mean):
# short mono di tri
# bg: 0.2551 1.526 0.4646 0.07071
# nuc: 6.538 17.94 3.422 0.05819
# open: 5.016 17.47 6.897 2.121
We will ‘guess’ which hidden state is for the open region, which is the nucleosomal region, and which is the background. We compute from the HMM Emission matrix to pick the state with the highest sum of mean signals as the open state, the lowest as the backgound state, then the rest is the nucleosomal state. However it may not work in every case. In the above example, it may be tricky to call the second row as ‘nuc’ and the third as ‘open’. If the users want to exchange the state assignments of the ‘nuc’ and ‘open’, they can modify the state assignment in the HMM model file (e.g. test_model.json). For the above example, the model.json looks like this (we skipped some detail):
{"startprob": [...], "transmat": [...], "means": [...], "covars": [...],
"covariance_type": "full", "n_features": 4,
"i_open_region": 2, "i_background_region": 0, "i_nucleosomal_region": 1,
"hmm_binsize": 10}
We can modify the assignment of: "i_open_region": 2, "i_background_region": 0, "i_nucleosomal_region": 1, by assigning 1
to open, and 2 to nucleosomal as: "i_open_region": 1, "i_background_region": 0, "i_nucleosomal_region": 2, Then save the
HMM in a new model file such as new_model.json.
Then next, we can re-run macs3 hmmratac with the same parameters
plus an extra option for the HMM model file like macs3 hmmratac --model new_model.json