how to align the results and call site-level methylation frequency?

[PengNi]

Dear the authors,

For R10.4.1, the results in predictions.tsv only contains read_id pos prob, so are there any scripts for aligning the results to reference and calling site-level frequency?

Thank you!

Best,
Peng

[chenxi-bio]

The read_id and prob can be converted into MM and ML tags and written into the BAM file.

[rjmat97]

It seems like they have the scripts for everything already present in the repo.

I discovered this extract_ref_pos.py script that generates a read-position -> reference_position map tsv file. This has to be intersected with the predictios.tsv file that rock fish generated. And then the file is ready for aggregation using the 'calculate_site_freq.py' script command.

#Commands to run (after running 'rockfish inference'):

##step1: generate read map

conda activate rockfish
#same bam file as the one used for inferencing
python rockfish/scripts/extract_ref_pos.py --workers ${threads} ${bam_file} ${reference_file} > reference_map.tsv

##step2: intersecting the read map with the predictions file

#performing the intersection:
r = pl.read_csv('reference_map.tsv', separator='\t', has_header=False,  new_columns=['read_id', 'pos', 'ref', 'chrom'])
d = pl.read_csv('rpredictions.tsv', separator='\t')
d = d.join(r, how='inner', on=['read_id', 'pos'])
d = d.select(['read_id', 'chrom', 'ref', 'prob']).rename({'ref': 'pos', 'chrom': 'ctg'})
d.write_csv('predictions_mapped.tsv', separator='\t')

Although the calculate_site_freq.py script is available in the repository, I was unable to run it as is, when examining the code it doesnt seem like the syntax was correct, so I've taken the liberty to edit it inorder to make it work.

def parse_rockfish(rockfish, get_predicted=False, type_suffix=''):
    if get_predicted:
        value_col = pl.when(pl.col('prob') > 0.5).then(1).otherwise(0).alias(
            f'Rockfish{type_suffix}').cast(pl.UInt8)
    else:
        value_col = pl.col('prob').alias(f'prob{type_suffix}').cast(pl.Float32)
    rockfish = rockfish.select([
        pl.col('read_id'),
        pl.col('ctg').cast(pl.Categorical),
        pl.col('pos').cast(pl.UInt32), value_col
    ])
    return rockfish
def main(df):
    N_NOMOD_COL = (pl.col('prob') <= 0.5).sum().alias('N_NOMOD_COL')
    N_MOD_COL = (pl.col('prob') > 0.5).sum().alias('N_MOD_COL')
    FREQ_COL = (pl.col('N_MOD_COL') /
                (pl.col('N_NOMOD_COL') + pl.col('N_MOD_COL'))).alias('freq')
    END_COL = (pl.col('pos') + 1).cast(pl.UInt32).alias('end')
    read_level = parse_rockfish(df)
    read_level = read_level.with_columns((pl.col('prob') - 0.5).abs().alias('abs')) \
                           .sort('abs', descending=True) \
                           .group_by(['ctg', 'pos']) \
                           .head(n=1000) \
                           .group_by(['ctg', 'pos']) \
                           .agg([N_MOD_COL, N_NOMOD_COL]) \
                           .sort(['ctg', 'pos'])#.collect()
    read_level = read_level.select([
        'ctg',
        pl.col('pos').alias('start'), END_COL, FREQ_COL, 'N_MOD_COL', 'N_NOMOD_COL'
    ])
    return read_level
d = main(df)
d.write_csv('predictions_aggregated.tsv', separator='\t')

This output can then be used for downstream analysis