A computer-implemented and reference-free method for identifying variants in nucleic acid sequences

Abstract

There is provided a computer-implemented method for identifying of nucleic acid variants between two cells, such as a normal cell vs. a pathological cell of a patient, or a cell at two different stages of development. The method is alignment-free, as it does not depend on the use of a reference genome, and is based on the generation and comparison of polymorphic k-mers derived from the nucleotide sequence reads of both biological states. The invention accurately identifies all sorts of genetic variants, ranging from single nucleotide substitutions (SNVs) to large structural variants with great sensitivity and specificity. As a major novelty, it also identifies non-human insertions, such as those derived from retroviruses. Altogether, this invention allows the integration with specific hardware architectures in order to speed up the executions to an unprecedented level.

1 . A computer-implemented method for identifying of nucleic acid variants between two genomic states comprising the steps of:
 A) Inputting 2 sets of nucleic acid reads, which are sequences retrieved from a nucleotide sequencing method, wherein the first set of reads corresponds to cells representing a first test state, and the second set of reads corresponds to cells representing a second control state;   B) Filtering the reads, wherein the filtering comprises:
 B1) Keeping only the reads with at least a percentage X1 of their bases with a Phred quality score higher than 20, being X1 equal to or above 90%; 
 B2) Splitting the reads with an undefined nucleotide, giving one sequence before, and one sequence after the undefined nucleotide, the latter being discarded; and 
 B3) Discarding the sequence reads with less than X2 bases, wherein X2 is from 25 to 50; 
   C) Generating a hashtable structure comprising:
 C1) Generating a number of N-X2+1 new reads for each read of sequence length N, wherein the new N-X2+1 reads correspond to all k-mers with length X2 nucleotides; and 
 C2) Building a hashtable structure, which comprises all the k-mers generated in step C1) and further comprises the number of times each k-mer is observed in the two sets of reads corresponding to first and second states. 
   D) Detecting candidate variants in the sequence between first state and second state, wherein a k-mer of the hashtable structure is taken as a candidate breakpoint, which represents a variant between the first and second states, if it fulfills all the following requirements:
 D1) At least one inflection based on a k-mer's stem must have at least X3 reads with the same variation between first and second states, being X3 at least 2; 
 D2) The percentage of first state reads in second state reads is not over a threshold X4, to account for possible contamination of control state reads with test state reads, wherein X4 is at least 5%. 
   E) Clustering and filtering test and control reads derived from all candidate breakpoints accepted in step D to build blocks, by carrying out the steps:
 E1) Retrieving reads which contain the stem of at least one k-mer that represents the candidate breakpoint selected in step D); 
 E2) The reads of step E1) with at least X5 k-mer variants within a window of X6 nucleotides are taken as leading reads, wherein X5 is at least 7 and X6 is at least 10; 
 E3) Reads whose k-mers share at least one stem with a leading read are merged to give a block; and 
 E4) If the nucleic acid whose variant is being identified is a double stranded DNA, then both forward and reverse variants are taken into account when building the block. 
   F) Aligning blocks taking their leading reads as a reference:   F1) For each read in the block, take the leading read's stem and find the longest inflection or partial inflection between the read and the leading read. F2) Successively position each read so that its matching inflection or partial inflection is aligned against the leading read;   wherein:   “candidate breakpoint” means the candidate nucleotide position where the sequence changes, that is immediately flanking a sequence variant.   “inflection” means a fragment of length k that can be derived from extending a stem of a length lower than the length of a k-mer of length k; and   “stem” means a fragment of a k-mer of length k with S defined bases, where S<k, and k-S omitted (undefined) bases; a k-mer being understood as all possible substrings of length k that are contained in a string.   
     
     
         2 . The computer-implemented method according to  claim 1 , further comprising the step following F2:
 F3) Obtaining first state scores, second state scores, and global similarity scores of each position in the block by measuring a ratio of most frequent nucleotide in that position relative to the total number of nucleotides.   
     
     
         3 . The computer-implemented method according to  claim 1 , further comprising the step of:
 G) Cataloguing and annotating blocks according to the following:
 G1) If blocks between the first and second states only differ in one substituted nucleotide, the variant is catalogued as containing a single nucleotide variant and the single nucleotide variant is annotated; 
 G2) If blocks between the first and second states differ in more than one nucleotide but the whole difference in sequence is contained within the block, the variant is catalogued as a small structural variant, and the small structural variant is annotated; and 
 G3) If blocks between the first and second states differ in more than one nucleotide and the whole difference in sequence is not contained within the block, the variant is catalogued as a large structural variant, and the boundaries of all large structural variants are extended by retrieving blocks overlapping at least X2 nucleotides in an iterative process which ends when the extended sequence reaches 200 nucleotides or when an ambiguous path is found. 
   
     
     
         4 . The computer-implemented method according to  claim 1 , further comprising the step of:
 H) Filtering of the blocks, according to the following:
 H1) The percentage of second state reads in first state reads is not over a threshold X7, to account for possible contamination of test state reads with control state reads, wherein X7 is at least 20%; 
   
     
     
         5 . The computer-implemented method according to  claim 1  further comprising optionally mapping second state blocks, and subsequently mapping first state blocks, on a reference genome. 
     
     
         6 . The computer-implemented method according to  claim 1 , wherein X2 is from 25 to 40. 
     
     
         7 . The computer-implemented method according to  claim 1 , wherein X3 is equal to or above 3. 
     
     
         8 . The computer-implemented method according to  claim 1 , wherein threshold X4 is 5%. 
     
     
         9 . The computer-implemented method according to  claim 1 , wherein the threshold X5 is from 10 to 15 and the threshold X6 is from 12 to 20. 
     
     
         10 . The computer-implemented method according to  claim 1 , wherein the first set of reads corresponds to pathological cells of a patient and the second set of reads corresponds to non-pathological cells of the same patient. 
     
     
         11 . A computer program product comprising program instructions for causing a computer system to perform the method for identifying nucleic acid variants between two genomic states as defined in  claim 1 . 
     
     
         12 . The computer program product according to  claim 11  embodied on a storage medium. 
     
     
         13 . The computer program product according to  claim 11  carried on a carrier signal. 
     
     
         14 . (canceled) 
     
     
         15 . A computer system comprising a processor and a memory, wherein the memory stores computer executable instructions that, when executed by the processor, cause the system to perform the method for identifying nucleic acid variants between two genomic states as defined by  claim 1 . 
     
     
         16 . The computer-implemented method according to  claim 1 , further comprising the steps of:
 Following step F2, F3) Obtaining first state scores, second state scores, and global similarity scores of each position in the block by measuring a ratio of most frequent nucleotide in that position relative to the total number of nucleotides.   G) Cataloguing and annotating blocks according to the following:
 G1) If blocks between the first and second states only differ in one substituted nucleotide, the variant is catalogued as containing a single nucleotide variant and the single nucleotide variant is annotated; 
 G2) If blocks between the first and second states differ in more than one nucleotide but the whole difference in sequence is contained within the block, the variant is catalogued as a small structural variant, and the small structural variant is annotated; 
 G3) If blocks between the first and second states differ in more than one nucleotide and the whole difference in sequence is not contained within the block, the variant is catalogued as a large structural variant, and the boundaries of all large structural variants are extended by retrieving blocks overlapping at least X2 nucleotides in an iterative process which ends when the extended sequence reaches 200 nucleotides or when an ambiguous path is found; 
   H) Filtering of the blocks, according to the following:
 H1) The percentage of second state reads in first state reads is not over a threshold X7, to account for possible contamination of test state reads with control state reads, wherein X7 is at least 20%; and 
   optionally mapping second state blocks, and subsequently mapping first state blocks, on a reference genome.   
     
     
         17 . The computer-implemented method according to  claim 16  further comprising optionally mapping second state blocks, and subsequently mapping first state blocks, on a reference genome. 
     
     
         18 . The computer-implemented method according to  claim 1 , wherein threshold X2 is from 25 to 40, X3 is equal to or above 3, threshold X4 is 5%, threshold X5 is from 10 to 15 and the threshold X6 is from 12 to 20. 
     
     
         19 . The computer-implemented method according to  claim 16 , wherein threshold X2 is from 25 to 40, X3 is equal to or above 3, threshold X4 is 5%, threshold X5 is from 10 to 15 and the threshold X6 is from 12 to 20. 
     
     
         20 . The computer-implemented method according to  claim 16 , wherein the first set of reads corresponds to pathological cells of a patient and the second set of reads corresponds to non-pathological cells of the same patient.