Как извлечь целые куски текста из профиля. Скрытый файл Markov .hmm

Question

В настоящее время я использую скрытый марковский подход для обнаружения вирусов в данных метагеномики. Я использую профиль, созданный Институтом Пастера на основе vFAM Питера Скуес-Кокса и др., 2014.

После использования профиля с HMMer и предоставления переведенных контигов в каждом кадре считывания, HMM смогли идентифицировать ожидаемые вирусы в положительном контроле. Тем не менее, множество совпадений (с оценкой 10 ^ -10 или менее для условных и независимых) соответствуют бактериальным областям со 100% -ной идентичностью и ~ 98% охватом согласно BLAST.

Эти ложноположительные имеют что-то общее: согласно НММ они соответствуют эндогенным ретровирусам или белкам гигантских вирусов (пример: Zn-зависимая алкогольдегидрогеназа, 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] и // символы и совпадения с именами в моем списке (ИМЯ FAM006805 как в заголовке).

Я ценю любые предложения. Спасибо!

Франциско Итурральде-Мартинес

Answer 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