首页
学习
活动
专区
圈层
工具
发布
社区首页 >问答首页 >如何从配置文件隐马尔可夫.hmm文件中提取完整的文本块

如何从配置文件隐马尔可夫.hmm文件中提取完整的文本块
EN

Stack Overflow用户
提问于 2019-01-08 00:32:50
回答 1查看 47关注 0票数 2

我目前正在使用一种基于隐马尔可夫模型的方法来检测元基因组数据中的病毒。我使用的是巴斯德研究所基于Peter Skewes-Cox等人,2014年的vFAMs制作的个人资料。

在与HMMer一起使用profile并在每个阅读框架中提供翻译的重叠群后,HMMs能够在阳性对照中识别预期的病毒。尽管如此,根据BLAST,许多匹配(条件和独立的evalue都在10^-10或更低)与细菌区域匹配,具有100%的同一性和~98%的覆盖率。

这些假阳性有一些共同之处:根据它们与内源性逆转录病毒或巨型病毒蛋白(例如:依赖锌的酒精脱氢酶、ABC转运蛋白等)匹配的HMM。

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

我在这里复制了我的个人资料的一大块作为解释:

代码语言:javascript
复制
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 )相匹配。

我很感谢你的建议。谢谢!

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

EN

回答 1

Stack Overflow用户

发布于 2019-01-18 06:09:33

解析文件是一种选择:

代码语言:javascript
复制
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
票数 0
EN
页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持
原文链接:

https://stackoverflow.com/questions/54078268

复制
相关文章

相似问题

领券
问题归档专栏文章快讯文章归档关键词归档开发者手册归档开发者手册 Section 归档