I am trying to understand the sanger sequencing ABI/AB1 file format better, and extract base calls from given signal intensities over time.
As I understand, reading in a raw AB1/ABI file into python, I can access the channel corrected values for the different bases in ['DATA9'] to ['DATA12'], and understand which bases belong to a specific channel through ['FWO_1'] i.e either ordered GATC or ATGC etc.
However, the file I am looking at has identified 412 bases as the output from SnapGene, yet has 4950 data records for each channel in the raw ['DATA9'] to ['DATA12']. How do you convert those 4950 data records into the bases correctly? Is there some form of normalisation - e.g every 10 raw records gives data for 1 base position, or peak calling signal detection algorithms. And if so, do you take the average over 10 records, or the highest peak in the 10 records from the channels and the channel with the highest peak/average is the correct base? Do you start at the beginning or do you have some form of cutoff? Does it make sense to convert this figure into a quality score similar to NGS, and do you use the highest peaks in the 10 record interval or the average?
I've tried using scipy's signal.find_peaks() function (python) which can get me various peaks with a minimum height and a minimum distance to other peaks, currently using a minimum peak height of 150 and distance to other peaks of 5, but still getting around 20-30 more peaks than the bases identified.
Could anyone suggest any logic, methods, algorithms or resources that could help with translating the raw signal intensities for the 4 channels into a sequence of accurate base calls?
I know software exists to identify the bases, but that's not what I am looking for. I intend to do this from scratch to understand the workings better so I can write some custom processing algorithms.
I hope this question makes sense.
Thanks
