I’m seeking input from anyone on the Wolfram community forums who cares to read about an algorithm for finding “approximate” subsequences in high-quality long-read chromosome DNA assemblies. At present I am seeking improvements to total run time, but other advice is welcome.INTRODUCTIONThe data sets. • Specimen chromosomes. Of interest here are perennial plants, which fall into the higher genomic classification of eukaryotes. For quality control, these are chiefly from PacBio Sequel devices and have 30x or more Genome coverage (the assayed sequences overlap each chromosome 30 times or more). Eukaryote chromosomes can be diploid or higher polyploids. There are speculative assay and analysis tools for higher ploidy DNA but the present application is limited to diploid specimens. The “letters” in a DNA sequence are called “base pairs”, and the number of them in a sequence has units of “bp”. Chromosome sizes in perennial plants range from about 10 Mbp to over 60 Mbp. The number of chromosomes in a species is constant and ranges from a few to over 30. Example data: https://www.ncbi.nlm.nih.gov/datasets/genome/GCA_009761775.1/ (15 chromosomes, 14.8 Mbp to 30.6 Mbp in length). • PCR molecular markers. These come in sets of {forward, reverse} pairs called primers. The sets are often specific to a species, but some have utility at higher order (Genus, Family). The primer pairs are labeled {F, R}. In published data sets they are typically in what is termed 5’-3’ order. Specimen chromosome sequences often have the ordering matching F 5’-3’. In this case, the R primer must be placed in 3’-5’ order (with BioSequenceReverseComplement) to match its occurrence in the specimen chromosome. However, some markers have been obtained in such a way that the R is actually the F in the specimen chromosome, and the reverse complement of the F is the actual R. This situation can be verified by preliminary exact searches for the primers in the chromosome (of which only 1 or a few might be found). In the algorithm described here, the {F, R} primers are searched for individually. PCR primer lengths are often in the range 16 bp to 24 bp. The marker set sizes are most commonly about a dozen, but overall range from a few to a few thousand. Example data relevant to the chromosomes above: https://doi.org/10.1080/14620316.2005.11511920 p.219 (SSR markers), https://doi.org/10.3389/fpls.2021.755358 Table 1 (SNP markers). Note that the latter paper is based on a prior version of the chromosome data set listed above.Nature of the data. • Actual primer occurrences. PCR assay has revealed that actual occurrences range from: exact, exact length but a few altered bp, expanded length with a few inserted bp and possibly a few more incorrect bp, shorter than primer length with a few deleted bp and possibly a few more incorrect bp. In all cases, the correct bp present in the assay are in the sequence order of the primer. For example, marker AG-Marker-F might occur as: CAGGAGGAGGAACGCT (exact), CAGGAAGAAGAACGCT (same length, two G’s appear as A’s), CAGGAAAGGAGGAACGCT (two inserted A’s but otherwise exact), GGAGAGGAGCGCT (missing 1st two bp and 7th bp, and the A five positions from the end is a G). • The reality of genome assemblies. Two types of errors occur in chromosome level genome assemblies: local and regional (my terms). Local error refers to one or two incorrect bp in a short sequence. For example, AGTT is represented as AGGT. This is rare in high-coverage long-read sequencing due to the tremendous amount of overlapping reads. Regional errors are due to misassembly. This occurs when a significant length segment is misplaced or is not actually part of the chromosome (See https://doi.org/10.1111/tpj.14635 for more details). Regional errors cannot be accounted for by the algorithm under discussion here.The search criteria. • Every search returns a “candidate”. • Letters in a candidate are only valid when they appear in primer order. Missing primer letters are ok. • Two parameters are used to specify searches: and wA. is the maximum permissible incorrect bp in a candidate. wA is the “width allowance” in the candidate search. • I currently use = Max[4, Round[20% of primer length]]. This is based on observations of PCR assays of primers. • For any given there will be 2· + 1 width allowance values ranging from - to + (including zero). • Negative wA values specify that the candidate length is shorter than the primer. There are multiple ways this can occur: at the beginning of the candidate, in the middle, at the end, or some combination of the three. • Specifying the search criteria by a single RegularExpression and attempting to use it with StringPosition on an entire chromosome will quickly exceed the recursion limit of RegularExpression. • At present I multiplex the search criteria into many fixed-length regular expressions. For = 4, there are about 70 regex, for = 5, there are about 200 regex.
ϵ
accept
ϵ
accept
ϵ
accept
ϵ
accept
ϵ
accept
ϵ
accept
ϵ
accept
ϵ
accept
ϵ
accept
CURRENT ALGORITHM FOR SINGLE CPU COMPUTER
kernelCount=9;(*dependsonnumberofCPUcores,installedRAM,andWolframlicense*)Do[ chrSequence=Import[chrFileSpec]; Do[ regexSpecifications=multiplexer[primer]; (*typically~70,occasionally200,rarelyafew1000*) InitKernels[kernelCount]; ParallelDo[ regexString=Spec2String[regexSpec]; positions=StringPosition[chrSequence,RegularExpression[regexString]]; (*1000’swillbereturned*) vettedCandidates=OrderedPartsCalc[positions,regexSpec,chrSequence]; exportFile=ExportFileCalc[regexSpec]; Export[exportFile,vettedCandidates], {regexSpec,regexSpecifications},Method”CoarsestGrained” ]; CloseKernels[]; regexSpecifications=Null; Pause[1], {primer,primers}(*afewtoafew1000*) ]; chrSequence=Null, {chrFileSpec,chromosomeFiles}(*typicallyaboutadozen*)];