蛋白质家族分类Pfam¶

蛋白质家族分类的逻辑：

想象一个"姓氏数据库"：
  - "张"家有100万人 → 他们都有某些共同特征
  - "李"家有80万人 → 他们有另一些共同特征

Pfam做的事：
  - "激酶域"家族有10万个蛋白 → 它们都能磷酸化底物
  - "免疫球蛋白域"家族有5万个蛋白 → 它们都有类似的折叠

怎么识别家族成员？
  传统方法：BLAST比对（像用照片认人）
  Pfam方法：HMM模型（像用特征模型识别：身高+脸型+发色...）

HMM的优势：
  - 不是跟单个序列比，而是跟整个家族的"平均画像"比
  - 能识别远缘同源蛋白（BLAST找不到的）
  - 更低的假阳性率

Pfam数据（2024版本）：
  Pfam-A条目: ~20,795个蛋白质家族
  覆盖序列:    UniProtKB中80%+的蛋白质
  现已整合到InterPro中（维护由InterPro团队接管）

# HMMER是Pfam背后的搜索引擎
# 安装HMMER
conda install -c bioconda hmmer  # conda安装（推荐）

# 或者源码安装
wget http://eddylab.org/software/hmmer/hmmer-3.4.tar.gz  # 下载源码
tar -xzf hmmer-3.4.tar.gz  # 解压
cd hmmer-3.4  # 进入目录
./configure --prefix=/usr/local  # 配置安装路径
make  # 编译
make install  # 安装

# 验证安装
hmmsearch -h  # 查看帮助信息

# 下载Pfam-A HMM数据库
wget https://ftp.ebi.ac.uk/pub/databases/Pfam/current_release/Pfam-A.hmm.gz  # 下载HMM库
gunzip Pfam-A.hmm.gz  # 解压
hmmpress Pfam-A.hmm  # 建立索引（搜索前必须做！）

# ========== hmmscan：一条蛋白质查所有Pfam家族 ==========
# "这个蛋白质属于哪些家族？"

hmmscan \
    --domtblout results.domtab \  # 域表格输出（易解析）
    --tblout results.tbl \        # 序列表格输出
    -E 1e-5 \                     # 全序列E-value阈值
    --domE 1e-5 \                 # 域E-value阈值
    --cpu 4 \                     # 使用4个CPU
    Pfam-A.hmm \                  # Pfam HMM数据库
    query_protein.fasta           # 输入蛋白质序列

# 查看结果
head -5 results.domtab  # 查看域匹配结果

# 结果列说明（domtblout格式）：
# 列1:  target name     → Pfam家族ID（如PF00069）
# 列2:  target accession → Pfam条目号
# 列4:  query name      → 你的蛋白质名称
# 列7:  full E-value    → 全序列E-value
# 列12: domain E-value  → 域E-value（更重要）
# 列18: hmm from        → HMM模型匹配起始
# 列19: hmm to          → HMM模型匹配结束
# 列20: ali from        → 序列上的匹配起始
# 列21: ali to          → 序列上的匹配结束
# 列23: description     → 家族描述

# ========== hmmsearch：一个Pfam家族查哪些蛋白质包含它 ==========
# "哪些蛋白质属于激酶家族？"

# 先提取特定Pfam家族的HMM
hmmfetch Pfam-A.hmm PF00069 > kinase.hmm  # 提取蛋白激酶域的HMM

# 搜索蛋白质数据库
hmmsearch \
    --domtblout kinase_hits.domtab \  # 域表格输出
    --tblout kinase_hits.tbl \        # 序列表格输出
    -E 1e-10 \                        # E-value阈值
    --cpu 8 \                         # 8个CPU
    kinase.hmm \                      # 激酶域HMM模型
    all_proteins.fasta                # 蛋白质序列数据库

#!/usr/bin/env python3
"""解析HMMER/Pfam搜索结果"""

import pandas as pd  # 数据处理
import re  # 正则表达式

def parse_domtblout(filename):
    """解析hmmscan/hmmsearch的domtblout格式结果"""

    results = []
    with open(filename, 'r') as f:
        for line in f:
            if line.startswith("#"):  # 跳过注释行
                continue

            fields = line.split()  # 按空格分隔
            if len(fields) < 23:  # 字段数不够就跳过
                continue

            result = {
                "pfam_id": fields[1],  # Pfam条目号（如PF00069.29）
                "pfam_name": fields[0],  # Pfam家族名称
                "query": fields[3],  # 查询蛋白质名称
                "query_len": int(fields[5]),  # 蛋白质长度
                "full_evalue": float(fields[6]),  # 全序列E-value
                "full_score": float(fields[7]),  # 全序列得分
                "dom_evalue": float(fields[11]),  # 域E-value
                "dom_score": float(fields[13]),  # 域得分
                "hmm_from": int(fields[15]),  # HMM匹配起始
                "hmm_to": int(fields[16]),  # HMM匹配结束
                "seq_from": int(fields[17]),  # 序列匹配起始
                "seq_to": int(fields[18]),  # 序列匹配结束
                "description": " ".join(fields[22:])  # 家族描述
            }
            results.append(result)

    df = pd.DataFrame(results)
    return df

# ========== 使用示例 ==========
df = parse_domtblout("results.domtab")  # 解析结果
print(f"总命中数: {len(df)}")

# 统计最常见的Pfam家族
family_counts = df["pfam_name"].value_counts()
print("\n最常见的Pfam家族 TOP10：")
print(family_counts.head(10))

# 统计每个蛋白质的域数量
protein_domains = df.groupby("query")["pfam_name"].count()
print(f"\n每个蛋白质平均域数: {protein_domains.mean():.1f}")

# 导出域组成结构
domain_arch = df.groupby("query").apply(
    lambda x: " + ".join(x.sort_values("seq_from")["pfam_name"])
).reset_index()
domain_arch.columns = ["protein", "domain_architecture"]
print("\n蛋白质域组成结构：")
print(domain_arch.head(10).to_string(index=False))

# 保存结果
df.to_csv("pfam_results.csv", index=False)
print("\n结果已保存到 pfam_results.csv")

#!/usr/bin/env python3
"""通过InterPro API查询Pfam家族信息"""

import requests  # HTTP请求

# ========== 查询特定Pfam家族 ==========
def query_pfam_family(pfam_id):
    """查询Pfam家族详细信息"""

    # Pfam现已整合到InterPro，通过InterPro API查询
    url = f"https://www.ebi.ac.uk/interpro/api/entry/pfam/{pfam_id}"
    response = requests.get(url, headers={"Accept": "application/json"})

    if response.status_code == 200:
        data = response.json()
        metadata = data.get("metadata", {})

        print(f"\n===== Pfam家族: {pfam_id} =====")
        print(f"名称: {metadata.get('name', 'N/A')}")
        print(f"类型: {metadata.get('type', 'N/A')}")
        print(f"描述: {metadata.get('description', ['N/A'])[0] if metadata.get('description') else 'N/A'}")
        print(f"包含蛋白质数: {metadata.get('counters', {}).get('proteins', 'N/A')}")
        print(f"对应InterPro: {metadata.get('integrated', 'N/A')}")

        return data
    else:
        print(f"查询失败: HTTP {response.status_code}")
        return None

# ========== 查询蛋白质的Pfam域 ==========
def get_pfam_domains(uniprot_id):
    """获取蛋白质包含的所有Pfam域"""

    url = f"https://www.ebi.ac.uk/interpro/api/entry/pfam/protein/uniprot/{uniprot_id}"
    response = requests.get(url, headers={"Accept": "application/json"})

    if response.status_code == 200:
        data = response.json()
        results = data.get("results", [])

        print(f"\n蛋白质 {uniprot_id} 的Pfam域：")
        print(f"{'Pfam ID':<15}{'名称':<30}{'InterPro'}")
        print("-" * 60)

        for entry in results:
            meta = entry.get("metadata", {})
            print(f"{meta.get('accession', 'N/A'):<15}"
                  f"{meta.get('name', 'N/A'):<30}"
                  f"{meta.get('integrated', 'N/A')}")

        return results
    else:
        print(f"查询失败: HTTP {response.status_code}")
        return []

# ========== 使用示例 ==========
# 查询蛋白激酶域家族
query_pfam_family("PF00069")

# 查询TP53的所有Pfam域
get_pfam_domains("P04637")

# ========== 从多序列比对构建自己的HMM模型 ==========

# 第一步：准备多序列比对
# 假设你有一组同源蛋白质序列
mafft --auto homolog_sequences.fasta > aligned.afa  # 用MAFFT做多序列比对

# 第二步：构建HMM模型
hmmbuild \
    my_family.hmm \   # 输出HMM模型文件
    aligned.afa        # 输入多序列比对文件

# 第三步：用自建模型搜索
hmmsearch \
    --domtblout custom_hits.domtab \  # 输出文件
    -E 1e-10 \                        # E-value阈值
    my_family.hmm \                   # 自建HMM模型
    target_proteome.fasta             # 要搜索的蛋白质组

========================================
蛋白质家族分类Pfam 速查表
========================================

【HMMER核心命令】
hmmscan              → 用序列搜Pfam库（这蛋白属于哪些家族？）
hmmsearch            → 用HMM搜序列库（哪些蛋白属于这个家族？）
hmmbuild             → 从比对构建HMM模型
hmmpress             → 对HMM库建索引（搜索前必做）
hmmfetch             → 从HMM库提取特定模型

【常用参数】
-E                   → 全序列E-value阈值（推荐1e-5）
--domE               → 域E-value阈值（推荐1e-5）
--domtblout          → 域表格输出（最易解析）
--tblout             → 序列表格输出
--cpu                → CPU核数

【E-value判读】
< 1e-10              → 非常可靠的匹配
1e-10 ~ 1e-5         → 可靠的匹配
1e-5 ~ 0.01          → 需要进一步确认
> 0.01               → 可能是假阳性

【Pfam ID格式】
PF00069              → 蛋白激酶域
PF00071              → Ras GTPase家族
PF00076              → RNA识别域(RRM)
PF00400              → WD40重复

【数据库下载】
Pfam-A HMM           → ftp.ebi.ac.uk/pub/databases/Pfam/
InterPro整合         → www.ebi.ac.uk/interpro/

【面试考点】
Q: HMM和BLAST的区别？
A: BLAST是序列-序列比对，HMM是序列-模型比对，能识别更远的同源关系

Q: Pfam和InterPro的关系？
A: Pfam是InterPro整合的子库之一，Pfam维护已由InterPro团队接管

Q: 什么是域组成结构(domain architecture)？
A: 蛋白质上各功能域的排列顺序，如"SH2 + kinase + SH3"
========================================