epi2me-labs / pore-c-py

Other
12 stars 4 forks source link

An unknown error occurred #9

Open GLking123 opened 5 months ago

GLking123 commented 5 months ago

Dear author: Thank you very much for developing such a great software. For various reasons, I am unable to use the automated process of wf-pore-c. I just want to align pore-c data with the reference genome, such as first locating the restriction sites, and then using minimap2 for alignment, but I followed the method below and it threw an error?

`################################# workdir=/work/home/gonglei/data/new_new_minimap refdir=/work/home/gonglei/data/new_data/A.fasta fqdir=/work/home/gonglei/data/11.8porec1/pass.fq.gz ENZYME="NlaIII" OUTPUT="all"

source /work/home/gonglei/miniconda3/etc/profile.d/conda.sh conda activate samtools export PATH="/work/home/gonglei/miniconda3/envs/samtools/bin:$PATH" SOFTWARE_PATH="/work/home/gonglei/miniconda3/envs/samtools/bin/"

cd $workdir pore-c-py digest "${fqdir}" "${ENZYME}" \ | samtools fastq -T '*' \ | /work/home/gonglei/minimap2/minimap2 -ay -t 50 -x map-ont --split-prefix tmp "${refdir}" - \ | pore-c-py annotate - "${OUTPUT}" --monomers --stdout --summary --chromunity \ | tee "${OUTPUT}.ns.bam" \ | samtools sort --write-index -o "${OUTPUT}.cs.bam" - samtools index "${OUTPUT}.ns.bam"`

The log file is as follows: [19:00:47 - Digest ] Digesting concatemers from [PosixPath('/work/home/gonglei/data/11.8porec1/pass.fq.gz')]. [19:00:47 - AnntateBAM] Processing reads from - [19:00:47 - Digest ] Digesting unaligned sequences from <pysam.libcalignmentfile.AlignmentFile object at 0x2b912f639360> [M::mm_idx_gen::120.9921.51] collected minimizers [M::mm_idx_gen::140.2013.92] sorted minimizers [M::main::140.2033.92] loaded/built the index for 79 target sequence(s) [M::mm_mapopt_update::143.0413.86] mid_occ = 3923 [M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 79 [M::mm_idx_stat::144.8833.82] distinct minimizers: 110678714 (34.28% are singletons); average occurrences: 13.815; average spacing: 5.272; total length: 8061357846 [M::worker_pipeline::530.24419.09] mapped 1513164 sequences [M::worker_pipeline::23531.60846.98] mapped 1511410 sequences [01:33:40 - Digest ] Found 190532421 monomers in 15726183 concatemers. [01:33:40 - Digest ] Excluded 0 concatemers. [01:33:40 - Digest ] Finished digestion. [M::bam2fq_mainloop] discarded 0 singletons [M::bam2fq_mainloop] processed 190532421 reads [M::worker_pipeline::23710.69947.00] mapped 1518571 sequences [M::worker_pipeline::23783.78547.00] mapped 645849 sequences [M::mm_idx_gen::23916.58546.75] collected minimizers [M::mm_idx_gen::23926.30546.74] sorted minimizers [M::main::23926.30546.74] loaded/built the index for 69 target sequence(s) [M::mm_mapopt_update::23926.30546.74] mid_occ = 3923 [M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 69 [M::mm_idx_stat::23929.07346.73] distinct minimizers: 111733303 (34.23% are singletons); average occurrences: 14.019; average spacing: 5.273; total length: 8259067193 [M::mm_idx_gen::24057.93246.49] collected minimizers [M::mm_idx_gen::24069.93546.48] sorted minimizers [M::main::24069.93546.48] loaded/built the index for 88 target sequence(s) [M::mm_mapopt_update::24069.93546.48] mid_occ = 3923 [M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 88 [M::mm_idx_stat::24072.03446.47] distinct minimizers: 109953755 (34.14% are singletons); average occurrences: 13.838; average spacing: 5.273; total length: 8022362377 [M::mm_idx_gen::24202.89146.23] collected minimizers [M::mm_idx_gen::24215.63746.22] sorted minimizers [M::main::24215.63746.22] loaded/built the index for 96 target sequence(s) [M::mm_mapopt_update::24215.63746.22] mid_occ = 3923 [M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 96 [M::mm_idx_stat::24217.39446.21] distinct minimizers: 110845302 (33.84% are singletons); average occurrences: 13.841; average spacing: 5.273; total length: 8090588999 [M::mm_idx_gen::24348.48945.97] collected minimizers [M::mm_idx_gen::24359.96545.96] sorted minimizers [M::main::24359.96545.96] loaded/built the index for 118 target sequence(s) [M::mm_mapopt_update::24359.96545.96] mid_occ = 3923 [M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 118 [M::mm_idx_stat::24361.77445.96] distinct minimizers: 109911408 (33.83% are singletons); average occurrences: 13.819; average spacing: 5.271; total length: 8005691022 [M::mm_idx_gen::24492.50245.72] collected minimizers [M::mm_idx_gen::24501.01445.71] sorted minimizers [M::main::24501.01445.71] loaded/built the index for 158 target sequence(s) [M::mm_mapopt_update::24501.01445.71] mid_occ = 3923 [M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 158 [M::mm_idx_stat::24502.91245.71] distinct minimizers: 109785241 (34.04% are singletons); average occurrences: 13.898; average spacing: 5.273; total length: 8046380483 [M::mm_idx_gen::24636.64345.47] collected minimizers [M::mm_idx_gen::24649.40145.46] sorted minimizers [M::main::24649.40145.46] loaded/built the index for 711 target sequence(s) [M::mm_mapopt_update::24649.40145.46] mid_occ = 3923 [M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 711 [M::mm_idx_stat::24651.52445.45] distinct minimizers: 109316989 (33.82% are singletons); average occurrences: 13.963; average spacing: 5.273; total length: 8048423058 [M::mm_idx_gen::24665.29045.43] collected minimizers [M::mm_idx_gen::24665.60945.43] sorted minimizers [M::main::24665.60945.43] loaded/built the index for 9578 target sequence(s) [M::mm_mapopt_update::24665.60945.43] mid_occ = 3923 [M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 9578 [M::mm_idx_stat::24665.73945.43] distinct minimizers: 10233182 (68.18% are singletons); average occurrences: 10.239; average spacing: 5.333; total length: 558709858 [WARNING] failed to parse the FASTA/FASTQ record next to '�'. Continue anyway. minimap2: format.c:449: write_sam_cigar: Assertion `clip_len[0] < qlen && clip_len[1] < qlen' failed. [W::sam_parse1] empty query name [W::sam_read1_sam] Parse error at line 10900 Traceback (most recent call last): File "/work/home/gonglei/miniconda3/envs/samtools/bin/pore-c-py", line 8, in sys.exit(run_main()) File "/work/home/gonglei/miniconda3/envs/samtools/lib/python3.10/site-packages/pore_c_py/main.py", line 408, in run_main args.func(args) File "/work/home/gonglei/miniconda3/envs/samtools/lib/python3.10/site-packages/pore_c_py/main.py", line 316, in annotate_bam for concat_walk in annotate.annotate_alignments(inbam): File "/work/home/gonglei/miniconda3/envs/samtools/lib/python3.10/site-packages/pore_c_py/annotate.py", line 67, in annotate_alignments for concat_id, aligns in groupby( File "pysam/libcalignmentfile.pyx", line 2211, in pysam.libcalignmentfile.IteratorRowAll.next OSError: truncated file

For the above question, could you provide some debugging suggestions? Thank you for your valuable time and assistance. I sincerely look forward to your response!

sarahjeeeze commented 5 months ago

Hi, I am sure i have seen this error before check the pysam version in your environment is pysam>=0.21.0

GLking123 commented 5 months ago

Dear author: pysam version in my environment is pysam == 0.20.0 conda list pysam 0.20.0 pypi_0 pypi

Also, is the following code correct? ################################# workdir=/work/home/gonglei/data/new_new_minimap refdir=/work/home/gonglei/data/new_data/A.fasta fqdir=/work/home/gonglei/data/11.8porec1/pass.fq.gz ENZYME="NlaIII" OUTPUT="all" source /work/home/gonglei/miniconda3/etc/profile.d/conda.sh conda activate samtools export PATH="/work/home/gonglei/miniconda3/envs/samtools/bin:$PATH" SOFTWARE_PATH="/work/home/gonglei/miniconda3/envs/samtools/bin/" cd $workdir pore-c-py digest "${fqdir}" "${ENZYME}" | samtools fastq -T '*' | /work/home/gonglei/minimap2/minimap2 -ay -t 50 -x map-ont --split-prefix tmp "${refdir}" - | pore-c-py annotate - "${OUTPUT}" --monomers --stdout --summary --chromunity | tee "${OUTPUT}.ns.bam" | samtools sort --write-index -o "${OUTPUT}.cs.bam" - samtools index "${OUTPUT}.ns.bam

For the above question, could you provide some debugging suggestions? Thank you for your valuable time and assistance. I sincerely look forward to your response!

cjw85 commented 4 months ago

I suspect what has happened here is that minimap2 exhausted your machine's memory whilst aligning reads, this would cause it to produce incomplete output and the subsequent processes to fail. You can try performing the steps of the pipeline independently (producing temporary files) to check this assertion.

You may wish to add the options --cap-kalloc 100m --cap-sw-mem 50m to the minimap2 command.

GLking123 commented 4 months ago

refdir=/work/home/gonglei/data/new_data/A.fasta fqdir=/work/home/gonglei/data/11.8porec1/pass.fq.gz

minimap2 --split-prefix tmp -t 50 -x map-ont -a $refdir $fqdir

Directly run minimap2.

/work/home/gonglei/minimap2/minimap2 -ay -t 50 -x map-ont --split-prefix tmp "${refdir}"

After locating the restriction sites, running minimap2 failed. My Pore-C files end with .fq.gz, not with .bam files.

I'll try adding --cap-kalloc 100m --cap-sw-mem 50m to the minimap2 command.

sarahjeeeze commented 3 months ago

Did you get it run in the end? We have made a few updates including the above suggestion in the latest release so you might like to try the latest version 1.2.0

GLking123 commented 2 months ago

I'll give it a try now. Thank you for your suggestion.