623_宏基因组基因聚类¶

# === CD-HIT安装 ===
conda install -c bioconda cd-hit  # conda安装
cd-hit -h  # 验证安装

# === MMseqs2安装 ===
conda install -c bioconda mmseqs2  # conda安装
mmseqs version  # 查看版本（2026年最新为v16+）

# === 源码编译MMseqs2（可选，追求最新版）===
git clone https://github.com/soedinglab/MMseqs2.git  # 克隆仓库
cd MMseqs2 && mkdir build && cd build  # 创建编译目录
cmake -DCMAKE_BUILD_TYPE=RELEASE -DCMAKE_INSTALL_PREFIX=. ..  # 配置
make -j8 && make install  # 编译安装

# === 宏基因组预测基因聚类 ===
cd-hit-est \
    -i predicted_genes.fna \      # 输入：所有样本预测的基因核苷酸序列
    -o nr_genes_95.fna \          # 输出：非冗余基因集
    -c 0.95 \                     # 相似度阈值95%（同种水平）
    -aS 0.9 \                     # 较短序列被覆盖≥90%
    -n 10 \                       # word size（c≥0.95时用10）
    -G 0 \                        # 局部比对模式（宏基因组推荐）
    -M 0 \                        # 不限制内存
    -T 16 \                       # 16线程
    -d 0                          # 不截断序列名

# === 输出文件说明 ===
# nr_genes_95.fna       — 代表序列（非冗余基因集）
# nr_genes_95.fna.clstr — 聚类信息文件

# === 查看聚类统计 ===
grep -c "^>" nr_genes_95.fna  # 非冗余基因数
grep -c "^>Cluster" nr_genes_95.fna.clstr  # 聚类数

# === 蛋白序列聚类 ===
cd-hit \
    -i predicted_proteins.faa \   # 输入：蛋白序列
    -o nr_proteins_90.faa \       # 输出：非冗余蛋白集
    -c 0.90 \                     # 相似度阈值90%
    -n 5 \                        # word size（c≥0.7时用5）
    -G 0 \                        # 局部比对
    -aS 0.8 \                     # 较短序列覆盖≥80%
    -M 0 \                        # 不限制内存
    -T 16 \                       # 16线程
    -d 0                          # 不截断序列名

# word size (-n) 必须与相似度阈值 (-c) 匹配：
# c ≥ 0.95 → n = 10-11（核苷酸），n = 5（蛋白）
# c ≥ 0.90 → n = 8-9（核苷酸），n = 5（蛋白）
# c ≥ 0.85 → n = 7（核苷酸），n = 5（蛋白）
# c ≥ 0.80 → n = 5-6（核苷酸），n = 5（蛋白）
# c ≥ 0.70 → n = 4-5（核苷酸），n = 4（蛋白）
# c ≥ 0.60 → 核苷酸不支持，n = 4（蛋白）
# c ≥ 0.50 → 核苷酸不支持，n = 3（蛋白）
# c ≥ 0.40 → 核苷酸不支持，n = 2（蛋白）

import re  # 正则表达式
from collections import Counter  # 计数器

# 解析CD-HIT聚类文件
def parse_cdhit_clstr(clstr_file):
    """解析.clstr文件，返回每个聚类的成员列表"""
    clusters = {}  # 存储聚类信息
    current_cluster = None  # 当前聚类ID

    with open(clstr_file) as f:  # 打开聚类文件
        for line in f:  # 逐行读取
            if line.startswith(">Cluster"):  # 新聚类开始
                current_cluster = line.strip()  # 记录聚类名
                clusters[current_cluster] = []  # 初始化成员列表
            else:
                # 提取序列名
                match = re.search(r'>(\S+)\.\.\.', line)  # 正则匹配序列名
                if match:
                    seq_name = match.group(1)  # 获取序列名
                    is_rep = line.strip().endswith("*")  # 星号标记代表序列
                    clusters[current_cluster].append(
                        (seq_name, is_rep)
                    )

    return clusters  # 返回聚类字典

# 使用
clusters = parse_cdhit_clstr("nr_genes_95.fna.clstr")
sizes = [len(v) for v in clusters.values()]  # 每个聚类的大小

print(f"总聚类数: {len(clusters)}")  # 打印聚类数
print(f"最大聚类: {max(sizes)} 条序列")  # 最大聚类
print(f"单序列聚类: {sizes.count(1)} 个")  # 单例数
print(f"平均聚类大小: {sum(sizes)/len(sizes):.1f}")  # 平均大小

# === MMseqs2 easy-cluster：一条命令完成聚类 ===
mmseqs easy-cluster \
    predicted_genes.fna \         # 输入序列文件
    mmseqs_cluster \              # 输出前缀
    tmp/ \                        # 临时目录
    --min-seq-id 0.95 \           # 相似度阈值95%
    -c 0.9 \                      # 覆盖度阈值90%
    --cov-mode 1 \                # 覆盖模式：1=较短序列覆盖
    --threads 16                  # 16线程

# === 输出文件 ===
# mmseqs_cluster_rep_seq.fasta  — 代表序列
# mmseqs_cluster_cluster.tsv   — 聚类成员关系（制表符分隔）
# mmseqs_cluster_all_seqs.fasta — 所有序列（按聚类排序）

# === Linclust：线性时间聚类，适合超大数据集 ===
mmseqs easy-linclust \
    predicted_genes.fna \         # 输入序列文件
    linclust_result \             # 输出前缀
    tmp/ \                        # 临时目录
    --min-seq-id 0.95 \           # 相似度阈值95%
    -c 0.9 \                      # 覆盖度阈值90%
    --cov-mode 1 \                # 覆盖模式
    --threads 16                  # 16线程

# Linclust vs easy-cluster：
# Linclust是O(n)时间复杂度——序列数翻倍，时间也只翻倍
# easy-cluster是O(n*k)——更准确但更慢
# 数据量>1亿条序列时Linclust是唯一选择

# === 分步运行MMseqs2聚类 ===

# 第1步：创建数据库
mmseqs createdb predicted_genes.fna queryDB  # 创建MMseqs2数据库

# 第2步：聚类
mmseqs cluster \
    queryDB \                     # 输入数据库
    clusterDB \                   # 聚类结果数据库
    tmp/ \                        # 临时目录
    --min-seq-id 0.95 \           # 相似度阈值
    -c 0.9 \                      # 覆盖度
    --cov-mode 1 \                # 覆盖模式
    --threads 16                  # 线程数

# 第3步：提取代表序列
mmseqs createsubdb clusterDB queryDB repDB  # 创建代表序列子数据库
mmseqs convert2fasta repDB rep_seqs.fasta   # 转为FASTA格式

# 第4步：输出聚类成员关系
mmseqs createtsv queryDB queryDB clusterDB cluster.tsv  # 输出TSV

# --cov-mode 参数说明：
# 0 = 双向覆盖：query和target都要满足覆盖度阈值
# 1 = target覆盖：较短序列被覆盖≥阈值（最常用）
# 2 = query覆盖：query被覆盖≥阈值
# 3 = target序列长度被覆盖≥阈值
# 5 = 较短序列被覆盖≥阈值（类似CD-HIT的-aS）

# 宏基因组推荐：--cov-mode 1 -c 0.9
# 含义：较短的基因序列至少90%被较长的代表序列覆盖

#!/bin/bash
# === 完整的宏基因组非冗余基因集构建流程 ===

# 变量设置
GENE_DIR="prodigal_output"          # Prodigal预测结果目录
THREADS=32                          # 线程数
IDENTITY=0.95                       # 相似度阈值
COVERAGE=0.9                        # 覆盖度阈值

# === 第1步：合并所有样本的预测基因 ===
echo "合并所有样本的预测基因..."
cat ${GENE_DIR}/*_genes.fna > all_predicted_genes.fna  # 合并核苷酸
cat ${GENE_DIR}/*_proteins.faa > all_predicted_proteins.faa  # 合并蛋白

echo "总基因数: $(grep -c '^>' all_predicted_genes.fna)"

# === 第2步：去除短序列（<100bp）===
echo "过滤短序列..."
seqkit seq -m 100 all_predicted_genes.fna > filtered_genes.fna  # 过滤<100bp
echo "过滤后基因数: $(grep -c '^>' filtered_genes.fna)"

# === 第3步：用MMseqs2 Linclust聚类 ===
echo "开始MMseqs2 Linclust聚类..."
mmseqs easy-linclust \
    filtered_genes.fna \              # 输入
    gene_catalog \                    # 输出前缀
    tmp_mmseqs/ \                     # 临时目录
    --min-seq-id ${IDENTITY} \        # 相似度95%
    -c ${COVERAGE} \                  # 覆盖度90%
    --cov-mode 1 \                    # 较短序列覆盖模式
    --threads ${THREADS}              # 线程数

# === 第4步：统计结果 ===
NR_GENES=$(grep -c '^>' gene_catalog_rep_seq.fasta)
TOTAL_GENES=$(grep -c '^>' filtered_genes.fna)
REDUCTION=$(echo "scale=1; (1 - ${NR_GENES}/${TOTAL_GENES}) * 100" | bc)

echo "=== 聚类结果 ==="
echo "输入基因数: ${TOTAL_GENES}"
echo "非冗余基因数: ${NR_GENES}"
echo "冗余率: ${REDUCTION}%"

# === 第5步：清理临时文件 ===
rm -rf tmp_mmseqs/  # 删除临时目录
echo "完成！非冗余基因集: gene_catalog_rep_seq.fasta"

# === 用非冗余基因集做丰度定量 ===

# 第1步：建立索引
bwa index gene_catalog_rep_seq.fasta  # BWA建索引

# 第2步：将每个样本的reads比对到基因集
for sample in samples/*_R1.fastq.gz; do
    name=$(basename "$sample" _R1.fastq.gz)  # 提取样本名

    # BWA比对
    bwa mem \
        -t 16 \                    # 16线程
        gene_catalog_rep_seq.fasta \  # 参考基因集
        ${sample} \                # 正向reads
        ${sample/_R1/_R2} \        # 反向reads
    | samtools sort -@ 4 -o ${name}.bam  # 排序

    # 统计比对数
    samtools idxstats ${name}.bam > ${name}_counts.tsv  # 基因计数
done

# 第3步：合并所有样本的计数
python3 -c "
import pandas as pd  # 导入pandas
import glob  # 文件匹配

# 读取所有样本的计数文件
count_files = sorted(glob.glob('*_counts.tsv'))  # 找所有计数文件
counts = {}  # 存储计数

for f in count_files:  # 遍历每个文件
    name = f.replace('_counts.tsv', '')  # 提取样本名
    df = pd.read_csv(f, sep='\t', header=None,
                     names=['gene', 'length', 'mapped', 'unmapped'])
    counts[name] = df.set_index('gene')['mapped']  # 取比对数

# 合并为矩阵
abundance = pd.DataFrame(counts)  # 创建丰度矩阵
abundance = abundance.fillna(0)  # 缺失值填0

# TPM归一化
gene_lengths = pd.read_csv(count_files[0], sep='\t', header=None,
                           names=['gene', 'length', 'mapped', 'unmapped'])
lengths = gene_lengths.set_index('gene')['length']  # 基因长度

rpk = abundance.div(lengths, axis=0) * 1e3  # reads per kilobase
tpm = rpk.div(rpk.sum(axis=0), axis=1) * 1e6  # TPM归一化

tpm.to_csv('gene_abundance_tpm.tsv', sep='\t')  # 保存TPM矩阵
print(f'基因丰度矩阵: {tpm.shape[0]} 基因 x {tpm.shape[1]} 样本')
"