算法帮助!快速搜索字符串及其伴侣的算法

8
我正在寻找一种在巨大字符串(由数亿到数十亿个字符组成的生物基因组序列)中进行快速搜索的算法。
该字符串仅包含4个字符 {A,C,G,T},其中"A"只能与"T"配对,而"C"则与"G"配对。
现在我正在搜索可以反向配对的两个子字符串(两个子字符串的长度限制在{minLen,maxLen}之间,区间长度在{intervalMinLen,intervalMaxLen}之间)。
例如,字符串为:ATCAG GACCA TACGC CTGAT
约束条件:minLen = 4,maxLen = 5,intervalMinLen = 9,intervalMaxLen = 10
结果应该是:
1. "ATCAG"与"CTGAT"配对 2. "TCAG"与"CTGA"配对
谢谢!更新:我已经有了确定两个字符串是否可以配对的方法。唯一的问题是进行详尽搜索非常耗时。

4
听起来需要使用正则表达式 (System.Text.RegularExpressions.Regex)。这不是作业吧? - Mark Avenius
minLen、maxLen、intervalMinLen 和 intervalMaxLen 的典型值是什么? - ElKamina
intervalMinLen和intervalMaxLen控制什么? - Russell Borogove
@Mark:“你耳边的铃声是 Jamie Zawinski 说的‘有些人面对问题时会想:‘我知道了,我要使用正则表达式。’现在他们有两个问题了。’” - Tim Schmelter
1
你也可以在 http://biostar.stackexchange.com 上提问。 - Pierre
显示剩余5条评论
5个回答

4
我知道你不是在搜索子字符串,但我认为创建一个包含匹配的反转基因组字符串可能是值得的;任务将是在两个字符串中找到公共子串。

例子:

原始字符串

  ATCAG GACCA TACGC CTGAT

反转字符串:

  TAGTC CGCAT ACCAG GACTA

如果您将该字符串转换为其配对值(将T<->A和C<->G替换),则会得到有用的内容:
  ATCAG GCGTA TGGTC CTGAT

我知道这个预处理过程会很耗费时间和空间,但是之后你就可以使用标准的字符串算法,并且考虑到你要搜索的比较数量,这样做肯定是值得的。
当原始字符串和反向查找字符串相同时,我认为你的问题听起来与“最长公共子串”问题非常相似,该问题有很好的描述。你的第二次预处理将是构建后缀树以允许快速查找子字符串。
你最终会得到二次运行时间,但我怀疑你无法做得更好。

Edit 介绍了后缀树。 - faester

3

3
最简单的解决方案是从数据中构建一棵 Trie,最大高度为 maxLen。每个节点还应该有一个列表,其中包含找到特定字符串的位置(如果按顺序创建 Trie,则该列表将按升序排列)。
然后,在搜索时,只需反转补码搜索字符串并遍历 Trie。当找到匹配项时,请检查是否有一个匹配项位于正确的区间内,如果“是”,则输出字符串。
如果需要详细算法,请告诉我。

“夸奖字符串并遍历字典树”。从给定的示例来看,似乎两个字符串不需要按顺序匹配,而只需要具有相同数量的每个补充字符(AC可以匹配GT)。仅仅遍历字典树是无法实现这一点的,您认为呢? - efficiencyIsBliss
@efficiencyIsBliss 感谢您指出这一点。它应该“反向赞美字符串并遍历 trie”。 - ElKamina

2

我认为这是一个有趣的问题,因此我编写了一个基于“折叠”的程序,它从不同的“折叠点”向外扫描可能的对称匹配。如果N是核苷酸数量,M是“maxInterval-minInterval”,则运行时间应该为O(N*M)。我可能会错过一些边界情况,所以请小心使用代码,但它确实适用于提供的示例。请注意,我使用填充的中间缓冲区存储基因组,因为这减少了内部循环中所需的边界情况比较次数;这通过额外的内存分配来换取更好的速度。如果您进行任何更正或改进,请随意编辑帖子。

class Program
{
    public sealed class Pairing
    {
        public int Index { get; private set; }

        public int Length { get; private set; }

        public int Offset { get; private set; }

        public Pairing(int index, int length, int offset)
        {
            Index = index;
            Length = length;
            Offset = offset;
        }
    }

    public static IEnumerable<Pairing> FindPairings(string genome, int minLen, int maxLen, int intervalMinLen, int intervalMaxLen)
    {
        int n = genome.Length;
        var padding = new string((char)0, maxLen);
        var padded = string.Concat(padding, genome, padding);

        int start = (intervalMinLen + minLen)/2 + maxLen;
        int end = n - (intervalMinLen + minLen)/2 + maxLen;

        //Consider 'fold locations' along the genome
        for (int i=start; i<end; i++)
        {
            //Consider 'odd' folding (centered on index) about index i
            int k = (intervalMinLen+2)/2;
            int maxK = (intervalMaxLen + 2)/2;
            while (k<=maxK)
            {
                int matchLength = 0;
                while (IsPaired(padded[i - k], padded[i + k]) && (k <= (maxK+maxLen)))
                {
                    matchLength++;

                    if (matchLength >= minLen && matchLength <= maxLen)
                    {
                        yield return new Pairing(i-k - maxLen, matchLength, 2*k - (matchLength-1));
                    }
                    k++;
                }
                k++;
            }

            //Consider 'even' folding (centered before index) about index i
            k = (intervalMinLen+1)/2;
            while (k <= maxK)
            {
                int matchLength = 0;
                while (IsPaired(padded[i - (k+1)], padded[i + k]) && (k<=maxK+maxLen))
                {
                    matchLength++;

                    if (matchLength >= minLen && matchLength <= maxLen)
                    {
                        yield return new Pairing(i - (k+1) - maxLen, matchLength, 2*k + 1 - (matchLength-1));
                    }
                    k++;
                }
                k++;
            }
        }
    }

    private const int SumAT = 'A' + 'T';
    private const int SumGC = 'G' + 'C';
    private static bool IsPaired(char a, char b)
    {
        return (a + b) == SumAT || (a + b) == SumGC;
    }


    static void Main(string[] args)
    {
        string genome = "ATCAGGACCATACGCCTGAT";
        foreach (var pairing in FindPairings(genome, 4, 5, 9, 10))
        {
            Console.WriteLine("'{0}' pair with '{1}'",
                              genome.Substring(pairing.Index, pairing.Length),
                              genome.Substring(pairing.Index + pairing.Offset, pairing.Length));
        }
        Console.ReadKey();
    }


}

这是一个类似于折叠的问题。谢谢。 - Mavershang

1

我会考虑将字符串转换为16位长度的二进制:

A = 101
T = 010
C = 110
G = 001

这样每个16位单元最多可以容纳5个字符。与字符串搜索和比较相比,这应该非常快。

使用最高位来确定它是4序列单元还是5序列单元。

现在,您可以快速排序并将所有4序列和5序列单元分成自己的组(除非您需要找到部分匹配5序列单元的4序列单元)。

对于比较,您可以再次进行排序,然后您会发现以G开头的所有序列都会出现在以T、A和C开头的序列之前。然后,您可以取一个序列并将其与仅可能的排序数组部分进行比较-这应该是一个更短的问题空间。

此外,我选择这些三位编码的原因是您可以简单地将两个字符串异或在一起以查看它们是否匹配。如果连续得到15个1,则两者完全匹配。否则,它们不匹配。

你需要想办法解决反向并联的问题 - 我有一些想法,但这应该能让你开始,如果反向并联部分成为问题,请提出另一个问题。


这些内容可能有用,也可能没有用,因为我确定我错过了一些问题和基因测序的方面,但编码和比较至少应该使今天处理器上的比较非常快速。 - Adam Davis
1
为什么不使用A=00;T=01;C=10;G=11? - phoog
1
为什么不只用2位二进制(A:00,T:11,G:01,C:10)? - ElKamina
两个位的问题在于你不能将两个序列异或在一起并得到一个集合结果。请注意我的编码是对称的。这允许使用简单的逻辑操作来推断两个序列是否匹配 A-T 和 G-C。 - Adam Davis
我的也是。A.xor.T = 11,G.xor.C = 11。 - ElKamina
2
使用三个位的原因主要是因为OP遇到的问题之一是需要反转一组序列(反平行)。3个位可以清楚地表示所代表的字符是正向还是反向,这取决于所需的排序和比较操作。如果只有两个位,您无法确定它是C还是反向的G。 - Adam Davis

网页内容由stack overflow 提供, 点击上面的
可以查看英文原文,
原文链接