Improved method for identification of nucleotide modifications

Abstract

The invention relates to identifying one or more anomalous monomers in a polymer. This invention also relates to a method of preparing control data for use in a method of identifying anomalous monomers in a polymer. The invention also relates to an apparatus to perform the method and a system comprising this apparatus and a nucleic acid sequencer.

Description

[0001] Improved method for identification of nucleotide modifications

[0002] Field of the Invention

[0003] The present invention relates to a computer-implemented method for identifying one or more anomalous monomers in a polymer; and a computer program or computer-storage medium comprising instructions to perform the method. The invention also relates to an apparatus configured to perform the method and a system comprising this apparatus and a nucleic acid sequencing device.

[0004] Background

[0005] Cells are regulated by complex mechanisms including modifications to DNA, RNA, protein and other polymers to control their structure and function. These modifications play a critical part in many diseases.

[0006] Existing models for detection of polynucleotide modifications rely on knowing the modifications of interest which limits the identification of unknown modifications. However, it is believed that there are many modifications of these polymers which are currently unknown or impossible to observe in a given context.

[0007] Summary of the Invention

[0008] Current methods used in the art are focused on detecting specific modifications. The present invention instead detects any anomaly against a baseline of normal observation where, in context, an anomaly could be an as-yet unidentified modification or any observation that is differentiated from the canonical nucleotide call. This invention allows new, unknown modifications to be located in a polymer.

[0009] In a first aspect the invention provides a computer-implemented method for identifying one or more anomalous monomers in a polymer, the method comprising: a) obtaining a stream of input data comprising measurements for a plurality of consecutive monomers in the polymer; b) transforming the stream of input data in a) into a plurality of vector features; c) calculating the (a) distance between: i) each of the vector features in b); and ii) a plurality of vector features in a control data set, wherein the control data set was generated by transforming one or more streams of control data into one or more pluralities of vector features to form a control data set, wherein the control data comprises measurements for a plurality of consecutive monomers in one or more polymer(s); and d) based on the distances calculated in c), determining if a monomer associated with the one or more vector features in the polymer which is the input data is anomalous relative to the polymer(s) of the control data.

[0010] In a further aspect, the invention provides a computer-readable storage medium or a computer program comprising computer-executable instructions, which when executed by a computing system, are capable of causing the computing system to perform the method.

[0011] In a further aspect, the invention provides a computer-implemented method of preparing a control data set for use in a method of identifying one or more anomalous monomers in a polymer, the method comprising: a) obtaining one or more stream(s) of control data comprising measurements for a plurality of consecutive monomers in one or more polymer(s); and b) transforming the one or more stream(s) of control data into one or more pluralities of vector features to form a control data set, optionally wherein the control data set is normalised.

[0012] In a further aspect, the invention provides an apparatus comprising processing circuitry configured to perform the method for identifying one or more anomalous monomers in a polymer.

[0013] In a further aspect, the invention provides a system comprising: a) a nucleic acid sequencing apparatus; and b) the apparatus of claim 19, wherein the apparatus is operably connected to the sequencing apparatus and configured to receive a stream of input data from the sequencing apparatus.

[0014] Detailed description

[0015] Method for identifying one or more anomalous monomers in a polymer

[0016] The method is for identifying one or more anomalous monomers in a polymer. In this context, anomalous describes, and is defined to be, any monomer which is not found in the control data. That is, an anomaly, which is an occurrence of an anomalous monomer, is defined relative to the composition of the polymers used in the control data, and means a feature (which may be the presence or absence of a feature) not present in the control data. The control data may in particular contain modifications , i.e. non-canonical monomers. Then, what is searched for in the input data is anomaly compared to this control data, i.e. anomalies not present in the control data.

[0017] The anomaly may be a covalent modification. For example, an RNA / DNA modification.

[0018] Stream of data A stream of data is data that are emitted, or otherwise generated by some process, and which are capable of being measured or recorded either in a continuous or in a discrete, incremental manner.

[0019] Input data

[0020] The stream of input data comprises measurement for a plurality of consecutive measurements in a polymer. The stream of input data may comprise an electrical signal and / or an optical signal. For example, an ionic current, impedance or fluorescence.

[0021] The input may for example be the raw electric signal, for example at different sampling rates or the raw signal may be pre-processed, for example normalised and / or segmented as commonly applied in base calling software.

[0022] The method may further comprise the pre-step of sequencing the polymer (e.g. DNA or RNA) to provide the stream of input data or control data.

[0023] Control data

[0024] One or more streams of control data may be used. These streams are chosen to represent the variety of data that the user of the device regards as non-anomalous.

[0025] These are transformed into one or more pluralities (groups) of vector features to form a control data set. The measurements for the control data may be synthetic. By synthetic it is meant that it is not measured directly using physical, or real-world lab equipment but instead generated (e.g. in silica). That is, where the measurements are synthetic measurements, this data is synthetically derived measurements.

[0026] Or the control data used may have been as for the input data above, i.e. electrical signal and / or optical signal. Where actual measurements are taken, the polymers used for the control data measurements may be natural or synthetic polymers.

[0027] The one or more polymers used for the control data may be canonical polymers (standard or most common sequence). For example, where the polymer is RNA, the RNA polymers used for the control data may include no modifications (or additionally include exactly the modifications that are not of interest due to searching for further unknown modifications in the polymer, as described further below).

[0028] By unmodified means no additional modification, i.e. the monomer units are the canonical ones. For example, in the case of DNA, the unmodified monomer units are adenosine, thymidine, guanosine and cytidine. Modified forms may be nucleotides where any of the following covalent chemical modification is present, for example 5-methyl-cytosine (5mC ), 5-hydroxy methylcytosine (5hmC), and 6-methyl-adenosine (6mA). In the case of RNA, the unmodified monomer units are adenosine, uracil, guanosine and cytidine.

[0029] Alternatively, as the anomaly is relative to the polymer in the control data, to find unknown RNA modifications, polymer(s) comprising monomers with known modifications may be used as the control data. The method then compares the input data with this control data and locates monomers with new, unknown RNA modifications which are not in the control data.

[0030] Known RNA modifications include the following:

[0031] N6-methyladenosine (m6A)

[0032] Inosine (I)

[0033] 5-methylcytosine (m5C) pseudouridine (^P)

[0034] N6-dimethyladenosine (m6,2A)

[0035] 1-methylguanosine (m1 G)

[0036] 2’-0 methyladenosine (2-OMeA) would also be covered by 2’-0-methylation (Nm)

[0037] 7-methylguanosine (m7G)

[0038] 5-Methyluridine (m5U)

[0039] N1 -methyladenosine (m1A)

[0040] N4-acetylcytosine (ac4C)

[0041] 2-thiouridine (s2U) uridylation adenosine-to-inosine (A-to-l) RNA editing

[0042] Dihydrouridine (D)

[0043] 5-formylcytidine (C)

[0044] 2’-0-methylation (Nm)

[0045] Queuosine (Q)

[0046] 5-methylaminomethyl-2-thiouridine (mnm5s2U)

[0047] Wybutosine (yW)

[0048] 3- methylcytidine (m3C)

[0049] 5-methoxycarbonylmethyl-2-thiouridine (mcm5s2U)

[0050] 2-thiocytidine (s2C)

[0051] N2,N2-dimethylguanosine (m2,2G)

[0052] N6-isopentenyladenosine (i6A)

[0053] N6-threonylcarbamoyladenosine (t6A) 3-methyluridine (m3U)

[0054] Wyosine (imG)

[0055] 1 -methylinosine (ml I)

[0056] Monomers / polymer

[0057] The polymer comprises a sequence of monomers. The polymer may be a biopolymer. For example, the polymer may be a polynucleotide, for example DNA or RNA (for example mRNA or tRNA). Alternatively, the polymer may be a protein, where the monomer units are amino acids. Or the polymer may be a polysaccharide.

[0058] The monomers are individual units which make up the polymer. The stream of data is measured for a plurality of consecutive monomers. By consecutive is meant adjacent or next to each other in the polymer.

[0059] Vector features

[0060] The input data and the control data are transformed into vector features. By this is meant transforming the stream of measurements, e.g. electrical signal, into a plurality of vectors. This allows the signal to be represented by vectors. The transformation is chosen so that the vector features provide concise, informative and contextual information about the shape of the signal around the monomer location in the stream. This contextual information is represented in the form of vectors or tensors and includes features extracted from un-parametrized measurements. That is, the vectorisation tracks or follows the signal. The present invention therefore uses features extracted from the signature transformation or another differential equation controlled by the signal. That is, the vector features reflect the shape of a curve or stream derived from the input data. The features derived from the shape of a stream are features derived from an un-parameterized curve and are independent (largely) of any sampling of that curve. The vector features comprise systematic features derived from the local shape of the stream of input data. The signature transformation is therefore a signature based transform therefore describes using the entire signature transform or features extracted from the signature transformation, i.e. parts of the signature transform that are most informative about what constitutes 'normality’ and applying a distance to those. These can be extracted from the signal and distances between these vectors calculated.

[0061] The transforming steps therefore comprise a signature based transform rather than the use of other features from the signal, e.g. the use of trivial low dimensional summaries such as dwell time or median signal intensity. A signature transformation is applied to the stream of input data and the one or more streams of control data to form a plurality of vector features from the input data and the control data respectively, e.g. a signature based transform extracting certain features from the signature based transform. Vector features can then be extracted and compared between the input data and the control data, i.e. by calculating the distance between the input data and the control data.

[0062] By un-parameterized is meant that any two streams which are related by re-parametrization will produce identical vector features. By two streams being related by re-parametrization it is meant that each stream comprises exactly the same collection of values in exactly the same order but possibly asynchronously and possibly different rates, and such that one stream can be written in terms of the other by the composition of the other with a re-parameterisation function.

[0063] The input data and control data are transformed into vector features in a consistent manner. For example, one postprocessing step commonly applied to the raw nanopore current signal consists in dividing the raw data series into overlapping 5-mer-corresponding signal segments and transforming those segments into vector features. The resulting segmentation is depicted by the vertical dashed lines on Figure 2. The window in this segmentation may be called a position.

[0064] Distance

[0065] The distance is calculated by comparing each vector feature generated from the input data with vectors in the control data set. That is, with the set of features from the control dataset associated with the location of the feature from the input data. The distance may be large for some of these comparisons, where the measurement in the input data set is very different from a particular measurement in the control data set. However, if not an anomaly, there will be some data represented by vector features in the control data set which are very similar to those of the vector feature of the input data. On the other hand, if there is an anomaly at a particular point in the polymer, vector features transformed from this data will likely be significantly different from any vector features in the control data set and the difference calculated at c) in claim 1 will be larger, e.g. than the noise of natural variability as calculated by the thresholds. As a result, the monomer, e.g. nucleotide, at this position can be classified as an anomalous monomer, e.g. an RNA monomer with a modification not present in the control data set. This may be done by calculating the minimum distance for the vector feature (also referred to herein as the score).

[0066] This distance may be described using the Mahalanobis distance (also referred to as variance norm) to each element in the corpus (the distance from each vector generated from the input data to each vector in the control data set). This distance is invariant to the scale of the underlying streams.

[0067] Determining if the monomer is anomalous relative to the polymer of the control data

[0068] As explained above, the distance calculated in step c) of claim 1 can be used as a measure of the conformance of the vector in the input data to the corpus of normality calculated from the control data (represented by the vector features obtained from the control data). To calibrate the control data set, a calibration data set may be used. This may be a separate set from different control data, or it may be a subset of the control data. If a subset of the control data is to be used, this is extracted (i.e. split from) from the control data prior to the distance function being calculated.

[0069] The calibration data may be used to test how accurately the vector features of the control data set represent the corpus of normality. That is, how concentrated the vectors of the control data are relative to the distance and set the thresholds between normal and anomalous vector features.

[0070] The determination in step d) may comprise checking if the distances calculated are beyond thresholds for anomaly, wherein the thresholds for anomaly were determined by: i) obtaining calibration data which is representative of the control data; ii) ensuring the distance function in step c) is independent of the calibration data; iii) transforming the one or more stream(s) of calibration data into one or more pluralities of vector features; iv) calculating the distance between each of the vector features in iii) and the plurality of vector features in the control data set; v) using the distances calculated in iv) to set thresholds for anomaly for the input data; and vi) applying a threshold obtained from v) to the input data, optionally wherein the calibration data was extracted from the control data.

[0071] Calibration may be carried out by vectorising the calibration data in the same way as for the input data, i.e. a) transforming one or more stream(s) of calibration data into one or more pluralities of vector features;

[0072] Then calculating the distance between: b) each of the vector features in a) above; and a plurality of vector features in the control data set, wherein the control data set was generated by transforming one or more streams of control data into one or more pluralities of vector features to form a control data set, wherein the control data comprises measurements for a plurality of consecutive monomers in one or more polymer(s). A score may be determined for each vector feature formed from the calibration data (in the same way as carried out for the vectors of the input data). This score may be the minimum distance for that vector feature. Once determined, the collection of the scores across all points in the calibration data can be used to estimate the probability that, for a chosen numerical value, a vector sampled from the calibration data will have a score that exceeds this value,

[0073] This may be done by any of the following methods:

[0074] By fitting a parametric distribution to the calibration scores, and computing the probability mentioned above by evaluating the complementary cumulative distribution function at the chosen numerical value.

[0075] Compute p (probability) such that the empirical (1 - p)-quantile of the calibration scores is the input score. This is the same method as the first however the empirical cumulative distribution function is used (3).

[0076] In summary, the calibration is carried out by performing the same actions on the calibration data set as for the input data. This results in a set of distances for the calibration data (from the control data). The analysis of the distribution of these distances can tell us how the control data set represents normality. From this distribution the probability of any distance calculated being an anomaly can be estimated.

[0077] These normalised p-values can then be used to decide how unusual or exceptional the score is for the input data stream. Optionally further, p-values may be corrected for multiple testing using any method known in field, for example the false discovery rate.

[0078] Step v) of claim 11 may therefore further comprise: calculating the distribution of scores for the calibration data. Then using this distribution to calculate probability estimates for scores for the input data (i.e. distances, for example minimum distances). Then using these probability estimates to set thresholds for anomaly and comparing the input data distances (e.g. scores) with these thresholds..

[0079] This calibration step allows us to convert conformance scores (distances calculated from the calibration set from the control data set) into probability estimates. While the scores over different segments might not be directly comparable (a score of say 0.05 might be extreme with respect to the calibration set of a segment, but not with respect to another one), the probability estimates are. This threshold determination may be pre-determined, i.e. not calculated concurrently with performing the method of claim 1 . The control data set may also be predetermined. Therefore, the method of claim 1 is to compare the input data with these predetermined values in the software.

[0080] The threshold determinations (i.e. if any vector feature has a value beyond the threshold for anomaly) may be aggregated across one or more consecutive monomers.

[0081] The output from the method is therefore data indicating a position of one or more anomalous monomer(s) in the polymer.

[0082] Normalising the input data

[0083] The input data may be collected from a different source, e.g. sequencing machine, than the control data. The vector features of the input data may therefore be normalised so as to be on the same scale as those from the control data.

[0084] Nanopore

[0085] Nanopore technology uses a nanopore embedded in a membrane. Translocation of the polymer through a nanopore can be used to supply the consecutive measurements which form the data sets.

[0086] The membrane may be electrically resistant and a current may be applied across the nanopore with changes in electrical signal (e.g. current) detected from a polymer passing through a nanopore. Therefore, the data output by the nanopore sequencing data is an electrical signal which changes over time as different monomer units pass through the pore. Suitable nanopore systems include those described in Mackenzie and Argyropoulos (Micromachines 2023, 14, 459).

[0087] Preparing a control data set for use in a method of identifying anomalous monomers in a polymer

[0088] The features described above for claim 1 equally apply to claim 15, to preparation of the control data set. For example, the control data set can be normalised and calibrated as described above, i.e. by i) obtaining calibration data which is representative of the control data; ii) ensuring the distance function in step c) is independent of the calibration data; iii) transforming the one or more stream(s) of calibration data into one or more pluralities of vector features; iv) calculating the distance between each of the vector features in iii) and the plurality of vector features in the control data set; and v) using the distances calculated in iv) to set thresholds for anomaly.

[0089] Computer program and non-transitory media

[0090] By computer program is meant machine readable program instructions. These may be provided on a transitory medium such as a transmission medium or on a non-transitory medium such as a storage medium. Such machine readable instructions (computer program code) may be implemented in a high level procedural or object oriented programming language. However, the program(s) may be implemented in assembly or machine language, if desired. In any case, the language may be a compiled or interpreted language, and combined with hardware implementations. Program instructions may be executed on a single processor or on two or more processors in a distributed manner.

[0091] Therefore also included are one or more non-transitory computer readable media storing machine- readable instructions which, when executed, cause one or more processors to perform the method of any of claims 1 to 13, and as exemplified in Figure 7. A further description of a general process in accordance with the invention can also be found at Figure 9.

[0092] Apparatus

[0093] The processing circuitry may be general purpose processor circuitry configured by program code to perform specified processing functions. The circuitry may also be configured by modification to the processing hardware. The configuration of the circuitry to perform a specified function may be limited exclusively to hardware, limited exclusively to software, or a combination of hardware modification and software execution. Program instructions may be used to configure the logic gates of general purpose or special purpose processor circuitry to perform a processing function. The processing circuitry is described further below.

[0094] The apparatus comprises a control system 200 (figure 8). The control system comprises one or more processors (202) collectively configured to perform the methods for identifying the one or more anomalous monomers in a polymer. That is, to a) receive a stream of input data comprising measurements for a plurality of consecutive monomers in the polymer (204); b) transform the stream of input data in a) into a plurality of vector features; c) calculate the distance between: i) each of the vector features in b); and ii) a plurality of vector features in a control data set, wherein the control data set was generated by transforming one or more streams of control data into one or more pluralities of vector features to form a control data set, wherein the control data comprises measurements for a plurality of consecutive monomers in one or more polymer(s); and d) based on the distances calculated in c), determine if a monomer associated with the one or more vector features in the polymer which is the input data is anomalous relative to the polymer(s) of the control data, and output the determination (206).

[0095] The control system may comprise one or more controllers collectively comprising at least one electronic processor having an electrical input for receiving an input signal; and at least one memory device electrically coupled to the at least one electronic processor and having instructions stored therein; and wherein the at least one electronic processor is configured to access the at least one memory device and execute the instructions. The stream of input data may be stored locally on the apparatus, or on the sequencer (in the system), for example on the hard drive of either the apparatus or the sequencer. Alternatively, input data for the methods, apparatus and system may not be stored locally, e.g. cloud storage.

[0096] System

[0097] The system comprises the apparatus as described above and a nucleic acid sequencing device, for example a nanopore sequencing device. The apparatus is operably connected to the sequencer and configured to receive input data from the sequencer.

[0098] Throughout the specification, unless the context demands otherwise, the terms ‘comprise’ or ‘include’, or variations such as ‘comprises’ or ‘comprising’, ‘includes’ or ‘including’ will be understood to imply the method or kit includes a stated integer or group of integers, but not the exclusion of any other integer or group of integers.

[0099] Each document, reference, patent application or patent cited in this text is expressly incorporated herein in their entirety by reference, which means it should be read and considered by the reader as part of this text. That the document, reference, patent application or patent cited in the text is not repeated in this text is merely for reasons of conciseness. Reference to cited material or information contained in the text should not be understood as a concession that the material or information was part of the common general knowledge or was known in any country.

[0100] Description of the Figures

[0101] Figure 1 shows the oligonucleotide sequences. Figure 2 shows electronic current values for 3 (color-coded) oligonucleotides from the control dataset from position 20 to 40. For each position, the sub-sequence of 5 bases which are in the pore are indicated.

[0102] Figure 3 shows RNA modification detection in synthetic modified oligonucleotides. A Kolmogorov- Smirnov (KS) test is applied, position-wise, to compare the distributions of the anomaly scores in a modified dataset and in a reference dataset devoid of RNA modifications. The p-value (y-axis) is reported for each position (x-axis). The grey shaded areas denote positions that contain a modified base. The green dotted lines indicate the positions where the modified base is at the center of the pore. The p-values highlighted in red correspond to local maxima (smaller peaks are removed until all consecutive peaks are separated by at least nine p-values).

[0103] Figure 4 shows: In the top plot, raw signals from the control data of three different threads (or read names) are depicted. The bottom plot illustrates the signature projection of each thread per position, using UMAP for dimensionality reduction to visualize a high-dimensional vector in a single plot.

[0104] Figure 5 illustrates three examples of the molecule-by-molecule approach in Oligos3 dataset. A curve was fitted to the calibration scores in a specific position, predicting the score for a corresponding test molecule in that position. The resulting p-value (y-axis) is used to interpret modifications in each position (x-axis). Blue bars represent unmodified positions, while red bars indicate positions where a modification is expected. A horizontal red line is added to highlight that no unmodified positions exceed this threshold. It's important to note that missing positions are common at the molecular level.

[0105] Figure 6 shows output scores of three different reads from the same region (positions 1 ,000,776 to 1 ,000,826) on chromosome 20. According to Nanopolish, only the read in 6a) (blue) is methylated. The reads in 6b) and 6c) plots Nanopolish shows no methylation. The grey-shaded area are the positions where the CG bases are inside the nanopore.

[0106] Figure 7 shows a flow chart of an example method. The method 100 comprises, at block 102, obtaining a stream of input data comprising measurements fora plurality of consecutive monomers in the polymer. At block 104, transforming the stream of input data into a plurality of vector features. At block 106, calculating the distance between the vector features from the input data with a plurality of vector features in a control data set. And, at block 108, based on the distance calculated, identifying one or more anomalous monomers. The method may be implemented by the apparatus 200 described herein in reference to Figure 8, but is not limited thereto.

[0107] Figure 8 shows a block diagram of an example apparatus with control system 200.

[0108] Figure 9 shows a further flow chart of an example method. Overview of method

[0109] Prior to the examples, we provide a more detailed description of the method, to exemplify one way of performing the method. Any of the features in points 1-6 below may be combined with the general description above.

[0110] 1. Data Conversion: Sequencing data obtained from the polymer is transformed into an input dataset, comprising multiple vectors. These vectors capture information reflecting the unparametrized shape of the signal around each monomer in the polymer. The transformed vectors are extended to provide concise, informative contextual information about the evolution of the polymer's shape around each monomer. This contextual information is represented in the form of vectors or tensors. This vectorisation step is pivotal to the method as it effectively captures information reflecting for example the shape of each monomer in the polymer.

[0111] 2. Reference Dataset: We create a reference dataset known as the 'corpus of unmodified feature data,' which represents the output range achievable by applying the Data Conversion process to unmodified sequencing data. This reference dataset can be generated through simulation or by converting sequencing data from multiple unmodified polymers.

[0112] 3. Variation Assessment: The data from a single sequence or a collection of sequences is compared to the corpus of unmodified feature data. This comparison can be done at the level of individual vector features, specific locations, or monomers. It entails a comparison of individual vector features from the new sequence with equivalent vectors from the corpus of unmodified feature data. An innovation or conformance score is determined for each location, considering the variability within the reference dataset (also referred to here as the control data set). Calibration of this score is achievable.

[0113] 4. Calculating Conformance Scores: Using a corpus of unmodified feature data in vector or tensor format, we establish an intrinsic conformance scoring function. This function assesses a new instance of a vector feature and provides a score reflecting the degree of innovation compared to the corpus. The corpus is divided into two parts: one for creating the scoring function and the other for forming a set of inliers used for score calibration. Conformance scores should be somewhat invariant to linear rescaling of features and changes in measurement units and should depend on the proximity of the new vector feature to the corpus. For example, we may estimate the nearest neighbour distance from the new instance of the vector feature to the corpus of unmodified feature data using the variance norm defined by the corpus.

[0114] 5. Scoring the Locations: After identifying and calibrating conformance scores for vector features at specific locations or monomers, whether from individual sequence data or aggregated data from multiple sequences, and by considering neighbouring locations, we can flag segments of the sequence where modifications, groups of modifications, or other anomalies have occurred or are likely to occur in the case of multiple converted sequences.

[0115] 6. Modification Determination: The information about locations with a high likelihood of a defect or anomaly, as generated in the previous steps, is then analyzed to determine the nature of the modification, as needed for the specific application. The method identifies modifications based on their novelty relative to the example monomers and sequences in the reference dataset, allowing it to detect complex modifications not seen previously, such as modification pairs, in addition to standard modifications. This step can be combined with contextual information, to improve accuracy.

[0116] A further overview of the process can be found at Figure 9.

[0117] Examples

[0118] Aspects of the present invention will now be illustrated by way of example only and with reference to the following experimentation.

[0119] Example 1 : Preparing the data

[0120] Data

[0121] We used publicly available data that includes seven distinct RNA modifications in synthetic oligonucleotides. However, our method does not directly utilise the nature of these modifications and is applicable to any other types of sequencing data.

[0122] The designed oligonucleotides include different modifications across the same sequence.

[0123] Datasets were made from the oligonucleotides shown in Figure 1 : the first comprises Oligo 1 , which contains two m6A sites at distinct positions and 2,324 different molecules. The second dataset, named Oligo2, include I, m5C, and Pseudoll (^P) modifications with a total of 1 ,289 different molecules. The third dataset, Oligo3, consists of m62A, m1G, and 2’-OMeA modifications and 1 ,982 molecules. Additionally, a fourth dataset Oligo4 was created with an unmodified sequence, serving as a control dataset with 2,205 different molecules.

[0124] For each dataset, namely Oligo 1 , Oligo2, Oligo3 and Oligo4 (control), multiple RNA of these types of oligonucleotides are sequenced by Oxford Nanopore Technologies (ONT) direct sequencing. For each molecule, a univariate electrical signal is measured by a sensor while the oligonucleotide translocates through the nanopore as shown in Figure 2. The DRS electrical signal, recorded discretely in time, is in the form of a univariate time series (or stream). We denote the space of time series over IR by:

[0125] Given a time series x = (x x^) e S(IR) , we call the length of x. A DRS dataset is a collection of such streams, with possibly varying lengths.

[0126] Next, we capture the shape information through using the signature transform from rough path theory to vectorize such time series. We map time series of possibly varying lengths into vectors of fixed dimensionality as described below.

[0127] Transformation

[0128] We transform each initial DRS stream x e S(IR) into a stream x e s(]Rd) with d > 1 , as explained in more detail below. Like this, we add a dimension (or dimensions) to the univariate data. This will allow us to extract more information from the stream before applying the signature transform by mapping one stream of data to another, aiming to extract relevant information for the specific problem at hand. These transformations and their combinations can be used for optimization purposes.

[0129] The transformations, offer various advantages. They can help reduce the dimension of the time series, preprocess the time series before applying the signature map to facilitate information extraction, exploit the signature's invariance to translation and / or reparametrization, and more as elaborated further in [1 , 2], While there are numerous pre-signature operations, we will provide some examples. As before, let's assume that we observe x. For any vector with increasing timestamps t e S(]R), we have a time augmentation by: which basically adds the timestamps as an extra coordinate, with a length of n. The lead-lag transformation of a d-dimensional stream of data, results in a 2t / -dimensional stream of data with a length of 2n - 1 , defined as follows:

[0130] The invisibility-reset transformation adds a dimension as: Vectorization

[0131] Rough path theory is a field in stochastic analysis that offers a rigorous mathematical framework for modelling the interaction between a stream and a physical control system. The signature transform, in particular, stands out for its remarkable properties, making it an effective feature mapping technique capturing the unparametrized shape of streamed data. For any continuous path x [0, T] -> ]Rdof bounded variation, and any s, t e [0, T], the signature of | [sis defined by that is, an infinite collection of (Riemann-Stieltjes) iterated integrals. The k01iterated integral can be viewed as a collection of dkscalars. We compute the first (say M e N with M > 1) iterated integrals of the signature of x |[s t], the piecewise linear interpolation of a stream x, with a dedicated highly optimized Python library. We obtain a vector of < / -i scalars. We note that computing the signature of a piecewise linear path with pieces has time complexity O(dM-e).

[0132] We choose the signature due to its ability to capture the sequential order of events and effectively model its effects without the need for high-dimensional recovery of individual data points. The components of the signature form our feature set, allowing for subsequent operations.

[0133] After illustrating the raw signal in Figure 2, Figure 4 shows the appearance of vectors after the signature transformation for some raw signals. It's essential to note that the signature generates a high-dimensional vector, requiring dimension reduction for plotting. We used UMAP https: / / umap- learn.readthedocs.io / en / latest / 1 allowing us to learn a manifold embedding that captures the signature features observed in the data. Figure 4 shows the raw clean data plots for three molecules along with the corresponding projection representation of the signature at the bottom.

[0134] Example 2: Detecting modifications in polynucleotide data

[0135] In the previous section, we have explained how to obtain a convenient vectorial representation of the DRS signals. We explain how the time intervals have been chosen for the experimental study.

[0136] Novelty score

[0137] Now, we explain how we assign a conformance score by computing the nearest neighbour Mahalanobis distance of the vectorized DRS signal y = S(xs t) around a monomer location, using a corpus of unmodified feature data around the same location.

[0138] We have access to unmodified data, like the control data shown in Figure 2 (also referred to as the “reference data set” above). We refer to this dataset as a corpus and denote it by Dn={y1, •••>yn}- Given a corpus, we compute two types of statistics, namely the sample mean and the sample covariance matrix

[0139] The Mahalanobis distance between an input y and an element ytof the corpus, is then defined as

[0140] Given an input y , we compute its Mahalanobis distance to every element of the corpus and use this collection of distances to construct a score a(y ; Z>n). Given that the Mahalanobis distance incorporates the covariance matrix in its definition, it serves as an indicator to normalise distances based on the variability to the corpus. The conformance score of y is defined as the nearest- neighbor (NN) distance, namely

[0141] We thereby have a data-driven notion of distance.

[0142] Determination of whether the difference represents an anomaly

[0143] Now, we explain how we build a one-class classifier for anomaly detection, in other words, how given a score we decide whether the signal is associated with a location that has been modified or not. To this aim, we use a portion of the modification-free data that was not used in the corpus, which we refer to as the calibration set and denote (similarly to the corpus) by

[0144] We precompute the scores of each element of the calibration set a(yi; Z>n), ..., a(yl'n; T)n). In this way we define what is clean in our setup. The same method is applied to the test set which might contain modifications, and the scores are computed for each element of the test set

[0145] Aggregated Result

[0146] We then analyse the calibration and test score distributions by conducting a statistical test to determine the difference between them, which is guantified by obtaining their respective p-values. This process involves comparing the two distributions, enabling us to assess their similarity through a rigorous statistical evaluation based on their distinct characteristics.

[0147] Given the precomputed set of calibration scores T>cal= {a(y'; Z>n) y'| e !),„} and the set of test scores T>test= {a(y*; Z>„) y*| e 2)*}, to compare their distributions, we perform a statistical test. Let P(Dcal) and P(Dtest) represent the probability distributions of the calibration and test datasets, respectively. To assess the difference between these two distributions we applied a Kolmogorov- Smirnov two-sample test. The test evaluates the null hypothesis that the two distributions are the same. The obtained p-value guantifies the level of similarity between the two distributions.

[0148] Molecule-by-molecule

[0149] We also possess the capability to detect modifications at the individual molecule level (based on each individual score a(y*; 2)„) ). Unlike the aggregated results that involve comparing two distributions, our molecule-by-molecule approach entails fitting a curve to the calibration scores (T>cal) distribution in a specific position. Subseguently, we make a prediction on one molecule test score in the same position based on this fit, resulting in a p-value. This process involves comparing the calibration score distribution to an individual test score, providing an understanding of the probability of the calibration score being greater than the test score in a single molecule.

[0150] Results:

[0151] The aggregated results are shown in Figure 3.

[0152] The molecule-by-molecule results are shown in Figure 5.

[0153] The modification detection results in the DRS signals are shown in Figure 3. Our scoring function, dependent on control data, is applied to every position and read within the dataset. Then, for each position, we compare the dataset scores distribution (Oligol , Oligo2, or Oligo3) to the calibration scores distribution.

[0154] The molecule-by-molecule results for DRS signals are depicted in Figure 5, showing three different molecules from the Oligos3 dataset. While this serves as an illustrative example of how we can sometimes correctly identify modifications at the molecular level, it's essential to note that this may not always be the case, given the substantial number of molecules — more than 2,000 per dataset.

[0155] This highlights the robustness of our aggregated results too.

[0156] By applying the scoring function per position we mean the following. One postprocessing step commonly applied to the ONT DRS signal consists in dividing the raw data series x = (x-i, . into 5-mer-corresponding signal segments. The resulting segmentation is depicted by the vertical dashed lines on Figure 2. We call a window in this segmentation, a position. When the current is recorded over that window / position, approximately 5 bases reside in the pore. For this experimental study, we compute the signature and score the signals over such windows.

[0157] The highest p-values consistently align with positions that contain at least one modified base. Notably, distinct peaks spread across the sequence, aligning with the initially modified positions. It is important to note that in the aggregated results no peaks are observed outside these specific positions, indicating an improvement over other methods which require exact knowledge of the types of modifications you are looking for. This shows the sensitivity of our scoring function to modifications.

[0158] By conducting a statistical test, we can identify and highlight regions where an anomaly has occurred, flagging them as potentially modified. This could lead to a more detailed analysis to determine the specific type of modification. At this point, we are not categorizing the anomalies according to a modification, just determining their location.

[0159] We present a semi-supervised approach that does not require prior knowledge of specific modifications and accurately pinpoints their locations. As a result, this method holds the potential to identify novel or previously unidentified modifications.

[0160] Example 3: The method works to identify anomalous monomers in various polymers

[0161] To demonstrate the technique works on various types of polymers, we used the tool Nanopolish to identify

[0162] We downloaded publicly available DNA data processed using Oxford Nanopore sequencing4. Specifically, we obtained the first dataset, FAB39088. The data was processed with the tool Nanopolish5, which includes a module called eventalign, that is widely used to align events or "squiggles" to a reference genome. Nanopore sequencing is sensitive to base modifications, and Nanopolish provides a step-by-step tutorial for detecting one such modification: DNA methylation. Methylation is relatively straightforward to detect because it typically occurs in a CpG context (kmers containing a C base followed by a G).

[0163] Using Nanopolish, we obtained labels for chromosome 20, covering positions 1 ,000,000 to 60,000,000, indicating for each read name and location a log-likelihood ratio signifying Nanopolish belief as to the positions in the read name where the nucleotide had been methylated. Nanopolish identifies modifications instead using a hidden Markov model and requires the modification to be known, as opposed to the use of signature based transformation as used in our method detailed in examples 1 and 2.

[0164] To validate our method, we identified a region where Nanopolish reported some read names as having a methylated modification and other read names as clean. We used the region 1 ,000,776 to 1 ,000,826 (shown in Figure 6). Three read names were called over this region and methylation detected in only one. We compiled a list of the kmers occurring in this region.

[0165] For each of these kmers, we considered all read names and matching clean kmers in the interval 1 ,000,000 to 60,000,000 that were not located in the short region.

[0166] We randomly selected 3,000 location-read name pairs from these clean positions. We used the associated Nanopore data as our corpus of normality for the given kmer. Note that these kmer read names pairs were sampled from diverse positions outside the interval being tested.

[0167] For locations in the identified region and for our three read names, we calculated the conformance of the Nanopore data at each position against the corpus of clean data associated for the kmer at that position.

[0168] Results:

[0169] The unique modified location (k-mer in the individual read-line) had a poor conformance, all other location in the region as well as the potentially modified location in the two clean read names showed conformance. This demonstrates the broad ability, beyond the RNA context, of the methodology to identify modifications against the corpus.

[0170] The shaded area in Figure 6 corresponds to the region where a CG base was inside the Nanopore, indicating potential methylation. We observe a high output score in this shaded area for the methylated read name identified by Nanopolish, and low output scores across the region in the absence of methylation.

[0171] This approach confirms that our method can localize methylation patterns in DNA data.

[0172] References

[0173] [1] James Morrill, Adeline Fermanian, Patrick Kidger, and Terry Lyons. A generalised signature method for time series. arXiv preprint arXiv:2006.00873, 2020 [2] Thomas Cochrane, Peter Foster, Terry Lyons, Imanol Perez Arribas. Anomaly detection on streamed data. arXiv preprint arXiv:2006.03487, 2020

[0174] [3] Laxhammar, Rikard. Conformal anomaly detection: Detecting abnormal trajectories in surveillance applications. Diss. University of Skovde, 2014. [4] https: / / github.com / nanopore-wgs-consortium / NA12878 / blob / master / Genome.md

[0175] [5] https: / / github.com / jts / nanopolish / tree / master

Claims

CLAIMS1. A computer-implemented method for identifying one or more anomalous monomers in a polymer, the method comprising: a) obtaining a stream of input data comprising measurements for a plurality of consecutive monomers in the polymer; b) transforming the stream of input data in a) into a plurality of vector features; c) calculating a distance between: i) each of the vector features in b); and ii) a plurality of vector features in a control data set, wherein the control data set was generated by transforming one or more streams of control data into one or more pluralities of vector features to form a control data set, wherein the control data comprises measurements for a plurality of consecutive monomers in one or more polymer(s); and d) based on the distances calculated in c), determining if a monomer associated with the one or more vector features in the polymer which is the input data is anomalous relative to the polymer(s) of the control data.

2. The method of claim 1 , wherein the measurements for the control data are: synthetic and / or non-synthetic measurements.

3. The computer-implemented method of any of claims 1-2, wherein the polymer is a biopolymer, optionally wherein the biopolymer is a polynucleotide, optionally wherein the polynucleotide is DNA or RNA.

4. The computer-implemented method of any of claims 1-3, wherein the method identifies polymers comprising monomers with more than one type of anomaly, optionally wherein the anomaly is a covalent modification.

5. The computer-implemented method of claim 3, wherein the polynucleotide is RNA and the anomaly is a covalent modification.

6. The computer-implemented method of any of claims 1-5 wherein the one or more polymer(s) without the anomaly used to generate the control data set comprise: a) polymers consisting of unmodified monomer; orb) polymers comprising one or more monomers comprising any one or more of the following: 5-methylcytosine, 5-hydroxymethyl-cytosine, Inosine, pseudouridine, N6- dimethyladenosine, 1 -methylguanosine, 2’-0 methyladenosine, 7-methylguanosine or N6- methyladenosine.

7. The computer-implemented method of any of claims 1-6 wherein the stream of data is measured from the polymer during translocation of the polymer through a nanopore.

8. The computer-implemented method of any of claims 1-7, wherein the measurements comprise an electrical signal or an optical property.

9. The computer-implemented method of claim 8, wherein the input data is an ionic current.

10. The computer-implemented method of any of claims 1-9, wherein the vector features obtained from the input data are normalised to the same scale as vector features of the control data set.

11. The computer-implemented method of claims 1-10, wherein the determination in step d) comprises checking if the distances calculated are beyond thresholds for anomaly, wherein the thresholds for anomaly were determined by: i) obtaining calibration data which is representative of the control data; ii) ensuring the distance function in step c) is independent of the calibration data; iii) transforming the one or more stream(s) of calibration data into one or more pluralities of vector features; iv) calculating the distance between each of the vector features in iii) and the plurality of vector features in the control data set; v) using the distances calculated in iv) to set thresholds for anomaly for the input data; and vi) applying a threshold obtained from v) to the input data, optionally wherein the calibration data was extracted from the control data.

12. The computer-implemented method of any of claim 11 , further comprising aggregating the threshold determinations associated with consecutive monomers to form new threshold determinations.

13. The computer-implemented method of any of claims 1-12, wherein vectorisation in b) for the input data is consistent with vectorisation of the control data in c).

14. The computer-implemented method of any of the previous claims, wherein steps b) and c) comprise extracting vector features from the signature based transformation of the input data and / or control data respectively.

15. A computer-implemented method for identifying one or more anomalous monomers in a polymer, the method comprising: a) obtaining a stream of input data comprising measurements for a plurality of consecutive monomers in the polymer, wherein the stream of input data is an ionic current measured from the polymer during translocation of the polymer through a nanopore; b) transforming the stream of input data in a) into a plurality of vector features; c) calculating a distance between: i) each of the vector features in b); and ii) a plurality of vector features in a control data set, wherein the control data set was generated by transforming one or more streams of control data into one or more pluralities of vector features to form a control data set, wherein the control data comprises measurements for a plurality of consecutive monomers in one or more polymer(s), and wherein steps b) and c) comprise extracting vector features from the signature based transformation of the input data and / or control data respectively; and d) based on the distances calculated in c), determining if a monomer associated with the one or more vector features in the polymer which is the input data is anomalous relative to the polymer(s) of the control data, wherein the polymer is RNA or DNA.

16. The computer-implemented method of any of the preceding claims, further comprising out-putting from the method data indicating the position of one or more anomalous monomer(s).

17. A computer-readable storage medium or a computer program comprising computer-executable instructions, which when executed by a computing system, are capable of causing the computing system to perform the method according to any of the preceding claims.

18. A computer-implemented method of preparing a control data set for use in a method of identifying one or more anomalous monomers in a polymer, the method comprising: a) obtaining one or more stream(s) of control data comprising measurements for a plurality of consecutive monomers in one or more polymer(s); and b) transforming the one or more stream(s) of control data into one or more pluralities of vector features to form a control data set, optionally wherein the control data set is normalised.

19. An apparatus comprising processing circuitry configured to perform the method of any one of claims 1-16.

20. A system comprising: a) a nucleic acid sequencing apparatus; and b) the apparatus of claim 19, wherein the apparatus is operably connected to the sequencing apparatus and configured to receive a stream of input data from the sequencing apparatus.

Images