脂质组学数据分析¶

# === 脂质组学数据分析流程 ===
import pandas as pd
import numpy as np
import seaborn as sns
import matplotlib.pyplot as plt

# 1. 读取脂质定量数据（已经过MS-DIAL处理）
lipids = pd.read_csv("lipid_quantification.csv", index_col=0)  # 行=样本，列=脂质
groups = pd.read_csv("sample_groups.csv", index_col=0)  # 样本分组信息

# 2. 按脂质大类汇总
def get_lipid_class(name):
    """从脂质名称提取大类，如 PC 34:1 -> PC"""
    return name.split(" ")[0]  # 取空格前的类别名

lipid_classes = [get_lipid_class(col) for col in lipids.columns]  # 提取每列脂质的大类

# 按大类求和
class_totals = pd.DataFrame(index=lipids.index)
for cls in set(lipid_classes):
    cols = [c for c, lc in zip(lipids.columns, lipid_classes) if lc == cls]
    class_totals[cls] = lipids[cols].sum(axis=1)  # 同大类的脂质浓度求和

# 3. 碳链长度与不饱和度分析
def parse_chain(name):
    """解析碳链信息，如 PC 34:1 -> (34, 1)"""
    parts = name.split(" ")
    if len(parts) >= 2 and ":" in parts[1]:
        carbon, double = parts[1].split(":")
        return int(carbon), int(double)
    return None, None

chain_info = pd.DataFrame([parse_chain(col) for col in lipids.columns],
                          columns=["Carbon", "DoubleBond"],
                          index=lipids.columns)

# 4. 绘制脂质类别占比饼图
class_mean = class_totals.mean()  # 各类脂质平均含量
plt.figure(figsize=(8, 8))
plt.pie(class_mean, labels=class_mean.index, autopct='%1.1f%%')  # 饼图
plt.title("Lipid Class Distribution")
plt.savefig("lipid_class_pie.png", dpi=150)

# 5. 差异脂质火山图
from scipy.stats import ttest_ind

disease = lipids.loc[groups["Group"] == "Disease"]  # 疾病组
control = lipids.loc[groups["Group"] == "Control"]  # 对照组

fold_changes = []
p_values = []
for col in lipids.columns:
    fc = np.log2(disease[col].mean() / control[col].mean())  # log2倍数变化
    _, pval = ttest_ind(disease[col], control[col])  # t检验p值
    fold_changes.append(fc)
    p_values.append(pval)

# 绘制火山图
plt.figure(figsize=(10, 8))
plt.scatter(fold_changes, -np.log10(p_values), alpha=0.5, s=20)  # 散点图
plt.axhline(-np.log10(0.05), color="red", linestyle="--")  # p=0.05线
plt.axvline(-1, color="blue", linestyle="--")   # FC=-2线
plt.axvline(1, color="blue", linestyle="--")    # FC=2线
plt.xlabel("log2(Fold Change)")
plt.ylabel("-log10(p-value)")
plt.title("Lipid Volcano Plot")
plt.savefig("lipid_volcano.png", dpi=150)