hhblits 是 HMM-HMM（Hidden Markov Model to Hidden Markov Model）比對方法的一部分，也是 HMMER 軟件套件中的工具之一。與 hhsearch 類似，hhblits 也用于進行高效的蛋白質序列比對，特別擅長于檢測遠緣同源性。

hh-suite 源碼?https://github.com/soedinglab/hh-suite

以下是 hhblits 的一些主要特點和用途：

1. HMM-HMM比對： hhblits 使用隱藏馬爾可夫模型（HMM）來表示蛋白質家族的模型，從而能夠更全面地考慮多序列信息。這使得它相對于傳統的序列-序列比對方法更能捕捉蛋白質家族的模式和結構。

2. 迭代比對： hhblits 支持迭代比對的策略，通過在每一輪中更新模型，逐漸提高比對的靈敏度。這種策略對于發現相對遠緣同源蛋白質非常有用。

3. 遠緣同源性檢測： 與其他 HMM-HMM比對方法類似，hhblits 被設計用于檢測遠緣同源蛋白質，即那些在序列上相對較遠但在結構和功能上相似的蛋白質。

4. 數據庫搜索： 用戶可以使用 hhblits 在大型蛋白質數據庫中搜索與給定蛋白質序列相似的蛋白質。

5. 靈敏度和特異性： hhblits 的設計目標是在保持高靈敏度的同時，降低假陽性的比對，以提高比對的特異性。

import glob
import os
import subprocess
from typing import Any, List, Mapping, Optional, Sequencefrom absl import logging
import contextlib
import shutil
import tempfile
import time_HHBLITS_DEFAULT_P = 20
_HHBLITS_DEFAULT_Z = 500@contextlib.contextmanager
def tmpdir_manager(base_dir: Optional[str] = None):"""Context manager that deletes a temporary directory on exit."""tmpdir = tempfile.mkdtemp(dir=base_dir)try:yield tmpdirfinally:shutil.rmtree(tmpdir, ignore_errors=True)@contextlib.contextmanager
def timing(msg: str):logging.info('Started %s', msg)tic = time.time()yieldtoc = time.time()logging.info('Finished %s in %.3f seconds', msg, toc - tic)class HHBlits:"""Python wrapper of the HHblits binary."""def __init__(self,*,binary_path: str,databases: Sequence[str],n_cpu: int = 4,n_iter: int = 3,e_value: float = 0.001,maxseq: int = 1_000_000,realign_max: int = 100_000,maxfilt: int = 100_000,min_prefilter_hits: int = 1000,all_seqs: bool = False,alt: Optional[int] = None,p: int = _HHBLITS_DEFAULT_P,z: int = _HHBLITS_DEFAULT_Z):"""Initializes the Python HHblits wrapper.Args:binary_path: The path to the HHblits executable.databases: A sequence of HHblits database paths. This should be thecommon prefix for the database files (i.e. up to but not including_hhm.ffindex etc.)n_cpu: The number of CPUs to give HHblits.n_iter: The number of HHblits iterations.e_value: The E-value, see HHblits docs for more details.maxseq: The maximum number of rows in an input alignment. Note that thisparameter is only supported in HHBlits version 3.1 and higher.realign_max: Max number of HMM-HMM hits to realign. HHblits default: 500.maxfilt: Max number of hits allowed to pass the 2nd prefilter.HHblits default: 20000.min_prefilter_hits: Min number of hits to pass prefilter.HHblits default: 100.all_seqs: Return all sequences in the MSA / Do not filter the result MSA.HHblits default: False.alt: Show up to this many alternative alignments.p: Minimum Prob for a hit to be included in the output hhr file.HHblits default: 20.z: Hard cap on number of hits reported in the hhr file.HHblits default: 500. NB: The relevant HHblits flag is -Z not -z.Raises:RuntimeError: If HHblits binary not found within the path."""self.binary_path = binary_pathself.databases = databasesself.n_cpu = n_cpuself.n_iter = n_iterself.e_value = e_valueself.maxseq = maxseqself.realign_max = realign_maxself.maxfilt = maxfiltself.min_prefilter_hits = min_prefilter_hitsself.all_seqs = all_seqsself.alt = altself.p = pself.z = zdef query(self, input_fasta_path: str) -> List[Mapping[str, Any]]:"""Queries the database using HHblits."""with tmpdir_manager() as query_tmp_dir:a3m_path = os.path.join(query_tmp_dir, 'output.a3m')db_cmd = []for db_path in self.databases:db_cmd.append('-d')db_cmd.append(db_path)cmd = [self.binary_path,'-i', input_fasta_path,'-cpu', str(self.n_cpu),'-oa3m', a3m_path,'-o', '/dev/null','-n', str(self.n_iter),'-e', str(self.e_value),'-maxseq', str(self.maxseq),'-realign_max', str(self.realign_max),'-maxfilt', str(self.maxfilt),'-min_prefilter_hits', str(self.min_prefilter_hits)]if self.all_seqs:cmd += ['-all']if self.alt:cmd += ['-alt', str(self.alt)]if self.p != _HHBLITS_DEFAULT_P:cmd += ['-p', str(self.p)]if self.z != _HHBLITS_DEFAULT_Z:cmd += ['-Z', str(self.z)]cmd += db_cmdprint("cmd:",' '.join(cmd))logging.info('Launching subprocess "%s"', ' '.join(cmd))process = subprocess.Popen(cmd, stdout=subprocess.PIPE, stderr=subprocess.PIPE)with timing('HHblits query'):stdout, stderr = process.communicate()retcode = process.wait()if retcode:# Logs have a 15k character limit, so log HHblits error line by line.logging.error('HHblits failed. HHblits stderr begin:')for error_line in stderr.decode('utf-8').splitlines():if error_line.strip():logging.error(error_line.strip())logging.error('HHblits stderr end')raise RuntimeError('HHblits failed\nstdout:\n%s\n\nstderr:\n%s\n' % (stdout.decode('utf-8'), stderr[:500_000].decode('utf-8')))with open(a3m_path) as f:a3m = f.read()raw_output = dict(a3m=a3m,output=stdout,stderr=stderr,n_iter=self.n_iter,e_value=self.e_value)return [raw_output]if __name__ == "__main__":hhblits_binary_path = "/home/zheng/software/hhsuite-3.3.0-SSE2-Linux/bin/hhblits"input_fasta_path = "/home/zheng/test/Q94K49.fasta"#hhsuite數據庫下載地址：https://wwwuser.gwdg.de/~compbiol/data/hhsuite/databases/hhsuite_dbs/  scop70_database_path = "/home/zheng/database/scop70_1.75_hhsuite3/scop70_1.75"pdb70_database_path = "/home/zheng/database/pdb70_from_mmcif_latest/pdb70"## 單一數據庫#hhblits_runner = HHBlits(binary_path = hhblits_binary_path,  #                         databases = [scop70_database_path])## 多個數據庫database_lst = [scop70_database_path, pdb70_database_path]hhblits_runner = HHBlits(binary_path = hhblits_binary_path,databases = database_lst) result = hhblits_runner.query(input_fasta_path)print(len(result))print(result[0]['a3m'])#print(result)