我目前正在使用一种基于隐马尔可夫模型的方法来检测元基因组数据中的病毒。我使用的是巴斯德研究所基于Peter Skewes-Cox等人,2014年的vFAMs制作的个人资料。
在与HMMer一起使用profile并在每个阅读框架中提供翻译的重叠群后,HMMs能够在阳性对照中识别预期的病毒。尽管如此,根据BLAST,许多匹配(条件和独立的evalue都在10^-10或更低)与细菌区域匹配,具有100%的同一性和~98%的覆盖率。
这些假阳性有一些共同之处:根据它们与内源性逆转录病毒或巨型病毒蛋白(例如:依赖锌的酒精脱氢酶、ABC转运蛋白等)匹配的HMM。
因此,我决定看看是否可以从配置文件中删除这些条目,以便减少误报,并列出了具有与逆转录病毒或巨型病毒相关的注释的所有家族。
我在这里复制了我的个人资料的一大块作为解释:
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 )相匹配。
我很感谢你的建议。谢谢!
弗朗西斯科·伊图拉尔德-马丁内斯
发布于 2019-01-18 06:09:33
解析文件是一种选择:
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=1https://stackoverflow.com/questions/54078268
复制相似问题