宏基因组组装后分析（Post-Assembly Analysis）¶

# 基本运行
prodigal -i assembled_contigs.fa \
         -o genes.gff \          # GFF格式基因坐标
         -a proteins.faa \       # 蛋白质序列
         -d genes.fna \          # 核苷酸CDS序列
         -p meta \               # 宏基因组模式（必须！）
         -f gff \                # 输出格式
         -q                      # 静默模式

# 统计预测结果
echo "Predicted genes: $(grep -c '>' proteins.faa)"
echo "Complete genes: $(grep -c 'partial=00' genes.gff)"
echo "Partial genes (5'): $(grep -c 'partial=10' genes.gff)"
echo "Partial genes (3'): $(grep -c 'partial=01' genes.gff)"
echo "Partial genes (both): $(grep -c 'partial=11' genes.gff)"

# 过滤短基因（可选，<100aa通常不可靠）
seqkit seq -m 100 proteins.faa > proteins_filtered.faa

# 多样本并行处理
for sample in sample1 sample2 sample3; do
  prodigal -i ${sample}_contigs.fa \
           -a ${sample}_proteins.faa \
           -d ${sample}_genes.fna \
           -p meta -f gff -q \
           -o ${sample}_genes.gff
done

# 合并所有样本的预测基因（添加样本前缀防止ID冲突）
for sample in sample1 sample2 sample3; do
  sed "s/^>/>${sample}_/" ${sample}_genes.fna >> all_genes.fna
  sed "s/^>/>${sample}_/" ${sample}_proteins.faa >> all_proteins.faa
done

# === 核苷酸水平聚类（推荐，保留定量能力）===
cd-hit-est -i all_genes.fna \
           -o nr_genes.fna \
           -c 0.95 \          # 95%序列一致性
           -aS 0.9 \          # 90%短序列覆盖度
           -n 10 \            # word size（c>0.9时用10）
           -d 0 \             # 描述信息全部保留
           -M 64000 \         # 内存限制(MB)
           -T 16 \            # 线程数
           -g 1               # 精确聚类模式

# === 蛋白质水平聚类 ===
cd-hit -i all_proteins.faa \
       -o nr_proteins.faa \
       -c 0.95 \
       -aS 0.9 \
       -n 5 \               # word size（c>0.9时用5）
       -d 0 \
       -M 64000 \
       -T 16 \
       -g 1

# === 大数据集分步聚类策略 ===
# 先以较低阈值粗聚，再精聚
cd-hit-est -i all_genes.fna -o step1.fna -c 0.99 -aS 0.9 -n 10 -M 64000 -T 16
cd-hit-est -i step1.fna -o nr_genes.fna -c 0.95 -aS 0.9 -n 10 -M 64000 -T 16

# === 使用MMseqs2（更快，适合超大数据集）===
mmseqs easy-linclust all_genes.fna nr_genes tmp \
  --min-seq-id 0.95 \
  -c 0.9 \
  --cov-mode 1 \
  --threads 32

# 统计结果
echo "Total input genes: $(grep -c '>' all_genes.fna)"
echo "Non-redundant genes: $(grep -c '>' nr_genes.fna)"
echo "Redundancy rate: $(echo "scale=2; (1-$(grep -c '>' nr_genes.fna)/$(grep -c '>' all_genes.fna))*100" | bc)%"

# 解析聚类文件
# 获取每个cluster的成员数
grep "^>" nr_genes.fna.clstr | sed 's/.*Cluster //' | sort -n | tail

# === 1. eggNOG-mapper（一站式注释）===
# 翻译非冗余核苷酸为蛋白质（如果之前是核苷酸聚类）
# 或直接使用蛋白质聚类结果
emapper.py -i nr_proteins.faa \
           -m diamond \
           --itype proteins \
           --data_dir /data/eggnog_db/ \
           --cpu 32 \
           --output nr_gene_emapper \
           --output_dir ./annotation/ \
           --override

# === 2. KofamScan（KEGG KO注释，比eggNOG的KEGG注释更准确）===
# 下载KOfam数据库
wget ftp://ftp.genome.jp/pub/db/kofam/ko_list.gz
wget ftp://ftp.genome.jp/pub/db/kofam/profiles.tar.gz
gunzip ko_list.gz; tar xzf profiles.tar.gz

# 运行
exec_annotation -o kofam_results.txt \
  --profile profiles/ \
  --ko-list ko_list \
  --cpu 32 \
  --tmp-dir /tmp/kofam/ \
  nr_proteins.faa

# 提取显著KO（带*号的）
grep "^\*" kofam_results.txt | awk '{print $2"\t"$3}' > significant_ko.tsv

# === 3. dbCAN（CAZyme注释）===
# 使用run_dbcan
run_dbcan nr_proteins.faa protein \
  --db_dir /data/dbCAN_db/ \
  --out_dir ./dbcan_results/ \
  --dia_cpu 16 \
  --hmm_cpu 16 \
  --tf_cpu 16

# 推荐取至少2种方法一致的结果
awk '$5 >= 2' dbcan_results/overview.txt > dbcan_consensus.txt

# === 4. CARD（抗性基因）===
rgi main -i nr_proteins.faa \
  -o rgi_results \
  -t protein \
  -a DIAMOND \
  --clean \
  -n 16

# === 5. VFDB（毒力因子）===
diamond blastp -d /data/VFDB/VFDB_setA_pro.dmnd \
  -q nr_proteins.faa \
  -o vfdb_results.txt \
  --outfmt 6 \
  -e 1e-5 \
  --id 40 \
  --query-cover 70 \
  -p 16

# ========== 方法1：Salmon（推荐）==========
# 1. 建立Salmon索引
salmon index -t nr_genes.fna -i salmon_index -p 16

# 2. 各样本定量
for sample in sample1 sample2 sample3; do
  salmon quant -i salmon_index \
    -l A \
    -1 ${sample}_clean_R1.fastq.gz \
    -2 ${sample}_clean_R2.fastq.gz \
    -o salmon_quant/${sample} \
    -p 16 \
    --meta \
    --validateMappings
done

# 3. 合并各样本的TPM/count
# 提取TPM
for sample in sample1 sample2 sample3; do
  tail -n +2 salmon_quant/${sample}/quant.sf | \
    awk -v s=${sample} 'BEGIN{OFS="\t"}{print $1, $4}' > salmon_quant/${sample}_tpm.txt
done

# 合并成矩阵
paste salmon_quant/*_tpm.txt | \
  awk 'BEGIN{OFS="\t"}{printf $1; for(i=2;i<=NF;i+=2) printf "\t"$i; print ""}' > gene_tpm_matrix.tsv

# ========== 方法2：Bowtie2 + SAMtools ==========
# 1. 建索引
bowtie2-build --threads 16 nr_genes.fna bt2_index

# 2. 比对
for sample in sample1 sample2 sample3; do
  bowtie2 -x bt2_index \
    -1 ${sample}_clean_R1.fastq.gz \
    -2 ${sample}_clean_R2.fastq.gz \
    -S ${sample}_mapped.sam \
    --no-unal \
    -p 16 \
    --very-sensitive

  samtools view -bS -@ 8 ${sample}_mapped.sam | \
    samtools sort -@ 8 -o ${sample}_mapped_sorted.bam
  samtools index ${sample}_mapped_sorted.bam
  rm ${sample}_mapped.sam
done

# 3. 统计每个基因的count
# 使用featureCounts（需要GTF/SAF格式的注释）
# 先将GFF转SAF
awk -F'\t' 'BEGIN{OFS="\t"; print "GeneID\tChr\tStart\tEnd\tStrand"} 
  /^[^#]/ && $3=="CDS" {
    match($9, /ID=([^;]+)/, arr); 
    print arr[1], $1, $4, $5, $7
  }' genes.gff > genes.saf

featureCounts -a genes.saf -F SAF \
  -o gene_counts.txt \
  -T 16 \
  -p --countReadPairs \
  sample1_mapped_sorted.bam \
  sample2_mapped_sorted.bam \
  sample3_mapped_sorted.bam

# ========== 方法3：CoverM ==========
# 计算基因覆盖度
coverm genome \
  --bam-files sample1.bam sample2.bam sample3.bam \
  --genome-fasta-files nr_genes.fna \
  --methods rpkm count covered_bases \
  --threads 16 \
  > coverm_results.tsv

# R语言构建KO丰度表
library(tidyverse)

# 读取基因丰度矩阵（TPM）
gene_tpm <- read_tsv("gene_tpm_matrix.tsv")

# 读取KO注释
ko_annotation <- read_tsv("significant_ko.tsv",
                          col_names = c("gene_id", "KO"))

# 合并
gene_ko_tpm <- gene_tpm %>%
  inner_join(ko_annotation, by = c("gene" = "gene_id"))

# 按KO汇总丰度（同一KO的基因TPM求和）
ko_abundance <- gene_ko_tpm %>%
  group_by(KO) %>%
  summarise(across(starts_with("sample"), sum, na.rm = TRUE)) %>%
  ungroup()

# 保存KO丰度表
write_tsv(ko_abundance, "ko_abundance_matrix.tsv")

# 同理构建COG丰度表
cog_annotation <- read_tsv("annotation/nr_gene_emapper.emapper.annotations",
                            comment = "##") %>%
  select(query, COG_category) %>%
  filter(COG_category != "-")

# KEGG Module完整性评估
# 读取KEGG Module定义（哪些KO组成一个Module）
# 计算每个样本中每个Module的完整度

# Python版本
import pandas as pd
import numpy as np

# 读取
gene_tpm = pd.read_csv("gene_tpm_matrix.tsv", sep="\t", index_col=0)
ko_anno = pd.read_csv("significant_ko.tsv", sep="\t", names=["gene_id", "KO"])

# 合并并汇总
merged = gene_tpm.merge(ko_anno, left_index=True, right_on="gene_id")
ko_abundance = merged.groupby("KO").sum(numeric_only=True)
ko_abundance.to_csv("ko_abundance_matrix.tsv", sep="\t")

# KEGG Pathway丰度
# 从KO映射到Pathway（使用KEGG REST API或本地文件）

# 使用DESeq2对功能丰度做差异分析
library(DESeq2)

# 将TPM转为近似count（如果用Salmon count输出则直接用）
# 或使用ALDEx2/MaAsLin2（更适合组成性数据）
library(ALDEx2)

# KO丰度矩阵（行=KO，列=样本）
ko_matrix <- as.matrix(read.csv("ko_count_matrix.csv", row.names = 1))
conditions <- c(rep("Healthy", 5), rep("Disease", 5))

# ALDEx2差异分析
aldex_result <- aldex(ko_matrix, conditions, mc.samples = 128,
                      test = "t", effect = TRUE, denom = "all")

# 筛选差异KO
sig_ko <- aldex_result[aldex_result$we.eBH < 0.05 & abs(aldex_result$effect) > 0.5, ]

# KEGG pathway富集分析
library(clusterProfiler)
# 将差异KO映射到pathway做富集

# ===== 完整宏基因组组装后分析流水线 =====

# 假设已有各样本的组装结果：sample*_contigs.fa
THREADS=32

# 1. 批量ORF预测
for s in sample1 sample2 sample3 sample4 sample5; do
  prodigal -i ${s}_contigs.fa -a ${s}_prot.faa -d ${s}_gene.fna \
    -p meta -f gff -q -o ${s}.gff &
done
wait

# 2. 合并所有基因（加样本前缀）
> all_genes.fna
for s in sample1 sample2 sample3 sample4 sample5; do
  sed "s/^>/>${s}__/" ${s}_gene.fna >> all_genes.fna
done

# 3. CD-HIT聚类去冗余
cd-hit-est -i all_genes.fna -o nr_genes.fna \
  -c 0.95 -aS 0.9 -n 10 -d 0 -M 64000 -T ${THREADS} -g 1

echo "Non-redundant gene catalog: $(grep -c '>' nr_genes.fna) genes"

# 4. 翻译为蛋白序列
seqkit translate nr_genes.fna > nr_proteins.faa

# 5. 功能注释
emapper.py -i nr_proteins.faa -m diamond --itype proteins \
  --data_dir /data/eggnog_db/ --cpu ${THREADS} \
  --output nr_emapper --output_dir ./annotation/ --override

# 6. KofamScan
exec_annotation -o kofam_results.txt \
  --profile /data/kofam/profiles/ --ko-list /data/kofam/ko_list \
  --cpu ${THREADS} nr_proteins.faa

# 7. 定量（Salmon）
salmon index -t nr_genes.fna -i salmon_idx -p ${THREADS}
for s in sample1 sample2 sample3 sample4 sample5; do
  salmon quant -i salmon_idx -l A \
    -1 ${s}_R1.fq.gz -2 ${s}_R2.fq.gz \
    -o quant/${s} -p 8 --meta --validateMappings
done

# 8. 合并定量结果
python3 -c "
import pandas as pd, os
samples = ['sample1','sample2','sample3','sample4','sample5']
dfs = []
for s in samples:
    df = pd.read_csv(f'quant/{s}/quant.sf', sep='\t', usecols=['Name','TPM'])
    df = df.rename(columns={'TPM': s})
    dfs.append(df.set_index('Name'))
result = pd.concat(dfs, axis=1)
result.to_csv('gene_tpm_matrix.tsv', sep='\t')
print(f'Gene abundance matrix: {result.shape[0]} genes x {result.shape[1]} samples')
"