查找出现两次的模式，并且每次允许 <=2 不匹配

Posted 2023-03-15

【中文标题】查找出现两次的模式，并且每次允许 <=2 不匹配

【英文标题】：Find pattern that is present twice and allow <=2 mismatches on each

【发布时间】：2021-04-28 21:52:42

【问题描述】：

我有一个 400,000 次读取的 fastq 文件（所以速度很重要）。在序列中集成了应该出现两次的条形码。给定一个条形码，我想找到条形码出现两次且

TATCTTGTGGAAAGGACGAAACACCGAACACAAAGCATAGATGCGTTTAAGAGCTATGCTGGAAAACAGCATAGCAAGTTTAAATAAGGCAAATCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGAC TATCTTGTGGAAAGGGACGAAACACCGGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGAC TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACCGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGAC TATCTTGTGGAAAGGACGAAACACCGAGTCCGAGCAGAAGAAGAAGTTTAAGAGCTATGCTGGAAACAGCATAGCAAGTTTAAATAAGGCTAGTCCGTTATCAACTTGAAAAAGTGGCACCGAGTCGGTGCTTTTTTTATTCGACGATAGGGGTGGCAGGGGAGGCCGAGGAGGAAGAAGGGGAGGTGGCAGATTCGACCGATAGGTGGCGTAACTAGATCTTGAGACAAA

请注意，第四个序列中的第一个条形码缺少一个字符。我尝试过使用 biopython 和正则表达式，但考虑到我有 5k 个条形码，它太慢了。我想知道在 python 或 grep、awk 或其他任何东西中是否有可用的快速解决方案。谢谢。

【问题讨论】：

几年前我曾使用 FlexBar (github.com/seqan/flexbar) 完成类似的任务 - 值得一看。您可能会在bioinformatics.stackexchange.com 上获得此类特定领域问题的更好/最新答案，例如bioinformatics.stackexchange.com/questions/5375/…

您是否尝试在 python 中编译正则表达式，然后在所有字符串上运行它？

agrep -2 ATTCGACCGATAGG file

这个问题已经在这个论坛上被问过几次了。谷歌它或搜索档案。

没有一个能解决我的问题。请先仔细阅读我的问题。我不是这里的新手。

【参考方案1】：

使用 GNU awk：

 awk ' for (i=1;i<=NF;i++)  fnd=0;subs=$i;while (match(subs,"ATTCGACCGATAGG"))  subs=substr(subs,RSTART+RLENGTH);if (RSTART>0)  fnd++;print fnd   if (fnd <=2)  print $i   ' file

解释：

 awk ' for (i=1;i<=NF;i++)                            # Loop on each space delimited field
         fnd=0;                                         # Initialise fnd variable/counter
         subs=$i;                                       # Initialise substring variable
         while (match(subs,"ATTCGACCGATAGG"))  
           subs=substr(subs,RSTART+RLENGTH);            # Check for multiple matches of "ATTCGACCGATAGG" in subs.
           if (RSTART>0)  
              fnd++;                                    # Increment fnd if string found in subs
            
          
         if (fnd <=2)  
            print $i                                    # If found twice or less than twice print the field
         
         
       ' file

【讨论】：

感谢您的回答。抱歉，我不太了解 awk，它在代码的哪一部分评估允许的

是的，fnd 代表匹配数。

那么，有没有办法允许

match函数匹配条码，fnd统计每个字符串的匹配数。