RawHash integration available?

[maximilianmordig]

Hi
Thank you for the nice paper.

How did you integrate the segmentations into RawHash/Sigmoni? Is the code available?

Is the internal R10.4.1 Zymo dataset available?

[BakicSara]

Hi @maximilianmordig,

We have adapted the RawHash and Sigmoni code to simply parse a file that contains predicted segmentation for the sake of evaluation. I can provide that code to you if you're still interested.

We also have the pod5 files corresponding to R10.4.1 Zymo available, but we couldn't deposit them anywhere. I can provide the files to you!

[canfirtina]

Just hijacking the issue here if you do not mind @BakicSara :-)

It is also now possible to provide external segmentation information as either 1) series of event values, 2) segmentation points (peaks), 3) or move table information as a tsv file in RawHash2. You can also extract the segmentation information that RawHash2 generates if you want to evaluate them externally.

There is also a new repository that we are currently improving to include these experimental features so that the community can easily test their techniques by using RawHash2 as a framework: https://github.com/STORMgroup/RawHash2/tree/main?tab=readme-ov-file#external-segmentation

[jelber2]

I would like to chime in that I made some changes to Campolina via my fork at https://github.com/jelber2/Campolina , and you can run it as (after following the other installation steps in this repo):

python3 ~/git/Campolina/inference.py --pod5_dir ./ --model_path ~/git/Campolina/weights/R10_model.pth \
--tgt_dir ./ --mode analysis --bs 64

output is test_multithread_events.parquet

this script

parquet_to_tsv.py

import pyarrow.parquet as pq
pf = pq.ParquetFile('test_multithread_events.parquet')
current_id = None
current_starts = []
with open('test_multithread_events.tsv', 'w') as f:
    for batch in pf.iter_batches(batch_size=500_000):
        rows = batch.to_pydict()
        for rid, start in zip(rows['read_id'], rows['event_start']):
            if rid != current_id:
                if current_id is not None:
                    f.write(current_id + '\t' + '\t'.join(str(s) for s in current_starts) + '\n')
                current_id = rid
                current_starts = [start]
            else:
                current_starts.append(start)
    # flush the last read
    if current_id is not None:
        f.write(current_id + '\t' + '\t'.join(str(s) for s in current_starts) + '\n')

can be run as:

python3 parquet_to_tsv.py

output is test_multithread_events.tsv

that can be read by rawhash2, but my problem is that when I add --no-event-detection --events-file test_multithread_events.tsv to my rawhash2 command, I only get a PAF file with all unaligned reads. I used a command such as:

~/git/rawhash2-2/bin/rawhash2 -x ava-sensitive --r10 \
-p ~/git/rawhash2/extern/kmer_models/dna_r10.4.1_e8.2_400bps/9mer_levels_v1.txt -t10  \
--no-event-detection --events-file test_multithread_events.tsv chr21_2.pod5.ind \
chr21_2.pod5 2> chr21_2.pod5.rawhash2_align.txt |pigz -p 10 > chr21_2.paf.gz &

You can download the files to test for yourself @canfirtina

# get pod5 file
wget 'https://www.dropbox.com/scl/fi/zaxffqhzz614las1gdsfc/chr21_2.pod5?rlkey=byct2ux4490und0xz4zocdp47&st=y7yi21oo' -O chr21_2.pod5

# get parquet file of events from Campolina
wget 'https://www.dropbox.com/scl/fi/b37z3ygc63n07tonocfe6/chr21_2.pod5.events.parquet?rlkey=alzp836p5zp9ulii5nemf0lil&st=3nu6ca8w' -O chr21_2.pod5.events.parquet

# rename for script compatibility
mv chr21_2.pod5.events.parquet test_multithread_events.parquet

# convert parquet to tsv
python3 parquet_to_tsv.py

# make rawhash2 index
rawhash2 -x ava-sensitive --r10 \
-p ~/git/rawhash2/extern/kmer_models/dna_r10.4.1_e8.2_400bps/9mer_levels_v1.txt \
-t10 --no-event-detection --events-file test_multithread_events.tsv \
-d chr21_2.pod5.ind chr21_2.pod5 2> chr21_2.pod5.rawhash2_index.txt

# run rawsamble to overlap reads
rawhash2 -x ava-sensitive --r10 \
-p ~/git/rawhash2/extern/kmer_models/dna_r10.4.1_e8.2_400bps/9mer_levels_v1.txt  \
-t10 --no-event-detection \
--events-file test_multithread_events.tsv chr21_2.pod5.ind chr21_2.pod5 \
2> chr21_2.pod5.rawhash2_align.txt |pigz -p 10 > chr21_2.paf.gz

Best,
Jean Elbers

[jelber2]

--events-file is now working with the converted output of campolina with #jelber2/RawHash2@664f4ed, but it may be that the converted output of campolina should be used with --peaks-file and perhaps I was confused.

[canfirtina]

Hi @jelber2,

Thanks for testing this out and sharing the detailed workflow!

Campolina's output should be used with --peaks-file, not --events-file as Campolina outputs segmentation boundary positions (e.g., 42, 128, 215, ...). This is actually what --peaks-file expects.

--events-file expects pre-computed event values (e.g., -0.52, 0.13, 1.24, ...). When your parquet_to_tsv.py passes the integer positions (42, 128, 215) to --events-file, RawHash2 interprets them as signal amplitudes, which is why you originally got all unaligned reads.

Your commit 664f4ed works around this by converting the float values back to positions inside the --events-file code path, but the cleaner approach is to simply use the right option:

# Index
rawhash2 -x ava-sensitive --r10 \
  -p ~/git/rawhash2/extern/kmer_models/dna_r10.4.1_e8.2_400bps/9mer_levels_v1.txt \
  -t10 --peaks-file test_multithread_events.tsv \
  -d chr21_2.pod5.ind chr21_2.pod5 2> chr21_2.pod5.rawhash2_index.txt

# Map
rawhash2 -x ava-sensitive --r10 \
  -p ~/git/rawhash2/extern/kmer_models/dna_r10.4.1_e8.2_400bps/9mer_levels_v1.txt \
  -t10 --peaks-file test_multithread_events.tsv \
  chr21_2.pod5.ind chr21_2.pod5 \
  2> chr21_2.pod5.rawhash2_align.txt | pigz -p 10 > chr21_2.paf.gz

Some more points:

• RawHash2 reads the raw signal itself, computes segment means between your peaks, and normalizes.
• sig-diff quantization filtering is automatically disabled when segmentation information is externally provided.
• --no-event-detection is not needed (RawHash2 skips its own peak detection and uses yours instead, but still runs event generation from the peaks)
• Providing external segmentations is an experimental feature. If the peaks list is too long (e.g., >100K peaks per read), there may be issues related to fgets() but we are planning to fix that as well.

Also if there are further RawHash2-related issues, feel free to create issues here: https://github.com/STORMgroup/RawHash2/issues.