PSP - HHblits 算法搜索 BFD 与 UniRef30 的结果分析 (bfd_uniref_hits.a3m)

2023-11-10 05:41
文章标签 算法 分析 搜索 bfd psp hits hhblits uniref30 uniref a3m

本文主要是介绍PSP - HHblits 算法搜索 BFD 与 UniRef30 的结果分析 (bfd_uniref_hits.a3m)，希望对大家解决编程问题提供一定的参考价值，需要的开发者们随着小编来一起学习吧！

欢迎关注我的CSDN：https://spike.blog.csdn.net/
本文地址：https://spike.blog.csdn.net/article/details/132047940

HHblits
MMseqs2 与 HHblits 的算法比较：

  • 蛋白质序列搜索算法 MMseqs2 与 HHblits 的搜索结果差异
  • HHblits 算法搜索 BFD 与 UniRef30 的结果分析

Paper: HHblits: lightning-fast iterative protein sequence searching by HMM-HMM alignment

  • 2011.12 Nature Methods

HHblits 是一种基于 HMM-HMM 对齐的迭代蛋白质序列搜索工具，可以快速、灵敏、准确地找到数据库中与查询序列相似的蛋白质序列。HHblits 的主要特点有：

  • 使用隐马尔可夫模型（HMM）来表示查询序列和数据库序列，而不是使用简单的氨基酸序列。HMM 是一种统计模型，可以捕捉序列的进化变化和保守区域，从而提高序列相似性的检测能力。
  • 使用一种离散化的预过滤器，可以快速地筛选出与查询序列有潜在相似性的数据库序列，从而减少了 HMM-HMM 对齐的计算量。这种预过滤器可以将搜索时间提高 2500 倍，而不损失搜索灵敏度。
  • 使用迭代的方法，可以利用上一轮搜索得到的多序列比对（MSA）来构建更精确的查询 HMM，来进行下一轮搜索。这样可以不断地扩展搜索范围，找到更多的相关序列。

在 bfd_uniref_hits.a3m 中，以 UniRef100 开头的序列描述，即来源于 UniRef30 数据库，即：

>UniRef100_A0A2X3SLP5 Uncharacterized protein n=2 Tax=Streptococcus equi TaxID=1336 RepID=A0A2X3SLP5_STRSZ
-TNtsSQEVDQVAQALELMFDNNVSTSNFKKYVNNNFSDSEIAIAELELESRISNSrsefrvawnemggCIAGKIRDEFFAMISVGTIVKYAQKKAWKELALVVLKFVKANGLKTNAIIVAGQLALWAVQCG

a3m 文件格式是用于表示蛋白质序列比对的文本格式，是 FASTA 格式的一种扩展。在 a3m 文件格式中，蛋白质序列用单个字母表示，其中，匹配用大写字母表示，删除用横线(-) 表示，插入用小写字母表示，去掉小写字母，只保留大写字母与横线(-) ，则长度与目标序列相同。

其他序列就是来源于 BFD，以 SRR 开头的数量较多，即：

>SRR5690606_14355373 
-----QIEELAAQLEFLMEEALiiENGERTfdfEKIENEFgkeVKDEIKMLTVDAQVWqvqpgaitlaanqPWKDCMVGAIKDHFGvALVTAaleGGLWAYLEKKAYKEAAKLLVKF----AVGTNAVGIAGTLIYYGGKCT
...
>SoimicmetaTmtLPB_FD_contig_61_1212631_length_209_multi_1_in_0_out_0_1 # 2 # 76 # 1 # ID=3042713_1;partial=10;start_type=Edge;rbs_motif=None;rbs_spacer=None;gc_cont=0.693
---NEEIEQLAADLEFLMEEAAiydEKGKVVnfdfDLLEERFgYVLELEMLKEEIEAYnattegd--NDeiqlfswksCMISALKGHFGvALIEValtGGLWSYLEKKAYKEAAKLLLKI----GIEGNVIGLTAFLTWYSVDCI
...
>APFre7841882654_1041346.scaffolds.fasta_scaffold441692_2 # 291 # 482 # 1 # ID=441692_2;partial=01;start_type=ATG;rbs_motif=AGGAG;rbs_spacer=5-10bp;gc_cont=0.609
-EIPLEAqgisiFGANHCDGareserliahelAHQWFgnsvtakrwrhiwlhegFACYAEWLWSEHSGdrsaDEWAhhfHEKLASSPQDLLLADPGPRDMFddrvykrgalTLHVLRRTLGDENFFALLKDWTSRHRHGSAVTD------DFTGLAANYTDQSLRPLWDAWLYS--TEVPALDAESX-------------------------------------------------------------------------------------------------------------------
...

因此，使用 MMseqs2 分别搜索 BFD 和 UniRef30，再合并，与使用 HHblits 一起搜索两个库的效果相同，同时也可以计算不同搜索算法的结果差异，即：

  • 测试序列来源于 CASP15 提供的官方序列。
  • MMseqs2 的搜索参数，num_iterations 是 3，sensitivity 是 8，即 i3s8，其余默认。
  • HHblits 使用 AF2 默认的搜索工具，数据库也保持一致。

即：

Img

源码如下：

#!/usr/bin/env python
# -- coding: utf-8 --
"""
Copyright (c) 2022. All rights reserved.
Created by C. L. Wang on 2023/8/1
"""
import osfrom myutils.project_utils import read_file
from root_dir import DATA_DIRclass BfdUniref30Overlap(object):"""计算 MMseqs2 与 HHblits 搜索结果之间的重叠度"""def __init__(self):pass@staticmethoddef count_seqs(data_lines):return len(data_lines) // 2 - 1def split_hhblits(self, data_lines):"""拆分 Uniref30 与 BFD 数据"""n = len(data_lines)n = n // 2 * 2u_desc_list, u_seq_list = [], []b_desc_list, b_seq_list = [], []for i in range(2, n, 2):desc = data_lines[i][1:]seq = data_lines[i+1]seq = seq.replace("-", "")if desc.startswith("UniRef100"):desc = desc.split(" ")[0]u_desc_list.append(desc)u_seq_list.append(seq)else:items = desc.split("|")if len(items) > 1:desc = items[1]b_desc_list.append(desc)b_seq_list.append(seq)unique_u_seq = len(list(set(u_seq_list)))unique_u_desc = len(list(set(u_desc_list)))print(f"[Info] unique uniref30 desc: {unique_u_desc} / {len(u_desc_list)}")print(f"[Info] unique uniref30 seqs: {unique_u_seq} / {len(u_seq_list)}")unique_b_seq = len(list(set(b_seq_list)))unique_b_desc = len(list(set(b_desc_list)))print(f"[Info] unique bfd desc: {unique_b_desc} / {len(b_desc_list)}")print(f"[Info] unique bfd seqs: {unique_b_seq} / {len(b_seq_list)}")return u_desc_list, u_seq_list, b_desc_list, b_seq_listdef process(self, uniref30_path, bfd_path, hhblits_path):assert os.path.isfile(uniref30_path) and os.path.isfile(bfd_path) and os.path.isfile(hhblits_path)print(f"[Info] mmseqs uniref30: {uniref30_path}")print(f"[Info] mmseqs bfd: {bfd_path}")print(f"[Info] hhblits bfd and uniref30: {hhblits_path}")uniref30_lines = read_file(uniref30_path)bfd_lines = read_file(bfd_path)hhblits_lines = read_file(hhblits_path)print(f"[Info] uniref30: {self.count_seqs(uniref30_lines)}, bfd: {self.count_seqs(bfd_lines)}, "f"hhblits: {self.count_seqs(hhblits_lines)}")# ---------- 统计 MMseqs2 Uniref30 ---------- #u_desc_list, u_seq_list = [], []n = len(uniref30_lines)n = n // 2 * 2for i in range(2, n, 2):u_name = uniref30_lines[i][1:].split("\t")[0]u_seq = uniref30_lines[i+1].replace("-", "")assert u_name.startswith("UniRef100")u_desc_list.append(u_name)u_seq_list.append(u_seq)assert len(u_desc_list) == self.count_seqs(uniref30_lines)unique_u_seq = len(list(set(u_seq_list)))unique_u_desc = len(list(set(u_desc_list)))print(f"[Info] unique uniref30 desc: {unique_u_desc} / {len(u_desc_list)}")print(f"[Info] unique uniref30 seqs: {unique_u_seq} / {len(u_seq_list)}")# ---------- 统计 MMseqs2 Uniref30 ---------- ## ---------- 统计 BFD Uniref30 ---------- #b_desc_list, b_seq_list = [], []n = len(bfd_lines)n = n // 2 * 2for i in range(2, n, 2):b_name = bfd_lines[i][1:].split("\t")[0]b_seq = bfd_lines[i+1].replace("-", "")b_desc_list.append(b_name)b_seq_list.append(b_seq)assert len(b_desc_list) == self.count_seqs(bfd_lines)unique_b_seq = len(list(set(b_seq_list)))unique_b_desc = len(list(set(b_desc_list)))print(f"[Info] unique bfd desc: {unique_b_desc} / {len(b_desc_list)}")print(f"[Info] unique bfd seqs: {unique_b_seq} / {len(b_seq_list)}")# ---------- 统计 BFD Uniref30 ---------- ## ---------- 统计 HHblits BFD and Uniref30 ---------- #hu_desc_list, hu_seq_list, hb_desc_list, hb_seq_list = self.split_hhblits(hhblits_lines)assert len(hu_desc_list) + len(hb_desc_list) == self.count_seqs(hhblits_lines)# ---------- 统计 HHblits BFD and Uniref30 ---------- ## ---------- 统计 交集 ---------- #num_u_desc_is = len(set(u_desc_list).intersection(set(hu_desc_list)))print(f"[Info] uniref30 desc intersection: {num_u_desc_is}")num_u_seqs_is = len(set(u_seq_list).intersection(set(hu_seq_list)))print(f"[Info] uniref30 seqs intersection: {num_u_seqs_is}")num_b_desc_is = len(set(b_desc_list).intersection(set(hb_desc_list)))print(f"[Info] bfd desc intersection: {num_b_desc_is}")num_b_seqs_is = len(set(b_seq_list).intersection(set(hb_seq_list)))print(f"[Info] bfd seqs intersection: {num_b_seqs_is}")# ---------- 统计 交集 ---------- ## ---------- 统计 总数 ---------- #unique_mmseqs_seqs = list(set(u_seq_list + b_seq_list))unique_hhblits_seqs = list(set(hu_seq_list + hb_seq_list))print(f"[Info] unique_mmseqs_ses: {len(unique_mmseqs_seqs)}")print(f"[Info] unique_hhblis_ses: {len(unique_hhblits_seqs)}")# ---------- 统计 总数 ---------- #def main():fasta_name = "T1104-D1_A117"# fasta_name = "T1137s8-D1_A251"# fasta_name = "T1188-D1_A573"# fasta_name = "T1157s1_A1029"uniref30_path = os.path.join(DATA_DIR, "overlap", fasta_name, f"{fasta_name}-uniref30-i3s8.a3m")bfd_path = os.path.join(DATA_DIR, "overlap", fasta_name, f"{fasta_name}-bfd-i3s8.a3m")hhblits_path = os.path.join(DATA_DIR, "overlap", fasta_name, "bfd_uniref_hits.a3m")buo = BfdUniref30Overlap()buo.process(uniref30_path, bfd_path, hhblits_path)if __name__ == "__main__":main()

这篇关于PSP - HHblits 算法搜索 BFD 与 UniRef30 的结果分析 (bfd_uniref_hits.a3m)的文章就介绍到这儿，希望我们推荐的文章对编程师们有所帮助！


http://www.chinasem.cn/article/380690。

相关文章

不懂推荐算法也能设计推荐系统

不懂推荐算法也能设计推荐系统

本文以商业化应用推荐为例，告诉我们不懂推荐算法的产品，也能从产品侧出发， 设计出一款不错的推荐系统。 相信很多新手产品，看到算法二字，多是懵圈的。 什么排序算法、最短路径等都是相对传统的算法（注：传统是指科班出身的产品都会接触过）。但对于推荐算法，多数产品对着网上搜到的资源，都会无从下手。特别当某些推荐算法 和 “AI”扯上关系后，更是加大了理解的难度。 但，不了解推荐算法，就无法做推荐系
阅读更多...
康拓展开（hash算法中会用到）

康拓展开（hash算法中会用到）

康拓展开是一个全排列到一个自然数的双射（也就是某个全排列与某个自然数一一对应） 公式： X=a[n]*(n-1)!+a[n-1]*(n-2)!+...+a[i]*(i-1)!+...+a[1]*0! 其中,a[i]为整数,并且0<=a[i]<i,1<=i<=n。（a[i]在不同应用中的含义不同）； 典型应用： 计算当前排列在所有由小到大全排列中的顺序，也就是说求当前排列是第
阅读更多...
认识、理解、分类——acm之搜索

认识、理解、分类——acm之搜索

普通搜索方法有两种：1、广度优先搜索；2、深度优先搜索； 更多搜索方法： 3、双向广度优先搜索； 4、启发式搜索（包括A*算法等）； 搜索通常会用到的知识点：状态压缩（位压缩，利用hash思想压缩）。
阅读更多...
hdu1240、hdu1253（三维搜索题）

hdu1240、hdu1253（三维搜索题）

1、从后往前输入，（x，y，z）; 2、从下往上输入，（y , z, x）; 3、从左往右输入，（z，x，y）； hdu1240代码如下： #include<iostream>#include<algorithm>#include<string>#include<stack>#include<queue>#include<map>#include<stdio.h>#inc
阅读更多...
性能分析之MySQL索引实战案例

性能分析之MySQL索引实战案例

文章目录 一、前言二、准备三、MySQL索引优化四、MySQL 索引知识回顾五、总结 一、前言 在上一讲性能工具之 JProfiler 简单登录案例分析实战中已经发现SQL没有建立索引问题，本文将一起从代码层去分析为什么没有建立索引？ 开源ERP项目地址：https://gitee.com/jishenghua/JSH_ERP 二、准备 打开IDEA找到登录请求资源路径位置
阅读更多...
csu 1446 Problem J Modified LCS （扩展欧几里得算法的简单应用）

csu 1446 Problem J Modified LCS （扩展欧几里得算法的简单应用）

这是一道扩展欧几里得算法的简单应用题，这题是在湖南多校训练赛中队友ac的一道题，在比赛之后请教了队友，然后自己把它a掉 这也是自己独自做扩展欧几里得算法的题目 题意：把题意转变下就变成了：求d1*x - d2*y = f2 - f1的解,很明显用exgcd来解 下面介绍一下exgcd的一些知识点：求ax + by = c的解 一、首先求ax + by = gcd(a,b)的解 这个
阅读更多...
综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

综合安防管理平台LntonAIServer视频监控汇聚抖动检测算法优势

LntonAIServer视频质量诊断功能中的抖动检测是一个专门针对视频稳定性进行分析的功能。抖动通常是指视频帧之间的不必要运动，这种运动可能是由于摄像机的移动、传输中的错误或编解码问题导致的。抖动检测对于确保视频内容的平滑性和观看体验至关重要。 优势 1. 提高图像质量 - 清晰度提升：减少抖动，提高图像的清晰度和细节表现力，使得监控画面更加真实可信。 - 细节增强：在低光条件下，抖
阅读更多...
【数据结构】——原来排序算法搞懂这些就行，轻松拿捏

【数据结构】——原来排序算法搞懂这些就行，轻松拿捏

前言：快速排序的实现最重要的是找基准值，下面让我们来了解如何实现找基准值 基准值的注释：在快排的过程中，每一次我们要取一个元素作为枢纽值，以这个数字来将序列划分为两部分。 在此我们采用三数取中法，也就是取左端、中间、右端三个数，然后进行排序，将中间数作为枢纽值。 快速排序实现主框架： //快速排序 void QuickSort(int* arr, int left, int rig
阅读更多...
poj 3974 and hdu 3068 最长回文串的O（n）解法(Manacher算法)

poj 3974 and hdu 3068 最长回文串的O（n）解法(Manacher算法)

求一段字符串中的最长回文串。 因为数据量比较大，用原来的O（n^2）会爆。 小白上的O(n^2)解法代码：TLE啦~ #include<stdio.h>#include<string.h>const int Maxn = 1000000;char s[Maxn];int main(){char e[] = {"END"};while(scanf("%s", s) != EO
阅读更多...
秋招最新大模型算法面试，熬夜都要肝完它

秋招最新大模型算法面试，熬夜都要肝完它

💥大家在面试大模型LLM这个板块的时候，不知道面试完会不会复盘、总结，做笔记的习惯，这份大模型算法岗面试八股笔记也帮助不少人拿到过offer ✨对于面试大模型算法工程师会有一定的帮助，都附有完整答案，熬夜也要看完，祝大家一臂之力 这份《大模型算法工程师面试题》已经上传CSDN，还有完整版的大模型 AI 学习资料，朋友们如果需要可以微信扫描下方CSDN官方认证二维码免费领取【保证100%免费
阅读更多...

Copyright China编程 23002807@qq.com

沪ICP备2023033072号-2