如何从 Profile Hidden Markov .hmm 文件中提取整个文本块答案

【问题标题】：How to extract entire chunks of text from a Profile Hidden Markov .hmm file如何从 Profile Hidden Markov .hmm 文件中提取整个文本块
【发布时间】：2019-06-02 08:19:35
【问题描述】：

我目前正在使用基于隐马尔可夫模型的方法来检测宏基因组数据中的病毒。我使用了由 Peter Skewes-Cox 等人在 2014 年基于 vFAM 由巴斯德研究所制作的配置文件。

在使用 HMMer 配置文件并在每个阅读框中提供翻译的重叠群后，HMM 能够识别阳性对照中的预期病毒。尽管如此，根据 BLAST，许多匹配（条件和独立的 evalue 为 10^-10 或更小）与具有 100% 同一性和约 98% 覆盖率的细菌区域匹配。

这些误报有一些共同点：根据 HMM，它们与内源性逆转录病毒或巨型病毒蛋白相匹配（例如：锌依赖性乙醇脱氢酶、ABC 转运蛋白等）。

因此，我决定看看是否可以从配置文件中删除这些条目，以便减少误报，并列出所有具有与逆转录病毒或巨型病毒相关注释的家族。

我在这里复制了我的个人资料的一部分作为解释：

HMMER3/f [3.1b2 | February 2015]
NAME  FAM007957
LENG  1078
ALPH  amino
RF    no
MM    no
CONS  yes
CS    no
MAP   yes
DATE  Fri Oct 12 20:02:22 2018
NSEQ  7
EFFN  0.591309
CKSUM 134316360
STATS LOCAL MSV      -12.5867  0.69540
STATS LOCAL VITERBI  -13.9281  0.69540
STATS LOCAL FORWARD   -6.9899  0.69540
HMM          A        C        D        E        F        G        H        I        K        L        M        N        P        Q        R        S        T        V        W        Y   
            m->m     m->i     m->d     i->m     i->i     d->m     d->d
  COMPO   2.52786  4.09835  2.76055  2.58333  3.30703  2.91930  3.80486  2.88354  2.60376  2.56225  3.71312  2.89938  3.51565  3.18472  2.93829  2.53713  2.89512  2.66587  4.91819  3.50321
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.16684  3.93795  2.00858  0.61958  0.77255  0.00000        *
//
HMMER3/f [3.1b2 | February 2015]
NAME  FAM006805
LENG  283
ALPH  amino
RF    no
MM    no
CONS  yes
CS    no
MAP   yes
DATE  Fri Oct 12 20:20:45 2018
NSEQ  8
EFFN  0.714844
CKSUM 174391985
STATS LOCAL MSV      -11.1126  0.70178
STATS LOCAL VITERBI  -11.7648  0.70178
STATS LOCAL FORWARD   -5.4313  0.70178
HMM          A        C        D        E        F        G        H        I        K        L        M        N        P        Q        R        S        T        V        W        Y   
            m->m     m->i     m->d     i->m     i->i     d->m     d->d
  COMPO   2.58563  4.40070  2.84295  2.49411  3.55282  3.12077  3.71148  2.77600  2.56241  2.36701  3.54429  2.93369  3.66844  3.05176  2.79705  2.67258  2.87961  2.67320  4.73491  3.80457
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.02701  4.02100  4.74335  0.61958  0.77255  0.00000        *
      1   3.09160  4.61822  4.21161  3.81854  3.28069  3.94629  4.51938  2.47147  3.57779  1.85500  1.11955  4.07700  4.40970  3.95105  3.76521  3.45517  3.40087  2.49434  5.14000  3.91374      1 m - - -
          2.68618  4.42225  2.77519  2.73123  3.46354  2.40513  3.72494  3.29354  2.67741  2.69355  4.24690  2.90347  2.73739  3.18146  2.89801  2.37887  2.77519  2.98518  4.58477  3.61503
          0.02701  4.02100  4.74335  0.61958  0.77255  0.48576  0.95510
//

我的问题是，我怎样才能对 HMMER3/f [3.1b2 | 2015 年 2 月] 和 // 字符并与我列表中的名称匹配（标题中的名称 FAM006805）。

感谢任何建议。谢谢！

弗朗西斯科·伊图拉尔德-马丁内斯

【问题讨论】：

如果文件始终具有相同的长度，您可以使用 sed 提取您感兴趣的行。例如，sed -n '8,10p' 将提取第 8 到第 10 行。

标签： bioinformatics hidden-markov-models

【解决方案1】：

解析文件是一种选择：

from __future__ import print_function
import re
IDs=['FAM006805']

with open('tp.hmm', 'rt') as inp:
  flag=0
  chunk=''
  with open('tp_mod.hmm', 'wt') as newfile:
    for line in inp:
      if re.match(r'^//', line) and flag==0:
        chunk+=line
        print(chunk, file=newfile)
        chunk=''
      elif re.match(r'^//', line) and flag==1:
        flag=0
        chunk=''

      chunk+=line
      if re.match(r'^NAME\s+', line):
        print(line)
        m = re.match(r'^NAME\s+(\w+)', line)
        tp_id=m.group(1).strip()
        print(tp_id)
        if tp_id in IDs:
          flag=1

【讨论】：

您可以更新IDs 列表以排除它们被写入新文件。