Write altered reads to a slow5 file

[kaltinel]

Hi,

I would like to take a subset of reads for certain positions (for ex. read 1, read 2 and read 3), and update their signal entries with the signals from these new range of positions. Then I would like to write these reads with altered signal entries to a slow5 file. For this, I am implementing the following:

s5 = pyslow5.Open(filename,'r')
readID = "e9ew1ea4-0589-4c0c-435v-8b6f92f524e"
read_1 = s5.get_read(readID, pA=True)

Then, I would like to assign these new subset of signals as their 'own' signal:

read_1['signal'] = read_1['signal'][21177:21283]  # 21177: signal start_idx and  21283: signal end_idx
# do the same for read2 and read3

However, I am unsure how to write them to a slow5 file (which would be consisting of these 3 reads and their headers etc) I was thinking either to append the reads to a df, and then writing this df to a slow5. For ex.:

df = pd.DataFrame.from_dict(read_1)
df = df.append(pd.DataFrame.from_dict(read2))
df = df.append(pd.DataFrame.from_dict(read3))
s5_new = slow5.Open(file, 'w')
s5_new.write_record(df)
s5_new.close()

But in this case, shall I append the headers every time before a read? Or, would this way of appending of reads work?

I would be glad to have your help on the matter.

Thank you very much for your time.

[Psy-Fer]

Hello,

I'd start by having a look at the pyslow5 api to see how to write a slow5 file

https://github.com/hasindu2008/slow5lib/blob/master/docs/pyslow5_api/pyslow5.md

To see an example of reading a slow5 file, modifying it, and writing it to a new file, along with headers, please see this issue where I have previously given example code

https://github.com/hasindu2008/slow5lib/issues/70#issuecomment-1594094571

Let me know if you have any further troubles, or something isn't working.

James

[Psy-Fer]

I've given a code snippet below that might help you get to where you want to be

Let me know if you have any troubles or want some more explanation.

import pyslow5
import numpy as np

s5 = pyslow5.Open(filename,'r')
readID = "e9ew1ea4-0589-4c0c-435v-8b6f92f524e"
# if you want to just take a slice of the raw signal and put it back that way, don't convert to pA
# including aux get's the full record data, not just the primary data fields
read_1 = s5.get_read(readID, aux='all')

header, end_reason_labels = s5.get_empty_header(aux=True)
header_original = s5.get_all_headers()

# this is to update the end_reason_labels because ONT changed them
new_end_reason_labels = ['unknown', 'partial', 'mux_change', 'unblock_mux_change', 'data_service_unblock_mux_change', 'signal_positive', 'signal_negative']

# this copies the original header to put into the new file.
for i in header_original:
    if i in header:
        if header_original[i] is None:
            continue
        else:
            header[i] = header_original[i]

# open a new file to write the reads into
W = pyslow5.Open('modified_reads.blow5','w')
# write the header. Note, it's not actually written to disk until the first read is written
ret = W.write_header(header, end_reason_labels=new_end_reason_labels)
# print the return value
print("ret: write_header(): {}".format(ret))

# get an empty record and aux
record, aux = W.get_empty_record(aux=True)
for i in read_1:
    if i in record:
        if i == "signal":
            # insert modified signal values
            record[i] = read_1[i][21177:21283] 
        else: 
            record[i] = read[i]
    if i in aux:
        aux[i] = read[i]
print(record)
print(aux)
ret = W.write_record_batch(records, aux)
# print return code
print("ret: write_record(): {}".format(ret))

s5.close()
W.close()

[kaltinel]

Hi @Psy-Fer ,

Thank you very much for the code and the comments. I will try your suggestion and update you with the outcome. I appreciate your kind help and time.

[Psy-Fer]

Just adding a note that we solved this in person.