r - 在具有多行的 data.frame 中识别部分匹配的字符串(DNA 序列)所需的解决方案

Question

我正在寻找以下问题的解决方案:

我确实有一个超过 600 万行的数据框，其中一行包含测序信息(DNA 序列)。根据报告数据集的方式，数据框中会有重复的行。但是:这些重复并不是完美匹配。让我用一个例子来说明这一点。

row 1: ATCTCAGCATCATACCAACTACTA
...
row 5: ATCTCAGCATCATA..........

前一个 block 在数据框的两个不同行中显示了两个序列。这些点只是为了可视化目的而显示(它们不是数据集的一部分)。

目标是:将这些序列标记为相同。 (最后，我的目标是为每一行分配一个序列 ID，因此这两行应该具有相同的序列 ID，因为第 5 行中的序列是第 1 行序列的一部分，因此序列可能相同。

我尝试使用 base R 的 match 函数或使用 grep 进行了一些尝试，但这些方法即使没有失败也非常非常慢。

我也尝试过类似 Biostring 的将模式字典与引用匹配 函数的方法，但我已经在创建字典的步骤中失败了——这似乎是因为行中的序列非常不同。

(来自 Biostring 的错误信息。)

Error in .Call2("ACtree2_build", tb, pp_exclude, base_codes, nodebuf_ptr,  : 
  element 2 in Trusted Band has a different length than first element

有没有人知道如何实现我想要实现的目标？ 同样，一个挑战是数据框的大小超过 600 万行，并且基本上是针对数据框中的每一行测试每一行。

非常感谢您的任何反馈! 非常感谢!

信息添加 如果以下假设成立，是否有一种可行的方法:只有字符串匹配开头才有意义，并且至少有一个字符串必须匹配完整的字符序列。换句话说:在一个或多个不同行的字符串的开头可以找到一行的完整序列。

Answer 1

这是我为更简单的问题(查找作为其他序列的初始子串的序列)整理的内容。更一般的问题可以类似地解决，它只是更困惑并且需要(更多)更长的时间。对于下一步，我计划模拟您的数据(创建一个包含 600 万个长度在您提到的范围内的字符串的文本文件)并测试解决方案以查看需要多长时间。然后，正如我所说，我将在 Oracle 数据库中尝试同样的事情，看看是否存在巨大差异。只有“简单问题”在合理的时间内运行，“一般问题”才是合理的项目。

我假设数据有某种序列标识符(id 在数据库中是自然的)。您将在我在下面显示的输入文件中看到它。您还将看到输出的格式——对于作为较长序列的初始子字符串的每个较短序列，我都会显示两个序列(及其 ID)。注意 - 较短的序列，如 ACTAGC，可以是多个较长字符串的初始子串，如 ACTAGCTA 和 ACTAGCAGCA。我的输出只显示一个较长的序列，而不是所有较长的序列。

原则上，该算法很简单。按字母顺序对所有字符串进行排序，然后仅将每个字符串与下一个字符串进行比较。如果它不是子字符串，则它不能是数据集中任何其他字符串的子字符串。其余部分在 bash 中实现。

对于长度最多为 k 的 n 个序列，按字母顺序对所有序列进行排序是 O(kn log n) 并且检查每个字符串与下一个字符串的关系是 O(kn) - 这就是为什么这有机会处理您的 600 万行。

输入文件:

$ cat input_file
10010 ACATAAGAGTGATGATAGATAGATGCAGATGACAGATG
10011 ATAGAGATGAGACAGATGACAGAAGATAGATAGAGCAGATAG
10013 ATAGAGTAGAGAGAGAGTACAGATAGAGGAGAGAGATAGAC
10015 ACAGATAGCAGATAGACAGA
10016 ACAGATGACAGAAGATAGATAGA
10018 TAAGAGTGATGATAGATAGATGCAGA
10023 ATCACCGTTACAGATCG
10024 GTGATGATAGATAGATGCAGATGACAGATG
10025 ATAGAGTAGAGAGAGAGT
10030 TAAGAGTGATGATAGATAG
10044 TAAGAGTGATGATAGATAGATGCAATGA

编辑 - 下面的 BASH 脚本相当脑残，向所有人道歉。在本回答的最后，我将展示正确的方法；使用 SED，而不是 Shell 循环和读取命令

BASH 脚本 文件名:dupes.sh

#!/bin/bash
sort -k 2 input_file | 
{
    read key1 seq1
    while read key2 seq2
    do
        if [[ $(expr substr $seq2 1 ${#seq1}) == $seq1 ]]
        then
            echo ""
            echo "$key1 $seq1"
            echo "$key2 $seq2"
        fi  
        key1=$key2
        seq1=$seq2
    done
}

(我使用 echo 进行输出；您可能希望改为重定向到一个文件。)

调用和输出

$ ./dupes.sh

10025 ATAGAGTAGAGAGAGAGT
10013 ATAGAGTAGAGAGAGAGTACAGATAGAGGAGAGAGATAGAC

10030 TAAGAGTGATGATAGATAG
10044 TAAGAGTGATGATAGATAGATGCAATGA

编辑 - 正如我上面所说，虽然这是正确的答案，但解决方案很糟糕。这是在 bash 中执行此操作的正确方法。对于相同数量的输入数据(而不是 80 分钟)，此解决方案只需不到一分钟 (!!)。

sort -k 2 dna_sequences | sed -nE '{N; /^[^ ]+ ([^ ]+)\n[^ ]+ \1/p; D}'

输出可以重定向到一个文件，或者可以进一步处理(例如，我没有在每个匹配对之后添加一个换行符；这可以在输出的进一步处理中完成，或者通过其他方式，如果必要的)。