methylpy
methylpy copied to clipboard
yupenghe
binom-test failed
Hi yupeng, thanks for this powerful package. I ran methylpy 1.4.6 on a plant genome. Everything went smoothly except that binom-test failed.The error was like this:
Perform binomial test Fri Sep 9 03:23:16 2022
Traceback (most recent call last):
File "/home/many/software/methylpy-1.2.9/bin/methylpy", line 5, in
However, another sample was done successfully, and the two samples shared the same pipeline and reference genome. I checked their allc files, and there is no difference in fields[4]. They look like this:
test.zip
So can you figure out what is wrong with my sample? Appreciate for your help.
Hey, what command did you run for getting the error? Also, doing binomial test requires unmethylated control such as unmethylated lambda sequence. Please make sure that you have such control.
Got it. Do you have chrC in the reference genome? Are reads mapped to this chromosome?
binomial test is not required for DMRs and most analyses. So if you don't need it, you can do --binom-test False to get around this.
yes,the chloroplast genome is included in my reference genome. And the idxstats of bam file looks fine.
Hmm weird. Can you grep the chrC lines in allc files and share it with me for the problematic sample?
hi yupeng, I found that the last line of the problemmatic allc file is shifted, leaving a blank in fields[4]
And when I delete the last line and run methylpy test-allc, binom-test work.
So how can I prevent this from happening? 3 out of 4 samples had such problem. I attached the screen shots of their allc files.
Thank you so much.
Can you send me the last line in a file? If possible, can you also send me the reads (in bam) that overlap the last CpG of chrC? I believe you can do something like
samtools view -b BAM_FILE chrC:156341-156345 > TEST.bam
Hi, the test files are enclosed in the zip file. And for problematic samples, 2 BAM_FILEs were generated by methylpy paired-end-pipeline:
no_clonal.bam and no_clonal.bam.read2flipped.bam.
I did what you told me for both of the bamfiles.
And I am so sorry that I had reckressly discarded the last line of my incorrect allc file. So I only put the right one in the file.
Thanks a lot!
One more thing needed to test this. Do you mind to send me the chrC sequence you are using?
In case this is helpful, here is a command to get sequence of just one chromosome.
samtools faidx GENOME.fa chrC > chrC.fa
That is weird. I was not able to reproduce the error. If you do --unmethylated-control C: instead of chrC, will it help?
methylpy paired-end-pipeline
--read1-files cleaned_DH-1_R1.fastq.gz
--read2-files cleaned_DH-1_R2.fastq.gz
--sample cleaned_DH-1
--forward-ref AcDH_f
--reverse-ref AcDH_r
--ref-fasta clean.AcDH.fasta
--num-procs 52
--remove-clonal True
--path-to-output DH-1_output/
--sort-mem 50G
--unmethylated-control C:
--binom-test True
--path-to-picard /home/many/software