603 宏转录组分析流程（Metatranscriptomics）¶

========== 实验设计要点 ==========
1. RNA降解极快 → 采样后立即冻存在RNAlater中或液氮速冻
2. 建议每个样本 ≥ 40M paired-end reads（去rRNA后还能剩够用的mRNA reads）
3. 读长推荐 PE150（Illumina NovaSeq）
4. rRNA去除：建议湿实验用rRNA去除试剂盒（Ribo-Zero）+ 干实验用SortMeRNA双重去除
5. 最好同时做宏基因组，这样可以做DNA标准化
6. 生物学重复 ≥ 3

# ========== 步骤1：原始数据质控 ==========

# 用FastQC检查原始数据质量
fastqc \
    raw_data/*_R1.fastq.gz \   # 正向读段
    raw_data/*_R2.fastq.gz \   # 反向读段
    -o fastqc_raw/ \           # 输出目录
    -t 8                       # 使用8个线程

# 用Trimmomatic去接头和低质量碱基
trimmomatic PE \
    -threads 16 \                           # 使用16个线程
    raw_R1.fastq.gz raw_R2.fastq.gz \       # 输入：正反向原始数据
    trimmed_R1.fastq.gz unpaired_R1.fastq.gz \  # 输出：正向（配对+未配对）
    trimmed_R2.fastq.gz unpaired_R2.fastq.gz \  # 输出：反向（配对+未配对）
    ILLUMINACLIP:adapters.fa:2:30:10 \      # 去接头：允许2个错配
    LEADING:3 \                              # 去掉5'端质量<3的碱基
    TRAILING:3 \                             # 去掉3'端质量<3的碱基
    SLIDINGWINDOW:4:20 \                     # 滑动窗口4bp，平均质量<20就切
    MINLEN:50                                # 最短保留50bp的读段

# 也可以用fastp（更快，一步到位）
fastp \
    -i raw_R1.fastq.gz \            # 输入：正向原始数据
    -I raw_R2.fastq.gz \            # 输入：反向原始数据
    -o clean_R1.fastq.gz \          # 输出：正向质控后数据
    -O clean_R2.fastq.gz \          # 输出：反向质控后数据
    --detect_adapter_for_pe \       # 自动检测双端接头
    --cut_front --cut_tail \        # 两端切低质量
    --cut_window_size 4 \           # 滑窗大小4
    --cut_mean_quality 20 \         # 滑窗平均质量阈值
    --length_required 50 \          # 最短50bp
    --thread 16 \                   # 16个线程
    -h fastp_report.html            # HTML质控报告

# ========== 步骤2：去除宿主（人类）序列 ==========
# 宏转录组中可能包含宿主的RNA，需要去除

# 用KneadData一步完成去宿主
kneaddata \
    --input clean_R1.fastq.gz \        # 输入：质控后正向
    --input clean_R2.fastq.gz \        # 输入：质控后反向
    --reference-db /db/human_genome \  # 参考：人类基因组索引
    --output kneaddata_output/ \       # 输出目录
    --threads 16 \                     # 16个线程
    --bypass-trim                      # 跳过修剪（已经用fastp做过了）

# 或者手动用Bowtie2比对去宿主
bowtie2 \
    -x /db/human_genome/hg38 \                # 人类基因组Bowtie2索引
    -1 clean_R1.fastq.gz \                    # 正向读段
    -2 clean_R2.fastq.gz \                    # 反向读段
    --un-conc-gz host_removed_%.fastq.gz \    # 输出：未比对上的（非宿主）序列
    --threads 16 \                            # 16个线程
    -S /dev/null                              # SAM输出丢弃（不需要）

# 检查去宿主效果
echo "原始读段数:"
zcat clean_R1.fastq.gz | wc -l | awk '{print $1/4}'  # 计算原始读段数
echo "去宿主后读段数:"
zcat host_removed_1.fastq.gz | wc -l | awk '{print $1/4}'  # 计算剩余读段数

# ========== 步骤3：去除rRNA序列（宏转录组特有步骤） ==========
# 这是宏转录组与宏基因组最大的区别
# 即使做了湿实验ribo-depletion，仍然会有大量rRNA残留

# 方法1：SortMeRNA（推荐，准确率高）
sortmerna \
    --ref /db/sortmerna/smr_v4.3_default_db.fasta \  # rRNA参考数据库
    --reads host_removed_1.fastq.gz \                 # 输入：正向读段
    --reads host_removed_2.fastq.gz \                 # 输入：反向读段
    --paired_in \                                      # 双端模式：两端都比对上才算rRNA
    --out2 \                                           # 分别输出正反向
    --aligned rRNA_reads \                             # 输出：rRNA读段（扔掉）
    --other non_rRNA_reads \                           # 输出：非rRNA读段（保留！）
    --fastx \                                          # 输出FASTQ格式
    --threads 16 \                                     # 16个线程
    --workdir sortmerna_tmp/                           # 临时文件目录

# 检查rRNA去除效果
echo "去宿主后总读段:"
zcat host_removed_1.fastq.gz | wc -l | awk '{print $1/4}'
echo "rRNA读段数:"
cat rRNA_reads_fwd.fq | wc -l | awk '{print $1/4}'
echo "非rRNA（mRNA）读段数:"
cat non_rRNA_reads_fwd.fq | wc -l | awk '{print $1/4}'
# 通常mRNA只占10-20%，如果>50%说明ribo-depletion效果很好

# 方法2：RiBase（2025年benchmarking显示整体表现更好）
# 但SortMeRNA更成熟，文献引用更多

# ========== 步骤4：用HUMAnN3做功能注释 ==========
# HUMAnN3内部会自动调用MetaPhlAn做物种分类

# 合并双端读段（HUMAnN3需要单文件输入）
cat non_rRNA_reads_fwd.fq non_rRNA_reads_rev.fq > merged_mRNA.fq  # 合并正反向

# 运行HUMAnN3
humann \
    --input merged_mRNA.fq \                          # 输入：去rRNA后的mRNA读段
    --output humann_output/ \                         # 输出目录
    --threads 16 \                                    # 16个线程
    --nucleotide-database /db/chocophlan/ \           # 核酸数据库（ChocoPhlAn）
    --protein-database /db/uniref90/ \                # 蛋白数据库（UniRef90）
    --metaphlan-options="--bowtie2db /db/metaphlan/"  # MetaPhlAn数据库路径

# HUMAnN3会输出3个核心文件：
# 1. *_genefamilies.tsv  — 基因家族丰度（RPK）
# 2. *_pathabundance.tsv — 代谢通路丰度（RPK）
# 3. *_pathcoverage.tsv  — 代谢通路覆盖度（0-1）

# 查看基因家族结果
head -20 humann_output/*_genefamilies.tsv
# 格式：UniRef90_ID|物种名\t丰度值

# ========== 步骤5：多样本合并 ==========

# 合并所有样本的基因家族丰度
humann_join_tables \
    --input humann_output/ \                    # 输入：所有样本的输出目录
    --output merged_genefamilies.tsv \          # 输出：合并的基因家族表
    --file_name genefamilies                    # 匹配文件名

# 合并所有样本的通路丰度
humann_join_tables \
    --input humann_output/ \                    # 输入
    --output merged_pathabundance.tsv \         # 输出：合并的通路丰度表
    --file_name pathabundance                   # 匹配文件名

# ========== 步骤6：标准化 ==========

# CPM标准化（counts per million）
humann_renorm_table \
    --input merged_genefamilies.tsv \           # 输入：合并的基因家族表
    --output genefamilies_cpm.tsv \             # 输出：CPM标准化后
    --units cpm \                               # 标准化单位：CPM
    --update-snames                             # 更新样本名

# 重新分组为GO/KO/EC编号
humann_regroup_table \
    --input genefamilies_cpm.tsv \              # 输入：CPM标准化的基因家族
    --output genefamilies_ko.tsv \              # 输出：KO编号分组
    --groups uniref90_ko                        # 分组方式：UniRef90到KEGG KO

# 重命名为人类可读的名称
humann_rename_table \
    --input genefamilies_ko.tsv \               # 输入：KO编号
    --output genefamilies_ko_named.tsv \        # 输出：带名称的KO
    --names kegg-orthology                      # 使用KEGG命名

# ========== 步骤7：用DNA数据标准化RNA数据 ==========
# 这一步需要同一批样本的宏基因组数据
# RNA/DNA比值反映"真正的表达活性"，消除了基因拷贝数的影响

# 假设已经对宏基因组数据也跑过HUMAnN3
# DNA结果：dna_genefamilies_cpm.tsv
# RNA结果：rna_genefamilies_cpm.tsv

# 用Python做DNA标准化
python3 << 'EOF'
import pandas as pd  # 导入pandas
import numpy as np   # 导入numpy

# 读取RNA和DNA的基因家族丰度
rna = pd.read_csv("rna_genefamilies_cpm.tsv", sep="\t", index_col=0)  # RNA数据
dna = pd.read_csv("dna_genefamilies_cpm.tsv", sep="\t", index_col=0)  # DNA数据

# 对齐：只保留两者都有的基因家族
common_genes = rna.index.intersection(dna.index)  # 找共同的基因家族
rna_aligned = rna.loc[common_genes]  # RNA取共同基因
dna_aligned = dna.loc[common_genes]  # DNA取共同基因

# 计算RNA/DNA比值（加伪计数避免除以0）
pseudocount = 1e-6  # 伪计数
ratio = rna_aligned / (dna_aligned + pseudocount)  # RNA/DNA比值

# 保存结果
ratio.to_csv("rna_dna_ratio.tsv", sep="\t")  # 保存RNA/DNA比值表
print(f"共{len(common_genes)}个基因家族完成DNA标准化")  # 打印统计
EOF

# ========== 步骤8：差异表达分析 ==========
library(DESeq2)  # 差异表达分析包

# 读取基因家族计数矩阵（未标准化的原始计数）
counts <- read.delim("merged_genefamilies.tsv",  # 读取合并的基因家族表
                      row.names = 1,               # 第一列为行名
                      check.names = FALSE)          # 不修改列名

# 去除未分类的行和物种级别的行（只保留群落水平汇总）
counts <- counts[!grepl("\\|", rownames(counts)), ]  # 去掉含"|"的行（物种分层）
counts <- counts[rownames(counts) != "UNMAPPED", ]    # 去掉UNMAPPED
counts <- counts[rownames(counts) != "UNGROUPED", ]   # 去掉UNGROUPED

# 转换为整数（DESeq2要求）
counts <- round(counts)  # 四舍五入为整数

# 准备样本信息
sample_info <- data.frame(
  condition = factor(c(rep("Control", 5), rep("T2D", 5))),  # 分组信息
  row.names = colnames(counts)  # 样本名作为行名
)

# 构建DESeq2对象
dds <- DESeqDataSetFromMatrix(
  countData = counts,         # 计数矩阵
  colData = sample_info,      # 样本信息
  design = ~ condition        # 实验设计：按condition分组
)

# 过滤低表达基因（至少在3个样本中count>10）
keep <- rowSums(counts(dds) >= 10) >= 3  # 过滤条件
dds <- dds[keep, ]                        # 应用过滤

# 运行差异表达分析
dds <- DESeq(dds)  # DESeq2分析

# 提取结果
res <- results(dds,                           # 从DESeq对象提取结果
               contrast = c("condition", "T2D", "Control"),  # 比较：T2D vs Control
               alpha = 0.05)                  # 显著性阈值

# 查看显著差异基因
sig_genes <- subset(res, padj < 0.05 & abs(log2FoldChange) > 1)  # 筛选显著差异
sig_genes <- sig_genes[order(sig_genes$padj), ]  # 按P值排序
cat("显著差异基因数:", nrow(sig_genes), "\n")     # 打印数量

# 火山图
library(EnhancedVolcano)  # 火山图专用包
EnhancedVolcano(res,
                lab = rownames(res),       # 标签
                x = 'log2FoldChange',      # X轴：log2倍数变化
                y = 'padj',                # Y轴：校正后P值
                title = 'T2D vs Control',  # 标题
                pCutoff = 0.05,            # P值阈值线
                FCcutoff = 1)              # FC阈值线

# ========== 可选：基于组装的分析（适合新物种多的样本） ==========

# 用MEGAHIT组装
megahit \
    -1 non_rRNA_R1.fq \     # 输入：去rRNA后正向
    -2 non_rRNA_R2.fq \     # 输入：去rRNA后反向
    -o megahit_assembly/ \  # 输出目录
    -t 32 \                 # 32个线程
    --min-contig-len 500    # 最短contig 500bp

# 基因预测（Prodigal）
prodigal \
    -i megahit_assembly/final.contigs.fa \  # 输入：组装结果
    -o genes.gff \                           # 输出：GFF格式基因注释
    -a proteins.faa \                        # 输出：蛋白序列
    -p meta \                                # 宏基因组模式
    -f gff                                   # 输出GFF格式

# 功能注释（eggNOG-mapper）
emapper.py \
    -i proteins.faa \                 # 输入：蛋白序列
    -o eggnog_results \               # 输出前缀
    --cpu 16 \                        # 16个CPU
    -m diamond \                      # 使用DIAMOND比对
    --data_dir /db/eggnog/            # eggNOG数据库路径

# 定量：将读段比对回组装的基因
salmon quant \
    -t genes.fa \                    # 参考：预测的基因序列
    -l A \                           # 自动检测文库类型
    -1 non_rRNA_R1.fq \             # 正向读段
    -2 non_rRNA_R2.fq \             # 反向读段
    -o salmon_output/ \             # 输出目录
    --threads 16                    # 16个线程

========== 标准宏转录组分析流程 ==========
原始数据 → FastQC/fastp质控 → 去宿主(Bowtie2/KneadData)
    → 去rRNA(SortMeRNA) → HUMAnN3(功能注释)
    → 标准化(CPM/DNA标准化) → 差异分析(DESeq2)

========== 核心工具一览 ==========
质控     ：fastp / Trimmomatic / FastQC
去宿主   ：KneadData / Bowtie2
去rRNA   ：SortMeRNA（最常用）/ RiBase（2025新推荐）
功能注释 ：HUMAnN3（assembly-free推荐）
组装     ：MEGAHIT / Trinity
基因预测 ：Prodigal / FragGeneScan
物种分类 ：MetaPhlAn4（HUMAnN内置）
差异分析 ：DESeq2 / edgeR / ALDEx2
整合流程 ：SAMSA2 / MetaPro / MetaTrans

========== 关键数据库 ==========
ChocoPhlAn ：核酸参考数据库（HUMAnN用）
UniRef90   ：蛋白参考数据库（HUMAnN用）
MetaPhlAn DB：物种标记基因数据库
SILVA      ：rRNA参考（SortMeRNA用）
KEGG / MetaCyc：通路注释

========== 测序建议 ==========
读长     ：PE150
深度     ：≥40M reads/样本（去rRNA前）
rRNA去除 ：湿实验(Ribo-Zero) + 干实验(SortMeRNA)
存储     ：RNAlater或液氮速冻
重复     ：≥3个生物学重复