homopolish modpolish Index error

[giriarteS]

Can modpolish be used after racon+medaka? I am trying to polish my assembly (flye+6 rounds of racon+medaka) with modpolish but I get the following error: "python3 /mnt/DATA/git_repos/homopolish/homopolish.py modpolish -a medaka/consensus.fasta -q /mnt/Backup/Genomes_test/raw_dna/minion/F.verticillioides/MRC826E/MRC826E_trim_fastq.gz -s /mnt/DATA/git_repos/homopolish/fungi.msh -o homopolish -t 50

Downloaded GCA_003316995.2_ASM331699v2_genomic.fna.gz Downloaded GCA_003316975.2_ASM331697v2_genomic.fna.gz Downloaded GCA_011035275.1_ASM1103527v1_genomic.fna.gz Downloaded GCA_011036225.1_ASM1103622v1_genomic.fna.gz Downloaded GCA_000149555.1_ASM14955v1_genomic.fna.gz Downloaded GCF_000149555.1_ASM14955v1_genomic.fna.gz Downloaded GCA_003317015.2_ASM331701v2_genomic.fna.gz Downloaded GCA_011036015.1_ASM1103601v1_genomic.fna.gz Downloaded GCA_011036905.1_ASM1103690v1_genomic.fna.gz

Reference = [/mnt/Backup/Genomes_test/assembly/flye/F.verticillioides/MRC826E/homopolish/debug/homologous_sequences/GCA_003316995.2_ASM331699v2_genomic.fna.gz.fna.gz, /mnt/Backup/Genomes_test/assembly/flye/F.verticillioides/MRC826E/homopolish/debug/homologous_sequences/GCA_003316975.2_ASM331697v2_genomic.fna.gz.fna.gz, /mnt/Backup/Genomes_test/assembly/flye/F.verticillioides/MRC826E/homopolish/debug/homologous_sequences/GCA_011035275.1_ASM1103527v1_genomic.fna.gz.fna.gz, /mnt/Backup/Genomes_test/assembly/flye/F.verticillioides/MRC826E/homopolish/debug/homologous_sequences/GCA_011036225.1_ASM1103622v1_genomic.fna.gz.fna.gz, /mnt/Backup/Genomes_test/assembly/flye/F.verticillioides/MRC826E/homopolish/debug/homologous_sequences/GCA_000149555.1_ASM14955v1_genomic.fna.gz.fna.gz, /mnt/Backup/Genomes_test/assembly/flye/F.verticillioides/MRC826E/homopolish/debug/homologous_sequences/GCF_000149555.1_ASM14955v1_genomic.fna.gz.fna.gz, /mnt/Backup/Genomes_test/assembly/flye/F.verticillioides/MRC826E/homopolish/debug/homologous_sequences/GCA_003317015.2_ASM331701v2_genomic.fna.gz.fna.gz, /mnt/Backup/Genomes_test/assembly/flye/F.verticillioides/MRC826E/homopolish/debug/homologous_sequences/GCA_011036015.1_ASM1103601v1_genomic.fna.gz.fna.gz, /mnt/Backup/Genomes_test/assembly/flye/F.verticillioides/MRC826E/homopolish/debug/homologous_sequences/GCA_011036905.1_ASM1103690v1_genomic.fna.gz.fna.gz] Query = [medaka/consensus.fasta] Kmer size = 16 Fragment length = 3000 Threads = 1 ANI output file = homopolish/debug/ANI.txt

INFO [thread 0], skch::main, Count of threads executing parallel_for : 1 INFO [thread 0], skch::Sketch::build, window size for minimizer sampling = 24 INFO [thread 0], skch::Sketch::build, minimizers picked from reference = 30851510 INFO [thread 0], skch::Sketch::index, unique minimizers = 4853959 INFO [thread 0], skch::Sketch::computeFreqHist, Frequency histogram of minimizers = (1, 1098191) ... (7088, 1) INFO [thread 0], skch::Sketch::computeFreqHist, consider all minimizers during lookup. INFO [thread 0], skch::main, Time spent sketching the reference : 26.7354 sec INFO [thread 0], skch::main, Time spent mapping fragments in query #1 : 75.2743 sec INFO [thread 0], skch::main, Time spent post mapping : 0.0381457 sec INFO [thread 0], skch::main, ready to exit the loop INFO, skch::main, parallel_for execution finished

load Homo

Reference = [homopolish/debug/All_homologous_sequences.fna.gz] Query = [medaka/consensus.fasta] Kmer size = 16 Fragment length = 3000 Threads = 1 ANI output file = homopolish/debug/truth_ANI.txt

INFO [thread 0], skch::main, Count of threads executing parallel_for : 1 INFO [thread 0], skch::Sketch::build, window size for minimizer sampling = 24 INFO [thread 0], skch::Sketch::build, minimizers picked from reference = 30851510 INFO [thread 0], skch::Sketch::index, unique minimizers = 4853959 INFO [thread 0], skch::Sketch::computeFreqHist, Frequency histogram of minimizers = (1, 1098191) ... (7088, 1) INFO [thread 0], skch::Sketch::computeFreqHist, consider all minimizers during lookup. INFO [thread 0], skch::main, Time spent sketching the reference : 26.885 sec INFO [thread 0], skch::main, Time spent mapping fragments in query #1 : 75.9234 sec INFO [thread 0], skch::main, Time spent post mapping : 0.0205193 sec INFO [thread 0], skch::main, ready to exit the loop INFO, skch::main, parallel_for execution finished [M::mm_idx_gen::0.8991.00] collected minimizers [M::mm_idx_gen::1.0612.05] sorted minimizers [M::main::1.0612.05] loaded/built the index for 83 target sequence(s) [M::mm_mapopt_update::1.1341.98] mid_occ = 100 [M::mm_idx_stat] kmer size: 19; skip: 19; is_hpc: 0; #seq: 83 [M::mm_idx_stat::1.1761.94] distinct minimizers: 4248851 (99.24% are singletons); average occurrences: 1.018; average spacing: 9.991 [M::worker_pipeline::19.96514.14] mapped 49912 sequences [M::main] Version: 2.17-r941 [M::main] CMD: minimap2 -cx asm5 --cs=long -t 50 medaka/consensus.fasta homopolish/debug/All_homologous_sequences.fna.gz [M::main] Real time: 20.085 sec; CPU: 282.508 sec; Peak RSS: 4.439 GB 127.97561192512512 Traceback (most recent call last): File "/mnt/DATA/git_repos/homopolish/homopolish.py", line 102, in main() File "/mnt/DATA/git_repos/homopolish/homopolish.py", line 77, in main getPos(fixData,FLAGS.debug) File "/mnt/DATA/git_repos/homopolish/modules/mod_polish.py", line 80, in getPos H_misAry,H_AllAry = getPileUpAry(fixData,fixData.output_dir+"debug",fixData.output_dir+"debug/All_homologous_sequences.fna.gz")
File "/mnt/DATA/git_repos/homopolish/modules/mod_polish.py", line 146, in getPileUpAry misAry,posData_ary = MismatchPileup(fixData.output_dir+"debug/truth.paf",len(fixData.seq))#get SNP position ATCG array File "/mnt/DATA/git_repos/homopolish/modules/mod_polish.py", line 440, in MismatchPileup coverage[start_pos] += 1 IndexError: index 2852431 is out of bounds for axis 0 with size 2852431"

Gloria

[ythuang0522]

Hi @giriarteS,

We haven't tested modpolish on fungi yet. Have you ever successfully run the original homopolish polish? If so, can you give us your fasta genome and fastq reads for debugging? modpolish is a new feature for genomes extensively edited by modifications untrained in the ONT basecalling. If your genomes are not extensively modified, you don't need to run modpolish.

[giriarteS]

I run homopolish polish but I got the following for each contig: python3 /mnt/DATA/git_repos/homopolish/homopolish.py polish -a medaka/consensus.fasta -m R9.4.pkl -s /mnt/DATA/git_repos/homopolish/fungi.msh -o homopolish -t 50

[2022/05/25 17:24] INFO: Stage: Select closely-related genomes TIME Select closely-related genomes: 0 MINS 0 SECS. This contig contig_75 closely-related genome is less than 5, not to polish... [2022/05/25 17:24] INFO: RUN-ID: contig_197

[giriarteS]

I successfully run homopolish polish: homopolish.py polish -a medaka/consensus.fasta -m R9.4.pkl -l reference/F.verticillioides/GCF_000149555.1_ASM14955v1_genomic.fna.gz -o homopolish -t 25

[ythuang0522]

We found the error was due to an insufficient number of related genomes (<5), a bug in modpolish but not original homopolish. It looks to us there are six in GenBank but only one in RefSeq. Therefore, as you have done, the local db option -l can specify the related genomes downloaded on your own. However, I would suggest you concatenate all the six related genomes into one fasta not just provided only one.

[giriarteS]

I polished my genomes with several rounds of racon, then medaka, then homopolish. If I want to use illumina reads to polish with Pilon, should homopolish be run before or after pilon?

[ythuang0522]

It was trained on Medaka-polished genomes. We would recommend running homopolish right after Medaka and then followed by pilon.