methylpy add-methylation-level issue about methylpy HOT 15 CLOSED

It is weird. Allc file should be sorted. I can look into this if you can share some files to reproduce the issue. Also, it would be helpful to check whether the bam file created is sorted.

from methylpy.

Thanks for the response.
Looks like the "*processed_reads_no_clonal.bam " file is properly sorted.
What kind of data can I provide you?

from methylpy.

It will be great if you can provide the processed_reads_no_clonal.bam file that was used to generate the unsorted allc files.

from methylpy.

How about this?

subSample.bam.zip

from methylpy.

Thanks. Were you able to produce unsorted allc file from this bam file? Can you also point me to the genome (FASTA)?

from methylpy.

Were you able to produce unsorted allc file from this bam file?
Yes. I simply used all the example code on the webpage.
The genome is downloaded from here
ftp://ftp.ensemblgenomes.org/pub/plants/release-46/fasta/solanum_lycopersicum/dna/Solanum_lycopersicum.SL3.0.dna.toplevel.fa.gz

from methylpy.

I don't think so. The allc file I got from the bam file is sorted. Here is the command I used

methylpy call-methylation-state --input-file subSample.bam --sample test --ref-fasta Solanum_lycopersicum.SL3.0.dna.toplevel.fa --paired-end True

from methylpy.

That is good to know. Here is what I used. Is there anything wrong that I did not realize with this?
I looped each sample for paried end pipeline.

methylpy paired-end-pipeline
--read1-files ./_1.fq.gz --read2-files ./_2.fq.gz --sample ${folder[$indi]}
--merge-by-max-mapq True
--binom-test True
--unmethylated-control chloroplast
--min-cov 3
--forward-ref $fref --reverse-ref $rref --ref-fasta $reff
--num-procs 20 --sort-mem 35000000000
--path-to-output $directout
--remove-clonal True
--path-to-picard="/softwares/picard/"
--aligner-options "-p 8"
--trim-read False
1>$stdoutfile 2>$errfile

And used this function add-methylation-level

methylpy add-methylation-level
--input-tsv-file methylpy_wgenome.tsv
--output-file $outfile
--allc-files allc_ACA.tsv.gz allc_ACB.tsv.gz
--samples ACA-st ACB-st
--mc-type $mctype
--num-procs 4
1>$stdoutfile1 2>$errfile1

Then I do can call DMR on the allc file properly.

methylpy DMRfind
--allc-files allc_ACA.tsv.gz allc_ACB.tsv.gz allc_pCMT3-RNAiA.tsv.gz allc_pCMT3-RNAiB.tsv.gz
--samples ACA ACB pCMT3-RNAiA pCMT3-RNAiB
--mc-type $mctype
--sample-category wild wild RNAi RNAi
--min-cluster 2
--chroms 1 2 3 4 5 6 7 8 9 10 11 12
--num-procs 100
--dmr-max-dist 1000
--sig-cutoff 0.05
--output-prefix $prefixname1
1>$stdoutfile1 2>$errfile1

from methylpy.

The commands look good to me. Did you get unsorted allc files and end up with some error in the DMR finding step?

from methylpy.

Right. I got unsorted allc files. I could do DMR finding. But I got errors when I try to add methylation levels using the allc files. In the end, I just unzipped the allc file, sorted them, and then feed the files back to add methylation level, it worked. Anyway, as long as I did not get any errors or warnings with the DMR finding. That should be fine, correct?

from methylpy.

That is weird. I would recommend rerunning DMR finding on the sorted allc files to make sure things are working fine.

from methylpy.

I have a follow-up question- what does it mean if my output has NA?

from methylpy.

It means that there is no reads to estimate methylation level in that region in that sample.

from methylpy.

Hi
I went back and checked the output of this and it says this

samtools sort: couldn't allocate memory for bam_mem
[bam_sort_core] merging from 79 files and 1 in-memory blocks...
INFO 2020-09-21 20:19:00 MarkDuplicates

I guess my bam and allc file not sorted has to do with my settings of the number of processors and the sort memory being nonproportional to each other?
I am attaching the whole output here. Should I consider rerun?
https://www.dropbox.com/s/2vbkmiexmqqdzzc/ACA.log?dl=0

from methylpy.

Ah, that is good to know. You may want to drop --sort-mem 35000000000 option or replace it with something like --sort-mem 1G.

You don't have to rerun from scratch. Manually doing samtools sort to get sorted bam file and then regenerating allc files using methylpy call-methylation-state will do.

from methylpy.