Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Q: How to use HMMRATAC to give us a narrow peak file and have signalValue, p and q-value ? #648

Open
MinaMyr opened this issue Jun 13, 2024 · 1 comment

Comments

@MinaMyr
Copy link

MinaMyr commented Jun 13, 2024

Use case
Hello I want to do peak calling on my ATAC seq data, and for that I use, callpeak and hmmratac.

Describe the problem
My problem here is with hmmratac;

I use this command macs3 hmmratac -f BAM -i RES/atac/try_1/04-NobiasedRegions/ATAC20_D0_rep1_hg38_sort_dedup_biasedRegions.bam -u 15 -l 5 -c 100 -n RES/atac/try_1/06-PeakCalling/hmmratac_narrow/ATAC20_D0_rep1
I read in the doc that we should have a narrow peak file format but I have an *_accessible_regions.gappedPeak format which is not what I want. For the rest of my analysis I need a narrow peak format.

Describe the solution you tried
So I tried to create my narrow peak file from the gappedPeak, by hand, but the trouble I am now having is that for the signal value, p-value and q-value are equals to 0 (I think like this comment said it is intentional #592 (comment) ). Here is what my file look like;

track name="peak" description="peak" type=gappedPeak nextItemButton=on
20	208210	220320	peak_1	0	.	208450	220210	0	3	1240,8820,1550	240,1520,10450	0	0	0
20	220410	267710	peak_2	0	.	220490	267640	0	1	47150	80	0	0	0
20	279100	302420	peak_3	0	.	279110	302390	0	8	890,9080,6700,310,240,3640,1240,660	10,1010,10140,17070,17420,17720,21380,22630	0	0	0
20	302430	314800	peak_4	0	.	302470	314780	0	4	8730,450,2840,130	40,8880,9350,12220	0	0	0
20	314890	342820	peak_5	0	.	314920	342760	0	7	7230,4790,1170,3260,2190,4780,4170	30,7400,12210,13400,16670,18880,23700	0	0	0
20	342840	365480	peak_6	0	.	342880	365440	0	6	5620,10270,120,1640,2240,2450	40,5690,15990,16150,17830,20150	0	0	0
20	368100	402600	peak_7	0	.	368140	402540	0	4	2600,3260,230,28060	40,2650,5950,6380	0	0	0
20	443320	447340	peak_8	0	.	443330	447170	0	2	2860,910	10,2940	0	0	0
20	447420	479680	peak_9	0	.	447440	479630	0	2	15710,16460	20,15750	0	0	0
20	479740	480080	peak_10	0	.	479850	480010	0	1	160	110	0	0	0

I don't really know what to do now, since signal Value p-value and q-value are important for the rest of my analysis.
Thank you in advance !

@nickgladman
Copy link

nickgladman commented Jun 15, 2024

Hello and thanks for making this resource. I am commenting here to track this since I'm having the same issue: no narrowPeak output from the hmmratac using mapped .bam files with the following code:
macs3 hmmratac \ -i D*rep1_paired.sorted.bam \ *rep2_paired.sorted.bam \ --outdir directory \ -n output_name \ -u 20 \ -l 10 \ -c 2

I previously ran the with the --cutoff-analysis-only flag prior to entering the -u -l -c params. I'm using version 3.0.1 in a conda environment with a SGE scheduler. There are no errors in the output. Here are the last lines of the run output:

INFO @ 14 Jun 2024 17:06:54: [7795 MB] # Write the output... INFO @ 14 Jun 2024 17:09:30: [7795 MB] # Write accessible regions in a gappedPeak file: *accessible_regions.gappedPeak

Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment
Projects
None yet
Development

No branches or pull requests

2 participants