人类各染色体缺失/重复统计柱状图(matplotlib)

该博客介绍了如何利用Python的matplotlib库,绘制人类24条染色体的缺失和重复统计柱状图。首先,定义了读取数据文件的函数,然后合并绘图所需的数据,接着绘制柱状图,并最终展示结果。

绘制人类各条染色体缺失/重复统计柱状图。

读取数据文件函数

import pandas as pd
import matplotlib as plt

def read_file(file_path: str):
    if file_path.endswith(('.tsv', '.txt')):
        # 缺失值填充为NA
        return pd.read_csv(file_path, sep='\t').fillna('NA') if os.path.exists(file_path) else pd.DataFrame()
    elif file_path.endswith(('.xls', '.xlsx')):
        return pd.read_excel(file_path).fillna('NA') if os.path.exists(file_path) else pd.DataFrame()
    else:
        raise Exception("ERROR FILE FORMAT")

合并绘图数据

def merge_statistics(outdir='./'):
		# 读取文件
    dataframe1=read_file(file_path=outdir + "data1.tsv")
    dataframe2=read_file(
import pandas as pd import numpy as np from sklearn.model_selection import train_test_split from sklearn.metrics import classification_report, roc_auc_score import xgboost as xgb import shap import matplotlib.pyplot as plt import warnings # 忽略警告 warnings.filterwarnings('ignore') # 1. 数据读取与预处理 =========================================================== # 读取Excel文件(第二个工作表) df = pd.read_excel('附件.xlsx', sheet_name=1) print(f"数据维度: {df.shape}") print("列名:", df.columns.tolist()) # 2. 特征工程 ================================================================= # 2.1 构造标签(异常=1, 正常=0) df['label'] = df['染色体的非整倍体'].notnull().astype(int) # 2.2 特征构造 # GC偏差 df['GC_bias'] = (df['GC含量'] - 0.4).abs() # 读段比例 df['reads_ratio'] = df['唯一比对的读段数'] / df['原始读段数'] # Z值相关特征 z_columns = ['13号染色体的Z值', '18号染色体的Z值', '21号染色体的Z值'] df['max_Z'] = df[z_columns].abs().max(axis=1) df['sum_Z_sq'] = df[z_columns].apply(lambda x: np.sum(np.square(x)), axis=1) # X染色体偏移 df['ZX_offset'] = df['X染色体的Z值'].abs() # 3. 特征/标签准备 ============================================================= # 特征列表(包含原始特征和构造特征) features = [ '13号染色体的Z值', '18号染色体的Z值', '21号染色体的Z值', 'X染色体的Z值', 'GC_bias', 'reads_ratio', '孕妇BMI', 'max_Z', 'sum_Z_sq', 'ZX_offset' ] # 处理缺失值(用0填充) X = df[features].fillna(0) y = df['label'] # 4. 数据划分 ================================================================= X_train, X_test, y_train, y_test = train_test_split( X, y, test_size=0.25, random_state=42, stratify=y # 分层抽样保持类别比例 ) print(f"\n训练集: {X_train.shape[0]} 样本, 测试集: {X_test.shape[0]} 样本") print(f"阳性样本比例: 训练集 {y_train.mean():.2%}, 测试集 {y_test.mean():.2%}") # 5. 模型训练 ================================================================= # 计算类别权重处理不平衡数据 scale_pos_weight = (y_train == 0).sum() / (y_train == 1).sum() model = xgb.XGBClassifier( n_estimators=300, max_depth=4, learning_rate=0.1, subsample=0.8, colsample_bytree=0.8, objective='binary:logistic', scale_pos_weight=scale_pos_weight, random_state=42, eval_metric='auc' ) # 训练模型 model.fit(X_train, y_train) # 6. 模型评估 ================================================================= # 测试集预测 y_pred = model.predict(X_test) y_proba = model.predict_proba(X_test)[:, 1] # 输出评估报告 print("\n" + "="*50) print("模型评估报告:") print("="*50) print(classification_report(y_test, y_pred)) print(f"AUC = {roc_auc_score(y_test, y_proba):.4f}") # 7. SHAP 解释 ================================================================ # 初始化解释器 explainer = shap.TreeExplainer(model) # 计算SHAP值(使用部分样本提高速度) sample_indices = np.random.choice(X_test.index, min(500, len(X_test)), replace=False) X_sample = X_test.loc[sample_indices] shap_values = explainer.shap_values(X_sample) # 7.1 特征重要性摘要图 plt.figure(figsize=(10, 6)) shap.summary_plot(shap_values, X_sample, plot_type="bar", show=False) plt.title("Feature importance(SHAP)", fontsize=14) plt.tight_layout() plt.savefig('SHAP_特征重要性.png', dpi=300) # 7.2 单样本force plot plt.figure(figsize=(12, 4)) abn_samples = X_test[y_test == 1] if not abn_samples.empty: sample_idx = abn_samples.index[0] sample_data = X_test.loc[[sample_idx]] sample_shap = explainer.shap_values(sample_data) shap.force_plot( explainer.expected_value, sample_shap[0], sample_data.iloc[0], matplotlib=True, show=False ) plt.title(f"SHAP explanation of sample {sample_idx} (anomalous sample)", fontsize=12) plt.tight_layout() plt.savefig('SHAP_异常样本解释.png', dpi=300) # 8. 生成可解释报告 =========================================================== def get_shap_contributions(row): """获取单样本的SHAP贡献值""" idx_val = row.name if idx_val in X_sample.index: idx_pos = list(X_sample.index).index(idx_val) return dict(zip(features, shap_values[idx_pos])) return {col: np.nan for col in features} # 创建报告 report = pd.DataFrame() report['样本ID'] = X_test.index report = report.set_index('样本ID') shap_df = report.apply(get_shap_contributions, axis=1, result_type='expand') # 合并结果 final_report = pd.concat([ X_test, shap_df.add_prefix('SHAP_'), pd.Series(model.predict(X_test), index=X_test.index, name='预测标签'), pd.Series(y_proba, index=X_test.index, name='预测概率'), y_test.rename('真实标签') ], axis=1) # 保存报告 final_report.to_csv('女胎异常判定_可解释报告.csv', encoding='utf-8-sig') print("\n已生成 '女胎异常判定_可解释报告.csv'") print("="*50) print("处理完成!") 生成图后存在乱码,怎么进行修改
09-08
import pandas as pd import numpy as np import matplotlib.pyplot as plt import seaborn as sns from scipy import stats from sklearn.cluster import KMeans from sklearn.linear_model import LinearRegression from sklearn.preprocessing import PolynomialFeatures from sklearn.pipeline import Pipeline import warnings warnings.filterwarnings('ignore') # 读取清洗后的数据 df = pd.read_csv('cleaned_data.csv') # 确保只处理男胎数据:Y染色体浓度非空 df_male = df[df['Y染色体浓度'].notna()].copy() print(f"男胎数据量: {len(df_male)}") # 计算孕周数(如果尚未计算) def extract_weeks(s): if isinstance(s, str): if 'w' in s: parts = s.split('w') weeks = float(parts[0]) if '+' in parts[1]: days = float(parts[1].split('+')[1].replace('d', '')) weeks += days / 7 return weeks return np.nan if '孕周数' not in df_male.columns: df_male['孕周数'] = df_male['检测孕周'].apply(extract_weeks) # 移除孕周数或Y染色体浓度为NaN的行 df_male = df_male.dropna(subset=['孕周数', 'Y染色体浓度', '孕妇BMI']) # 定义BMI分组边界(初始分组基于临床经验) bmi_bins = [20, 28, 32, 36, 40, np.inf] bmi_labels = ['20-28', '28-32', '32-36', '36-40', '40+'] df_male['BMI分组'] = pd.cut(df_male['孕妇BMI'], bins=bmi_bins, labels=bmi_labels, right=False) # 检查每组数据量 print("各BMI分组数据量:") print(df_male['BMI分组'].value_counts()) # 对于每个BMI分组,找到Y染色体浓度达到4%的孕周数 # 由于每个孕妇可能有多次检测,我们首先对每个孕妇找到浓度达标的最早孕周 # 分组按孕妇代码和BMI分组 df_male['达标'] = df_male['Y染色体浓度'] >= 4.0 # 对于每个孕妇,找到达标的最早孕周 earliest_weeks = df_male[df_male['达标']].groupby('孕妇代码')['孕周数'].min().reset_index() earliest_weeks.rename(columns={'孕周数': '达标孕周'}, inplace=True) # 合并回原数据,获取每个孕妇的BMI分组 patient_bmi = df_male.groupby('孕妇代码')['孕妇BMI'].first().reset_index() patient_bmi['BMI分组'] = pd.cut(patient_bmi['孕妇BMI'], bins=bmi_bins, labels=bmi_labels, right=False) earliest_weeks = earliest_weeks.merge(patient_bmi, on='孕妇代码', how='left') # 对于每个BMI分组,计算达标孕周的分布 grouped = earliest_weeks.groupby('BMI分组')['达标孕周'] descriptive_stats = grouped.describe(percentiles=[0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95]) print("各BMI分组达标孕周的描述统计:") print(descriptive_stats) # 确定最佳NIPT时点:对于每组,选择第95百分位数的达标孕周,以确保95%的孕妇在该时点已达标 best_timing = grouped.quantile(0.95).reset_index() best_timing.rename(columns={'达标孕周': '最佳NIPT时点()'}, inplace=True) print("各BMI分组的最佳NIPT时点(第95百分位数):") print(best_timing) # 可视化达标孕周分布 plt.figure(figsize=(10, 6)) sns.boxplot(x='BMI分组', y='达标孕周', data=earliest_weeks) plt.axhline(y=12, color='g', linestyle='--', label='12周(早期风险低)') plt.axhline(y=27, color='r', linestyle='--', label='27周(中期风险高)') plt.title('各BMI分组Y染色体浓度达标孕周分布') plt.ylabel('孕周(周)') plt.legend() plt.savefig('各BMI分组达标孕周分布.png') plt.show() # 风险分析:比较最佳时点与风险周数 best_timing['风险等级'] = best_timing['最佳NIPT时点()'].apply( lambda x: '低风险' if x <= 12 else '高风险' if x <= 27 else '极高风险' ) print("最佳NIPT时点的风险等级:") print(best_timing[['BMI分组', '最佳NIPT时点()', '风险等级']]) # 检测误差分析 # 假设Y染色体浓度测量误差为±0.5%(根据临床典型误差),孕周估计误差为±0.5周 # 使用蒙特卡洛模拟评估误差对达标时间估计的影响 np.random.seed(42) n_simulations = 1000 error_results = [] for bmi_group in bmi_labels: group_data = earliest_weeks[earliest_weeks['BMI分组'] == bmi_group] if len(group_data) == 0: continue # 模拟误差 simulated_weeks = [] for _, row in group_data.iterrows(): base_week = row['达标孕周'] # 模拟孕周误差和浓度误差(浓度误差可能影响达标时间,但这里我们直接模拟达标孕周误差) # 由于达标孕周是基于浓度≥4%定义的,误差会同时影响浓度和孕周 # 简化:假设达标孕周估计有误差,误差分布为正态(均值0,标准差0.5周) simulated_week = np.random.normal(base_week, 0.5, n_simulations) simulated_weeks.extend(simulated_week) simulated_weeks = np.array(simulated_weeks) # 计算模拟后的第95百分位数 sim_95perc = np.percentile(simulated_weeks, 95) error_results.append({ 'BMI分组': bmi_group, '原始最佳时点': best_timing[best_timing['BMI分组'] == bmi_group]['最佳NIPT时点()'].iloc[0], '模拟最佳时点': sim_95perc, '误差': sim_95perc - best_timing[best_timing['BMI分组'] == bmi_group]['最佳NIPT时点()'].iloc[0] }) error_df = pd.DataFrame(error_results) print("检测误差对最佳NIPT时点的影响:") print(error_df) # 可视化误差影响 plt.figure(figsize=(10, 6)) x_pos = np.arange(len(bmi_labels)) plt.bar(x_pos, error_df['原始最佳时点'], width=0.4, label='原始最佳时点') plt.bar(x_pos + 0.4, error_df['模拟最佳时点'], width=0.4, label='模拟最佳时点(含误差)') plt.xticks(x_pos + 0.2, bmi_labels) plt.xlabel('BMI分组') plt.ylabel('最佳NIPT时点(周)') plt.title('检测误差对最佳NIPT时点的影响') plt.legend() plt.savefig('检测误差影响.png') plt.show() # 输出最终分组和最佳时点 print("\n最终BMI分组和最佳NIPT时点:") result_df = best_timing.copy() result_df['BMI区间'] = result_df['BMI分组'].apply( lambda x: '[20,28)' if x == '20-28' else '[28,32)' if x == '28-32' else '[32,36)' if x == '32-36' else '[36,40)' if x == '36-40' else '40+') result_df = result_df[['BMI分组', 'BMI区间', '最佳NIPT时点()', '风险等级']] print(result_df) 这段代码在运行后出现以下问题 C:\Study\anaconda3\envs\MCM\python.exe C:\Study\MCM\2.1.py 男胎数据量: 1082 各BMI分组数据量: BMI分组 28-32 539 32-36 412 36-40 93 20-28 19 40+ 18 Name: count, dtype: int64 各BMI分组达标孕周的描述统计: Empty DataFrame Columns: [count, mean, std, min, 5%, 10%, 25%, 50%, 75%, 90%, 95%, max] Index: [] 各BMI分组的最佳NIPT时点(第95百分位数): BMI分组 最佳NIPT时点() 0 20-28 NaN 1 28-32 NaN 2 32-36 NaN 3 36-40 NaN 4 40+ NaN 最佳NIPT时点的风险等级: BMI分组 最佳NIPT时点() 风险等级 0 20-28 NaN 极高风险 1 28-32 NaN 极高风险 2 32-36 NaN 极高风险 3 36-40 NaN 极高风险 4 40+ NaN 极高风险 检测误差对最佳NIPT时点的影响: Empty DataFrame Columns: [] Index: [] Traceback (most recent call last): File "C:\Study\MCM\2.1.py", line 126, in <module> plt.bar(x_pos, error_df['原始最佳时点'], width=0.4, label='原始最佳时点') File "C:\Study\anaconda3\envs\MCM\lib\site-packages\pandas\core\frame.py", line 4102, in __getitem__ indexer = self.columns.get_loc(key) File "C:\Study\anaconda3\envs\MCM\lib\site-packages\pandas\core\indexes\range.py", line 417, in get_loc raise KeyError(key) KeyError: '原始最佳时点' 进程已结束,退出代码为 1
09-07
为以下代码生成可视化图表:import pandas as pd from sklearn.model_selection import train_test_split, GridSearchCV from sklearn.ensemble import RandomForestClassifier from sklearn.metrics import accuracy_score, classification_report from sklearn.preprocessing import LabelEncoder from sklearn.feature_selection import SelectKBest, f_classif from imblearn.over_sampling import SMOTE # 读取文件 excel_file = pd.ExcelFile(r"C:\Users\admin\Desktop\fj.xlsx") # 获取指定工作表中的数据 df = excel_file.parse('女胎检测数据') # 删除染色体的非整倍体列中存在缺失值的行 df = df.dropna(subset=['染色体的非整倍体']) # 选取特征变量和目标变量 X = df[['13号染色体的Z值', '18号染色体的Z值', '21号染色体的Z值', 'X染色体的Z值', '13号染色体的GC含量', '18号染色体的GC含量', '21号染色体的GC含量', '孕妇BMI', '在参考基因组上比对的比例', '重复读段的比例', '唯一比对的读段数', '被过滤掉读段数的比例', 'GC含量']] y = df['染色体的非整倍体'] # 对目标变量进行编码 encoder = LabelEncoder() y = encoder.fit_transform(y) # 特征选择,选择 8 个最佳特征 selector = SelectKBest(score_func=f_classif, k=8) X = selector.fit_transform(X, y) # 使用 SMOTE 进行过采样,设置 k_neighbors 为 1 smote = SMOTE(random_state=42, k_neighbors=1) X_resampled, y_resampled = smote.fit_resample(X, y) # 划分训练集和测试集 X_train, X_test, y_train, y_test = train_test_split(X_resampled, y_resampled, test_size=0.2, random_state=42) # 定义要调整的参数范围 param_grid = { 'n_estimators': [100, 200, 300], 'max_depth': [None, 10, 20], 'min_samples_split': [2, 5, 10] } # 创建随机森林模型和网格搜索对象 rf = RandomForestClassifier(random_state=42) grid_search = GridSearchCV(estimator=rf, param_grid=param_grid, cv=3, scoring='accuracy') # 进行网格搜索 grid_search.fit(X_train, y_train) # 获取最优模型 best_rf = grid_search.best_estimator_ # 在测试集上进行预测 y_pred = best_rf.predict(X_test) # 评估模型性能 accuracy = accuracy_score(y_test, y_pred) print(f'模型准确率:{accuracy}') print('分类报告:') print(classification_report(y_test, y_pred)) # 这里假设原来是写错的地方,修改为 print print(y_test, y_pred)
09-08
评论
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

生信与基因组学

每一份鼓励是我坚持下去动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值