一、流程说明
1.批量统计多个榧属物种高质量蛋白序列长度，通过Mann-Whitney U检验与Cliff's Delta效应量，对比各物种与参考物种香榧的蛋白长度差异
2.绘制小提琴箱线图，横轴榧属不同物种，纵轴蛋白序列长度log10aa
[Cliff's Delta效应量]
非参数的效应量指标，用来比较两组数据分布差异到底有多大
使用条件：
1.两独立样本：Cliff's delta 主要用于两个独立组的比较（例如实验组和对照组），不能用于配对样本或单一样本，两组样本量可以不等，也没有最小样本量的严格限制
2.不要求正态分布，不要求方差齐性：可以是连续数据（如反应时、分数）、顺序数据（如等级排名、里克特量表），甚至是数据中存在极端异常值
3.Cliff's delta允许两组分布形状完全不同，比如一组严重偏态，一组均匀分布，当两组的分布形态差异过大时，虽然完全可以计算，但在解读时需要更谨慎。一个显著的delta值可能同时反映了“位置的偏移”和“形状的差异”，不一定单纯是总体水平的移动
4.当数据中有很多相同的数值时，Cliff's delta依然可用，但是此时效应量的绝对值可能会被轻微低估
分级：
二、环境配置
conda install -c conda-forge pandas numpy matplotlib seaborn scipy
conda install -c bioconda biopython三、核心代码
from pathlib import Path
from Bio import SeqIO
import pandas as pd
import seaborn as sns
import matplotlib.pyplot as plt
from scipy.stats import mannwhitneyu
FINAL_DIR = Path("/data/result/species_proteome")
target_species = ['fargesii', 'nucifera', 'taxifolia', 'jackii', 'Tgra']
OUTPUT_PDF = "/data/result/Protein_Length_Distribution.pdf"
OUTPUT_PNG = "/data/result/Protein_Length_Distribution.png"
OUTPUT_CSV = "/data/result/Effect_Size_Results.csv"
# 蛋白长度数据收集
def collect_length_data():
data = []
for sp in target_species:
pep_file = FINAL_DIR / sp / f"{sp}.final.filtered.pep.fa"
if not pep_file.exists():
print(f"[警告] 缺失文件：{pep_file}")
continue
for record in SeqIO.parse(pep_file, "fasta"):
data.append({"Species": sp, "Length (aa)": len(record.seq)})
return pd.DataFrame(data)
# 统计检验，U检验与效应量
def run_stat_test(df):
tgra_lengths = df[df['Species'] == 'Tgra']['Length (aa)'].values
results = []
for sp in ['fargesii', 'nucifera', 'taxifolia', 'jackii']:
sp_lengths = df[df['Species'] == sp]['Length (aa)'].values
# 曼惠特尼U检验
stat, p_val = mannwhitneyu(sp_lengths, tgra_lengths, alternative='two-sided')
# 计算Cliff's Delta效应量
cliffs_d = (2.0 * stat) / (len(sp_lengths) * len(tgra_lengths)) - 1.0
abs_d = abs(cliffs_d)
# 效应量分级
level = "Negligible" if abs_d < 0.147 else "Small" if abs_d < 0.330 else "Medium" if abs_d < 0.474 else "Large"
results.append({"Species": sp, "p_value": p_val, "Cliff_Delta": round(cliffs_d,4), "Effect_Level": level})
res_df = pd.DataFrame(results)
res_df.to_csv(OUTPUT_CSV, index=False)
return dict(zip(res_df['Species'], res_df.to_dict('records')))
#可视化绘图，小提琴箱线图
def plot_distribution(df, effect_dict):
plt.figure(figsize=(11, 6.5))
# 自定义色板
palette = {'fargesii':'#9499C8','nucifera':'#4E79A7','jackii':'#76B7B2','taxifolia':'#59A14F','Tgra':'#FC8D62'}
# 绘制小提琴图
sns.violinplot(x="Species", y="Length (aa)", data=df, palette=palette, cut=0, inner=None)
# 叠加箱线图
sns.boxplot(x="Species", y="Length (aa)", data=df, width=0.08, color="white", showcaps=False, showfliers=False)
# 对数坐标轴
plt.yscale('log')
# 标注中位数、效应量
medians = df.groupby('Species')['Length (aa)'].median()
for i, sp in enumerate(target_species):
plt.text(i-0.06, medians[sp], f"{medians[sp]:.0f}", va='center', ha='right')
# 图表美化
plt.title("Protein Length Distribution Comparison (Torreya & T. grandis)", fontweight='bold')
plt.xlabel("Species", fontweight='bold')
plt.ylabel("Protein Length (Log10 Scale)", fontweight='bold')
sns.despine()
plt.tight_layout()
# 保存图表
plt.savefig(OUTPUT_PDF, dpi=300)
plt.savefig(OUTPUT_PNG, dpi=300)
# 主程序
if __name__ == "__main__":
df = collect_length_data()
if not df.empty:
print(df.groupby("Species")["Length (aa)"].describe()[['count','mean','50%']].astype(int))
effect_dict = run_stat_test(df)
plot_distribution(df, effect_dict)
print("[完成] 数据统计与绘图结束！")四、结果输出
1.统计结果文件
/data/result/Effect_Size_Results.csv，包含榧属各物种与香榧对比的p 值,Cliff's Delta 值,效应量等级
2.分析图
PDF 格式：/data/result/Protein_Length_Distribution.pdf
PNG 格式：/data/result/Protein_Length_Distribution.png