计算允许不匹配的唯一句子的数量

Question

我有一个文本文件，里面有2亿句话。我想统计文件中特定类型句子的出现次数，并允许两个字符不匹配(这可能是插入重复字符或两个缺失字符(。字符总是A、G、C或T。不匹配字符的位置可能是随机的。我提供了一个小样本来说明我要解释的内容：

我有以下句子：

GTCGAGGTTCTCATCGCCTGGCCGCGCGTGTCTAGATCGGAAGAGCTCGT
GTTTAGGTTCTCATCGCCTGGCCGCGCGTGTCTAGATCGGAAGAGCTCGT
GTCGAGGTTCTCATCGCCTGGCCGCGCGTGTCTAGATCGGAAGAGCTAAT
TAACGTTCAGTTACGGCGTTGAGGTTTTACCTAAGATCGGAAGAGCTCGT
TCCGTAGCGCTCTGCTTCCAGTCGTGGCGGGGAGATCGGAAGAGCTCGTA
TACAAGACTTCATGAATAACGTGACTACGGAGATCGGAAGAGCTCGTATG
TAATGCCACTCCTCTCCCGACTGTTAACACTACTGGTTATATTGACCATG
CGACCTGGGTCAGCTCTGGAGTTTCGTTGAGTTAGATCGGAAGAGCTCGT
ATTTTGATAGTTTGACGGTTAATGCTGGTAATGGTGGTTTTCTTCATTGC
ACCCATGCCTACAGTATTGTTATCGGTAGCAAGCACATCACCTTGAATGC
GCAAGTTGCCATACAAAACAGGGTCGCCAGCAATATCGGTATAAGTCAAA
GAGTTCTAGTGTACGAGAGAGAGACGACGATGGAGATCGGAAGCGCTCTT
TGTTACTACAGGCATAATACGTGTTCCCGGATGAAGATCGGAAGAGCTCG
GACGACCAAAATTAGGGTCAACGCTACCTGTAGGAAGTGTCCGCATAAAG

例如，如果这是我看的第一句话

GTCGAGGTTCTCATCGCCTGGCCGCGCGTGTCTAGATCGGAAGAGCTCGT

然后是文件中的第二句话

GT**TT**AGGTTCTCATCGCCTGGCCGCGCGTGTCTAGATCGGAAGAGCTCGT

将与第一句相似，因为只有2个字符的差异。

然后是文件中的第三句话

GTCGAGGTTCTCATCGCCTGGCCGCGCGTGTCTAGATCGGAAGAGCT**AA**T

在这里，除了两个字符之外的最后一个字符都被修改了，但其他所有字符都与第一句相似。

因此，不变的句子可以是任何东西，然后将其余的句子与两个不匹配的句子进行比较，然后进行计数。差异可能是重复的字符、丢失的字符或新字符。最后，当你阅读文件时，输出将是一个特定的句子出现3次，允许不匹配。

GTCGAGGTTCTCATCGCCTGGCCGCGCGTGTCTAGATCGGAAGAGCTCGT 3 times

有没有办法让我用锥子？我试过了：

cat myfile.fq | 
awk '((NR-2)%4==0){character=$1;total++;count[character]++}END{for(character 
in count){if(!max||count[character]>max) 
{max=count[character];maxcharacter=character};if(count[character]==1){unique++}};print  total,unique}'

---

Ed Morton编辑-通过C美化程序运行上述awk脚本(https://codebeautify.org/c-formatter-beautifier)使其可读产生：

((NR - 2) % 4 == 0) {
character = $1;
total++;
count[character]++
}
END {
for (character in count) {
if (!max || count[character] > max) {
max = count[character];
maxcharacter = character
};
if (count[character] == 1) {
unique++
}
};
print total, unique
}

Answer 1

您的需求尚不明确，但我认为您可能会尝试做以下事情：

$ cat tst.awk
BEGIN {
tgtStr  = "APPLEISHEALTHY"
tgtLgth = length(tgtStr)
}
{
curStr  = $0
curLgth = length(curStr)
isMatch = 0
if ( curStr == tgtStr ) {
# curStr is tgtStr
# "APPLEISHEALTHY" vs "APPLEISHEALTHY"
isMatch = 1
}
else if ( curLgth == (tgtLgth-2) ) {
# curStr may be tgtStr minus 2 chars, e.g.
# "APPLEISHEALTHY" vs "APPLEISHEALT"
isMatch = 1
maxLgth = tgtLgth
curPos = tgtPos = 0
for (pos=1; pos<=maxLgth; pos++) {
curChar = substr(curStr,++curPos,1)
tgtChar = substr(tgtStr,++tgtPos,1)
if (curChar != tgtChar) {
if (curPos == tgtPos) {
# first char mismatch but curStr is 2 chars shorter
# than tgtStr so thats expected so advance tgtPos
# 1 char and back up curPos 1 char and continue.
curPos--
tgtPos++
}
else {
# still mismatching after first 2-char skip so fail
isMatch = 0
}
}
}
}
else if ( curLgth == tgtLgth ) {
# curStr may be tgtStr minus 2 chars plus 2 other chars, e.g.
# "APPLEISHEALTHY" vs "APPLEISHEALTXX"
}
else if ( curLgth == (tgtLgth+2) ) {
# curStr may be tgtStr plus 2 chars, e.g.
# "APPLEISHEALTHY" vs "APPLEISHEALTHYXX"
}
print curStr, (isMatch ? "is" : "is not"), "a match for", tgtStr
}

例如：

$ cat file
APPLEISHEALTHY
APPLEISALTHY
APPLEISXLTHY
$ awk -f tst.awk file
APPLEISHEALTHY is a match for APPLEISHEALTHY
APPLEISALTHY is a match for APPLEISHEALTHY
APPLEISXLTHY is not a match for APPLEISHEALTHY

你必须仔细思考以上逻辑，看看它是否正确，并为剩下的2种情况编写逻辑，但希望这能向你展示如何处理这个问题。

Answer 2

你似乎拥有的是测序机的"下一代短读">
就像我喜欢用awk破坏生物信息学问题以获取乐趣和利润一样。使用indels(插入/删除(，这非常适合使用正确的工具进行作业，这是局部对齐，如果你想让每个人都能尝试并真正复制结果，通常使用ncbi-blast
(除非这是家庭作业(

如果你在这个领域工作，安装blast并计算出实现目标的参数是你应该做的事情。