Q: HMMRATAC producing too many peaks

[igordot]

I ran HMMRATAC on two samples that should be biological replicates. The quality difference between them is high, although the quality is far from ideal for both. However, the low quality sample produced more than 10x as many peaks (<10k versus >100k) with the default settings.

For comparison, the Java version of HMMRATAC produced <1k peaks for the first sample and no peaks for the second sample (score cutoff of 30).

Here is a genome browser screenshot to illustrate the situation:

The called peaks for the first sample look like peaks. Not so much for the second sample. Is this expected? Are there some guidelines to help avoid this?

[taoliu]

One thing you can try is to increase the cutoff for calling candidate regions in the macs3 version of hmmratac using the option for prescan-cutoff -c. And perhaps first run the --cutoff-analysis in hmmratac to decide the OK parameter for -c. Also it is likely that the cutoff values for -u and -l for training HMM have to be adjusted as well. Check https://macs3-project.github.io/MACS/docs/hmmratac.html#choices-of-cutoff-values

Please let me know if this can help.

[igordot]

I saw that, but wasn't sure how to properly use the output to determine the optimal values. Do you have an example documented somewhere?

[taoliu]

The cutoff-analysis for hmmratac command is discussed here: https://macs3-project.github.io/MACS/docs/hmmratac.html#choices-of-cutoff-values

And a similar cutoff-analysis for callpeak command is discussed here: https://macs3-project.github.io/MACS/docs/callpeak.html#cutoff-analysis

And if you can generate the cutoff analysis report from hmmratac --cutoff-analysis-only and share with me here, we can analyze the report together in this thread.

[igordot]

I am sharing the cutoff report here (I had to add .txt because curiously GitHub does not allow .tsv):
test_cutoff_analysis.tsv.txt

It might even be easier to just paste it in whole:

score	npeaks	lpeaks	avelpeak
638.79	1	234	234.00
630.39	1	236	236.00
621.98	1	267	267.00
613.58	1	267	267.00
605.17	1	277	277.00
596.77	1	280	280.00
588.36	1	291	291.00
579.96	1	292	292.00
571.55	3	794	264.67
563.15	2	1972	986.00
554.74	2	2523	1261.50
546.34	3	2902	967.33
537.93	3	2913	971.00
529.53	3	2920	973.33
521.12	3	2925	975.00
512.71	3	2955	985.00
504.31	3	2989	996.33
495.90	3	3183	1061.00
487.50	3	3202	1067.33
479.09	3	3257	1085.67
470.69	3	3407	1135.67
462.28	3	3458	1152.67
453.88	3	3501	1167.00
445.47	3	3505	1168.33
437.07	3	3510	1170.00
428.66	3	3563	1187.67
420.26	3	3639	1213.00
411.85	3	3678	1226.00
403.45	3	3735	1245.00
395.04	3	4062	1354.00
386.64	3	4456	1485.33
378.23	3	4507	1502.33
369.83	3	4511	1503.67
361.42	3	4539	1513.00
353.02	3	4568	1522.67
344.61	3	4577	1525.67
336.21	3	4677	1559.00
327.80	3	4704	1568.00
319.40	3	5206	1735.33
310.99	2	6264	3132.00
302.58	3	6392	2130.67
294.18	3	6401	2133.67
285.77	3	6412	2137.33
277.37	3	6426	2142.00
268.96	3	6467	2155.67
260.56	3	6484	2161.33
252.15	3	6502	2167.33
243.75	3	6508	2169.33
235.34	4	6635	1658.75
226.94	5	6762	1352.40
218.53	5	6785	1357.00
210.13	6	6967	1161.17
201.72	6	6986	1164.33
193.32	6	7039	1173.17
184.91	6	7067	1177.83
176.51	7	7301	1043.00
168.10	7	7333	1047.57
159.70	8	7900	987.50
151.29	8	8338	1042.25
142.89	10	9056	905.60
134.48	13	9942	764.77
126.08	14	11087	791.93
117.67	18	12245	680.28
109.27	23	18535	805.87
100.86	27	25300	937.04
92.46	31	34402	1109.74
84.05	30	45057	1501.90
75.65	31	52352	1688.77
67.24	30	63203	2106.77
58.84	29	73885	2547.76
50.43	26	82657	3179.12
42.03	27	88598	3281.41
33.62	28	92070	3288.21
25.22	46	101435	2205.11
16.81	414	255592	617.37
8.40	14319	7995323	558.37
0.00	64045	2385081806	37240.72

[taoliu]

The fold change values for reasonable peaks are within the lower range. The default parameters for hmmratac are -u 20 -l 10 -c 1.2. It's possible that we need to increase the -c parameter a bit higher. I can't judge since the current npeak, lpeak, and avelpeak values of the score cutoff range 0.00-8.40 show huge difference and we don't have a finer resolution of the analysis between 0 and 8.4. You can try to use -u 20 -l 8 -c 5, and to see if the result can improve. We will try to add an option for fine-tuning this cutoff analysis to get better resolution in the future. At the meantime, we are adding a different emission model -- Poisson emission for HMM, as an addition to the current Gaussian emission model. In our internal test on some small dataset, it gives us more similar results as the java version HMMRATAC.

How is the current HMM from hmmratac? Could you share us the command line output as https://macs3-project.github.io/MACS/docs/hmmratac.html#tune-the-hmm-model?

[Jungal10]

I ran HMMRATAC on two samples that should be biological replicates. The quality difference between them is high, although the quality is far from ideal for both. However, the low quality sample produced more than 10x as many peaks (<10k versus >100k) with the default settings.

For comparison, the Java version of HMMRATAC produced <1k peaks for the first sample and no peaks for the second sample (score cutoff of 30).

Here is a genome browser screenshot to illustrate the situation:

The called peaks for the first sample look like peaks. Not so much for the second sample. Is this expected? Are there some guidelines to help avoid this?

Not a direct answer, but the scale on your IGV is quite different as well between the two

[igordot]

Not a direct answer, but the scale on your IGV is quite different as well between the two

Yes, each sample is auto-scaled independently since peak-calling is independent as well.

[taoliu]

Also related to #638, it's highly possible that the hidden markov model learned from the low quality data can't capture the expected chromatin structure around accessible sites. This is a common problem of "machine learning". In this case, perhaps a general peak calling process such as using "callpeak" function is more appropriate.

[igordot]

I am not surprised that machine learning is failing for these samples. I was primarily concerned about the large difference between the Java and MACS3 implementations.

[taoliu]

We are implementing a different emission model for HMM in MACS3 called the Poisson Emission model. During our internal tests, this simple model (compared with the Gaussian Emission model) generated results that were similar to those of the Java implementation. Stay tuned if you are interested in it :)