蛋白组DIA定量分析¶

# === DIA-NN分析流程（命令行版本） ===
# DIA-NN 1.9+ 安装：从 https://github.com/vdemichev/DiaNN 下载

# 1. Library-free模式直接分析DIA数据
diann \
  --f raw_data/sample1.raw \         # 输入原始数据（支持.raw/.mzML/.d）
  --f raw_data/sample2.raw \
  --f raw_data/sample3.raw \
  --lib "" \                          # 空字符串表示library-free模式
  --fasta human_uniprot.fasta \       # 蛋白质序列数据库
  --threads 8 \                       # 使用8个线程
  --out results/report.tsv \          # 输出报告文件
  --out-lib results/spectral_lib.tsv \# 生成的谱库（可用于后续分析）
  --qvalue 0.01 \                     # FDR阈值1%
  --matrices \                        # 输出蛋白/肽段矩阵
  --met-excision \                    # 考虑N端Met切除
  --missed-cleavages 1 \              # 允许1个漏切
  --min-pep-len 7 \                   # 最短肽段长度7个氨基酸
  --max-pep-len 30 \                  # 最长肽段长度30个氨基酸
  --min-pr-charge 2 \                 # 最小前体电荷2+
  --max-pr-charge 4                   # 最大前体电荷4+

# === DIA-NN结果下游分析 ===
import pandas as pd
import numpy as np
from scipy.stats import ttest_ind
from statsmodels.stats.multitest import multipletests

# 1. 读取DIA-NN输出的蛋白定量矩阵
protein_matrix = pd.read_csv("results/report.pg_matrix.tsv", sep="\t", index_col=0)
print(f"鉴定蛋白数: {protein_matrix.shape[0]}")  # 查看鉴定到的蛋白数量
print(f"缺失率: {protein_matrix.isnull().sum().sum() / protein_matrix.size:.2%}")

# 2. 过滤和填补
# 保留至少在70%样本中有定量值的蛋白
min_valid = int(0.7 * protein_matrix.shape[1])  # 至少70%样本有值
protein_filtered = protein_matrix.dropna(thresh=min_valid)  # 过滤

# 对数变换
protein_log = np.log2(protein_filtered)  # log2变换

# 3. 差异蛋白分析
group_disease = ["sample1", "sample2", "sample3"]  # 疾病组样本名
group_control = ["sample4", "sample5", "sample6"]   # 对照组样本名

results = []
for protein in protein_log.index:
    d = protein_log.loc[protein, group_disease].dropna()  # 疾病组数据
    c = protein_log.loc[protein, group_control].dropna()   # 对照组数据
    if len(d) >= 2 and len(c) >= 2:  # 至少2个有效值
        fc = d.mean() - c.mean()  # log2 fold change
        _, pval = ttest_ind(d, c)  # t检验
        results.append({"Protein": protein, "log2FC": fc, "pvalue": pval})

df_results = pd.DataFrame(results)
# FDR校正
df_results["padj"] = multipletests(df_results["pvalue"], method="fdr_bh")[1]
# 筛选差异蛋白
sig = df_results[(df_results["padj"] < 0.05) & (abs(df_results["log2FC"]) > 1)]
print(f"差异蛋白数: {len(sig)}")