当前位置：首页 > article >正文

mmseqs2进行pdb蛋白质序列聚类分析

article 2024/12/24 11:26:49

mmseqs2是一款用于搜索和聚类大规模蛋白质和核酸序列集的开源软件套件。

https://github.com/soedinglab/MMseqs2
本示例将运行mmseqs2对pdb数据库的蛋白质序列进行聚类分析。

1. 安装 MMseqs2

conda activate bioinfo

conda install -c bioconda mmseqs2

2. 下载pdb序列数据到指定文件夹

wget https://files.rcsb.org/pub/pdb/derived_data/pdb_seqres.txt.gz

gunzip pdb_seqres.txt.gz

3. 过滤只有蛋白质的序列

def extract_protein_sequences(input_file, output_file):
    """从FASTA文件中提取蛋白质序列并保存到新文件中"""
    with open(input_file, 'r') as infile, open(output_file, 'w') as outfile:
        write_seq = False  # 标记是否要写序列
        for line in infile:
            if line.startswith(">"):
                # 解析序列头部信息，检查是否是蛋白质
                if "mol:protein" in line:
                    write_seq = True
                    outfile.write(line)  # 写入头部信息
                else:
                    write_seq = False
            elif write_seq:
                outfile.write(line)  # 写入序列内容

# 输入FASTA文件
input_fasta = "/Users/zhengxueming/test/test_input_sequence.fasta"
# 输出只包含蛋白质序列的FASTA文件
output_fasta = "/Users/zhengxueming/test/test_protein_sequence.fasta"

# 提取蛋白质序列
extract_protein_sequences(input_fasta, output_fasta)

print(f"蛋白质序列已提取并保存到 {output_fasta}")

数据为演示数据，注意改为正确的文件和路径

4. 创建数据库

mmseqs createdb pdb_protein_sequence.fasta pdb_sequencesDB

5. 执行聚类分析

mmseqs cluster pdb_sequencesDB pdb_clustersDB_0.3 ../tmp --min-seq-id 0.3

这里 ../tmp是临时文件夹，你可以随意指定一个路径。--min-seq-id 0.3 表示聚类时最小序列相似性为 30%。

聚类参数

--min-seq-id：控制序列的相似度阈值，数值范围为 0 到 1。
--cov-mode：设置覆盖率模式，例如 --cov-mode 1 表示至少有一个序列完全覆盖。
--cluster-mode：设置聚类模式，不同模式会影响聚类的方式（如单链、链式聚类等）。

6. 提取聚类结果

mmseqs createtsv pdb_sequencesDB pdb_sequencesDB pdb_clustersDB_0.3 result.tsv

7. 查看cluster

#查看多少cluster
awk '{print $1}' result.tsv|sort|uniq|wc -l
# 两列，前面一列是该cluster的代表序列pdb号

8. 聚类结果的进一步分析

你可以将聚类结果与其他工具结合使用，进一步分析聚类的特性，或者使用 mmseqs result2flat 命令将聚类序列提取出来，进行下游分析：

mmseqs result2flat pdb_sequencesDB pdb_sequencesDB pdb_clustersDB_0.3 cluster_result.txt

结果文件解读：

>4mwa_C：表示聚类的代表序列（聚类中心，通常是一个代表该聚类的 PDB 序列），4mwa_C 是 PDB 序列的标识符，其中 4mwa 是 PDB ID，C 是链 ID。
数字列表：这些是与聚类中心相对应的其他序列的索引号。每一个数字对应于聚类数据库（pdb_clustersDB_0.3）中的某一个序列，这些序列与聚类中心的相似度在 0.3 阈值下聚在一起。