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

Question

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

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

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

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

目标是：标记这些序列相同。（最后，我的目标是为每一行分配一个 序列 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 数据库中尝试同样的事情，看看是否存在巨大差异。只有 "easy problem" 在合理的时间内运行，"general problem" 才是一个合理的项目。

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

原则上，该算法是微不足道的。按字母顺序对所有字符串进行排序，然后仅将每个字符串与下一个字符串进行比较。如果它不是子字符串，则它不能是数据集中任何其他字符串的子字符串。其余的正在 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[^ ]+ /p; D}'

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

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

Solution needed for identifying partially matching strings (DNA sequences) in a data.frame with many rows

r

sed

dna-sequence