Error evaluating base qualities with `read_analysis.py`

[capoony]

Hi all,

when running read_analysis.py with the --fastq flag set, I get the following SyntaxWarnings and an error at the end, which apparently results in no output of the base-quality error profile files. Accordingly it is subsequently not possible to run the simulator.py script with the --fastq option since these files are required.

Thanks for your help

Best, Martin

(nanosim) (anaconda3) [mkapun@nhm-phylo2 results]$ read_analysis.py genome \
>     -i /media/inter/mkapun/projects/HAPLOTYPES-Benchmark/simulations/data/COX1.fastq \
>     -rg /media/inter/mkapun/projects/HAPLOTYPES-Benchmark/simulations/data/Stor1_cox1.fa \
>     -o /media/inter/mkapun/projects/HAPLOTYPES-Benchmark/simulations/data/COX1_training \
>     --fastq \
>     -t 100
/opt/anaconda3/envs/nanosim/bin/besthit_to_histogram.py:46: SyntaxWarning: invalid escape sequence '\*'
  for item in re.findall('(:[0-9]+|\*[a-z][a-z]|[=\+\-][A-Za-z]+)', cs_string):
/opt/anaconda3/envs/nanosim/bin/besthit_to_histogram.py:84: SyntaxWarning: invalid escape sequence '\d'
  cigar = re.findall('(\d+)([MIDSHX=])', cigar_str)  # Don't find (\d+)N since those are introns
/opt/anaconda3/envs/nanosim/bin/model_base_qualities.py:25: SyntaxWarning: invalid escape sequence '\*'
  for item in re.findall('(:[0-9]+|\*[a-z][a-z]|[=\+\-][A-Za-z]+)', cs_string):

Running the code with following parameters:

infile /media/inter/mkapun/projects/HAPLOTYPES-Benchmark/simulations/data/COX1.fastq
ref_g /media/inter/mkapun/projects/HAPLOTYPES-Benchmark/simulations/data/Stor1_cox1.fa
aligner minimap2
g_alnm 
prefix /media/inter/mkapun/projects/HAPLOTYPES-Benchmark/simulations/data/COX1_training
num_threads 100
model_fit True
chimeric False
homopolymer False
fastq True
2024-09-08 09:12:36: Read pre-process
2024-09-08 09:12:36: Alignment with minimap2
[M::mm_idx_gen::0.001*2.85] collected minimizers
[M::mm_idx_gen::0.012*3.86] sorted minimizers
[M::main::0.012*3.85] loaded/built the index for 1 target sequence(s)
[M::mm_mapopt_update::0.012*3.83] mid_occ = 10
[M::mm_idx_stat] kmer size: 15; skip: 10; is_hpc: 0; #seq: 1
[M::mm_idx_stat::0.013*3.82] distinct minimizers: 120 (100.00% are singletons); average occurrences: 1.000; average spacing: 5.383; total length: 646
[M::worker_pipeline::0.032*4.85] mapped 500 sequences
[M::main] Version: 2.28-r1209
[M::main] CMD: minimap2 --cs -ax map-ont -t 100 /media/inter/mkapun/projects/HAPLOTYPES-Benchmark/simulations/data/Stor1_cox1.fa /media/inter/mkapun/projects/HAPLOTYPES-Benchmark/simulations/data/COX1_training_processed.fasta.gz
[M::main] Real time: 0.033 sec; CPU: 0.158 sec; Peak RSS: 0.019 GB
2024-09-08 09:12:36: Processing alignment file: bam
2024-09-08 09:12:36: Aligned reads analysis
2024-09-08 09:12:36: Computing KDE
2024-09-08 09:12:36: Computing KDE for aligned region
2024-09-08 09:12:36: Computing KDE for aligned reads
2024-09-08 09:12:36: Computing KDE for unaligned length
2024-09-08 09:12:36: Computing KDE for ht ratio
2024-09-08 09:12:36: Unaligned reads analysis
2024-09-08 09:12:36: match and error models
2024-09-08 09:12:36: Model fitting
2024-09-08 09:12:36: Mismatch fitting start
Mismatch parameters: 0.04795745255153716 0.9772561838132541 0.9926567732360942 Residual 0.001109464694009299
2024-09-08 09:12:37: Mismatch fitting done
2024-09-08 09:12:37: Insertion fitting start
Insertion parameters: 1.9821603228471387 4.33796857852048 0.5833991237731329 0.0511711826739132 Residual 0.00818994018516872
2024-09-08 09:12:46: Insertion fitting done
2024-09-08 09:12:46: Deletion fitting start
Deletion parameters: 2.5109434849874015 6.7579121539864735 0.5943697519101852 0.08959092674359895 Residual 0.002620869492597655
2024-09-08 09:12:55: Deletion fitting done
2024-09-08 09:12:55: Finished!
2024-09-08 09:12:55: Analyzing base qualities and estimating model parameters
2024-09-08 09:12:55: Parsing alignment file for base qualities relative to matches     and each error type
Traceback (most recent call last):
  File "/opt/anaconda3/envs/nanosim/bin/read_analysis.py", line 879, in <module>
    main()
  File "/opt/anaconda3/envs/nanosim/bin/read_analysis.py", line 875, in main
    model_base_qual.model_base_qualities(prefix + "_primary." + alnm_ext, prefix, unaligned_base_qualities)
  File "/opt/anaconda3/envs/nanosim/bin/model_base_qualities.py", line 108, in model_base_qualities
    quals_per_type = analyze_aligned_base_qualities(alnm_file)
                     ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
  File "/opt/anaconda3/envs/nanosim/bin/model_base_qualities.py", line 66, in analyze_aligned_base_qualities
    quals = alnm.query_alignment_qualities.tolist()
            ^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^^
AttributeError: 'NoneType' object has no attribute 'tolist'

[lcoombe]

Thanks for the report, Martin.

I was able to reproduce the AttributeError on my end - I believe it was due to a bug where in genome mode, the reads aligned to the reference using minimap2 were in fasta format, despite you specifying --fastq.

I have pushed a quick fix in this branch: https://github.com/bcgsc/NanoSim/tree/fastq_fix This fixed the issue on my end (using equivalent parameters to your command), but I'd greatly appreciate if you could test to see if it fixed the error on your end as well!

Thanks! Lauren

[lcoombe]

Update - the fixes were merged to master branch in https://github.com/bcgsc/NanoSim/pull/222

[capoony]

Thanks a lot, works like a charm :-)