一：赛题背景
生物学年龄评价是一种通过测量和分析生物体内特定指标或生理过程的状态，以评估个体的生理年龄和健康状况的方法。与传统的日历年龄相比，生物学年龄可以提供更准确的健康评估和疾病风险预测。随着 AI 技术的不断发展，计算科学与生命科学的整合与发展将为健康管理应用开发、衰老机制研究、抗衰老药物研发提供前沿的思路与方法。基于上述背景，我们举办基于甲基化测量数据预测生物学年龄的比赛。
这项比赛将结合计算科学和生命科学的发展，为健康管理应用的开发、衰老机制的研究以及抗衰老药物的研发提供前沿的思路和方法。传统的日历年龄只提供了时间上的度量，而生物学年龄则更加关注个体的生理状态和功能。通过测量和分析生物体内的特定指标或生理过程的状态，我们能够更准确地评估个体的生理年龄和健康状况，这对于健康管理和疾病风险预测至关重要。参赛者可以利用甲基化测量数据来预测生物学年龄。本次比赛将为研究人员和开发人员提供机会，合作解决生物学年龄评估中的挑战，并推动相关领域的进步。最终的目标是提供更精确、可靠和个性化的生物学年龄评估工具，为人们的健康管理和抗衰老策略提供有力支持。
二：背景知识
DNA甲基化是一种常见的表观遗传学现象，它属于基因组核酸序列的功能性修饰，可以在不改变 DNA 序列的情况下影响基因的表达，其基本特点是在DNA分子中添加一个甲基基团，会对基因表达和细胞功能产生重要影响。甲基化通常发生在CpG位点，即一个C碱基与紧邻的G碱基相连的位置。该过程由DNA甲基转移酶催化，将甲基选择性地添加到特定碱基上。
随着个体年龄的增长，DNA甲基化的模式和水平也会发生变化。这种变化可以反映细胞和组织的衰老过程，在整体上与个体的年龄相关联，并在一定程度上解释了为何随着年龄的增长，个体更容易患上与年龄相关的疾病，是预测衰老、疾病与死亡的个体化生物时钟。
三：比赛任务
本次比赛提供健康人和年龄相关疾病患者的甲基化数据，期待选手通过分析甲基化数据的模式和特征建立预测模型，可以根据某个人的甲基化数据来预测其生物学年龄。
初赛
任务：构建一个生物学时钟能够评价所提供样本的生物学年龄，
对于健康样本，生物学年龄尽可能接近真实年龄。
对于患病样本，生物学年龄比真实年龄更大。
复赛
在不提供任何额外数据的前提下，
对于健康样本，生物学年龄尽可能接近真实年龄。
对于患病样本，生物学年龄比真实年龄更大。
预测测试集中哪些样本患有Alzheimer’s disease。
四：数据集介绍
初赛
公开数据包含10296个样本，其中7833个样本为健康样本。每一个样本提供485512个位点的甲基化数据、年龄与患病情况。抽取80%作为训练样本，20%作为测试样本。
以训练集为例，一共包括8233样本，其中健康样本6266个，其余为患病样本，共涉及Alzheimer’s disease，schizophrenia，Parkinson’s disease，rheumatoid arthritis，stroke，Huntington’s disease，Graves’ disease，type 2 diabetes，Sjogren’s syndrome等类型。
数据集格式如下:
- trainmap.csv : 每一行代表一个样本
- traindata.csv : 大小为(8233, 485512)，每一列代表一个样本，每一行代表一个甲基化位点。
测试集testdata.csv 格式与训练集traindata.csv 格式一致。testmap.csv仅提供sample_id与gender。
复赛
复赛测试集和初赛测试集格式一致，包含150个样本的甲基化位点数据以及性别信息。这里仅包含type 2 diabetes， stroke和 Alzheimer’s disease。
五：评价指标
六：baseline代码明细
- 导入可能用到的包
#导入模块
import numpy as np
import pandas as pd
from collections import defaultdict, Counter
# defaultdict 允许我们在定义字典时给所有不存在的key设置默认值，这样当取不存在的key时，就不会报错
# Counter类的目的是用来跟踪值出现的次数。它是一个无序的容器类型，以字典的键值对形式存储，
# 其中元素作为key，其计数作为value。计数值可以是任意的Interger（包括0和负数）
import xgboost as xgb
import lightgbm as lgb
from catboost import CatBoostRegressor
from sklearn.model_selection import StratifiedKFold, KFold, GroupKFold
# StratifiedKFlod：分层采样，训练集与测试集中各类别样本的比列与原始数据中相同：(分类问题适用)
# KFlod：分层采样，将数据分成训练集和测试集，不考虑训练集与测试集中各类别数据是否相同；(回归问题适用)
# GroupKFold：特征列如果也代表类别，我们希望将这个特征列相同类别划成一组,就像df.groupby一样意思。
from sklearn.metrics import mean_squared_log_error,mean_absolute_error
import sys,os,gc,argparse,warnings
# argparse是一个python模块，用途是：命令行选项、参数和子命令的解释
warnings.filterwarnings('ignore')
- 编写压缩数据内存函数，返回压缩后的Datafram
#压缩内存函数，本数据过大，算力有限情况下运行该函数可以显著降低内存（50% - 70%），但运行时间极长
# 该方法的思路为：根据字段的最大值和最小值来匹配合适的数据类型。比如df中某一列变量（即字段）的最大值为1，最小值为0，
# 那么我们将该字段的数据类型转换为占用内存更少的int16来存储这一变量。
def reduce_mem_usage(df, verbose=True):
# 表示是否要输出冗余结果（内存占用减少了多少）
numerics = ['int16', 'int32', 'int64', 'float16', 'float32', 'float64']
start_mem = df.memory_usage().sum() / 1024**2 # start内存占用情况
for col in df.columns:
col_type = df[col].dtypes
if col_type in numerics:
c_min = df[col].min()
c_max = df[col].max()
if str(col_type)[:3] == 'int':
if c_min > np.iinfo(np.int8).min and c_max < np.iinfo(np.int8).max:
df[col] = df[col].astype(np.int8)
elif c_min > np.iinfo(np.int16).min and c_max < np.iinfo(np.int16).max:
df[col] = df[col].astype(np.int16)
elif c_min > np.iinfo(np.int32).min and c_max < np.iinfo(np.int32).max:
df[col] = df[col].astype(np.int32)
elif c_min > np.iinfo(np.int64).min and c_max < np.iinfo(np.int64).max:
df[col] = df[col].astype(np.int64)
else:
if c_min > np.finfo(np.float16).min and c_max < np.finfo(np.float16).max:
df[col] = df[col].astype(np.float16)
elif c_min > np.finfo(np.float32).min and c_max < np.finfo(np.float32).max:
df[col] = df[col].astype(np.float32)
else:
df[col] = df[col].astype(np.float64)
end_mem = df.memory_usage().sum() / 1024**2
if verbose: print('Mem. usage decreased to {:5.2f} Mb ({:.1f}% reduction)'.format(end_mem, 100 * (start_mem - end_mem) / start_mem))
return df
- 加载数据并将数据转化为常见格式，重写索引
- 这里读取了50000个甲基化位点，比baseline的10000个结果要好
path = 'ai4bio'
# baseline取10000条数据
traindata = pd.read_csv(f'{path}/traindata.csv', nrows=50000)
# 读取50000个甲基化位点数据
trainmap = pd.read_csv(f'{path}/trainmap.csv')
testdata = pd.read_csv(f'{path}/ai4bio_testset_final/testdata.csv', nrows=50000)
testmap = pd.read_csv(f'{path}/ai4bio_testset_final/testmap.csv')
#可选择是否运行压缩内存函数
# traindata = reduce_mem_usage(traindata)
# testdata = reduce_mem_usage(testdata)
traindata = traindata.set_index('cpgsite')
traindata = traindata.T
traindata = traindata.reset_index()
traindata = traindata.rename(columns={'index':'sample_id'})
traindata.columns = ['sample_id'] + [i for i in range(50000)]
# 保留sample_id唯一标识，甲基化位点数值编码
traindata.to_pickle(f'{path}/traindata_3.pkl')
# 生成pickle文件对数据进行永久储存
testdata = testdata.set_index('cpgsite')
testdata = testdata.T
testdata = testdata.reset_index()
testdata = testdata.rename(columns={'index':'sample_id'})
testdata.columns = ['sample_id'] + [i for i in range(50000)]
testdata.to_pickle(f'{path}/testdata_3.pkl')
- 读取数据
trainmap = pd.read_csv(f'{path}/trainmap.csv')
testmap = pd.read_csv(f'{path}/ai4bio_testset_final/testmap.csv')
traindata = pd.read_pickle(f'{path}/traindata_3.pkl')
testdata = pd.read_pickle(f'{path}/testdata_3.pkl')
- 查看数据形式
trainmap.head()
traindata.head()
traindata.info()
# 可以看到所有的训练样本已加载进来
for i in range(20):
null_cnt = traindata[i].isnull().sum() / traindata.shape[0]
print(f'特征{i}，对应的缺失率为{null_cnt}')
# 可以看到大多数特征的缺失率都较高
traindata[[i for i in range(1000)]].corr()
# 这里选取1000条特征来做相关性分析
testmap
# 拼接数据集
traindata = traindata.merge(trainmap[['sample_id', 'age', 'gender', 'sample_type', 'disease']],on='sample_id',how='left')
testdata = testdata.merge(testmap[['sample_id', 'gender']],on='sample_id',how='left')
traindata.head()
- 数据编码
disease_mapping = {
'control': 0,
"Alzheimer's disease": 1,
"Graves' disease": 2,
"Huntington's disease": 3,
"Parkinson's disease": 4,
'rheumatoid arthritis': 5,
'schizophrenia': 6,
"Sjogren's syndrome": 7,
'stroke': 8,
'type 2 diabetes': 9
}
sample_type_mapping = {'control': 0, 'disease tissue': 1}
gender_mapping = {'F': 0, 'M': 1}
traindata['disease'] = traindata['disease'].map(disease_mapping)
traindata['sample_type'] = traindata['sample_type'].map(sample_type_mapping)
traindata['gender'] = traindata['gender'].map(gender_mapping)
testdata['gender'] = testdata['gender'].map(gender_mapping)
- 加入描述性特征
traindata['max'] = traindata[[i for i in range(10000)]].max(axis=1)
traindata['min'] = traindata[[i for i in range(10000)]].min(axis=1)
traindata['std'] = traindata[[i for i in range(10000)]].std(axis=1)
traindata['var'] = traindata[[i for i in range(10000)]].var(axis=1)
traindata['skew'] = traindata[[i for i in range(10000)]].skew(axis=1)
traindata['mean'] = traindata[[i for i in range(10000)]].mean(axis=1)
traindata['median'] = traindata[[i for i in range(10000)]].median(axis=1)
testdata['max'] = testdata[[i for i in range(10000)]].max(axis=1)
testdata['min'] = testdata[[i for i in range(10000)]].min(axis=1)
testdata['std'] = testdata[[i for i in range(10000)]].std(axis=1)
testdata['var'] = testdata[[i for i in range(10000)]].var(axis=1)
testdata['skew'] = testdata[[i for i in range(10000)]].skew(axis=1)
testdata['mean'] = testdata[[i for i in range(10000)]].mean(axis=1)
testdata['median'] = testdata[[i for i in range(10000)]].median(axis=1)
# 入模特征选择
cols = [i for i in range(10000)] + ['gender','max','min','std','var','skew','mean','median']
- 模型训练
def catboost_model(train_x, train_y, test_x, seed = 2023):
folds = 5
kf = KFold(n_splits=folds, shuffle=True, random_state=seed)
oof = np.zeros(train_x.shape[0])
# 样本数
test_predict = np.zeros(test_x.shape[0])
cv_scores = []
for i, (train_index, valid_index) in enumerate(kf.split(train_x, train_y)):
print('************************************ {} ************************************'.format(str(i+1)))
trn_x, trn_y, val_x, val_y = train_x.iloc[train_index], train_y[train_index], train_x.iloc[valid_index], train_y[valid_index]
params = {'learning_rate': 0.1,
'depth': 5,
'bootstrap_type':'Bernoulli',
'random_seed':2023,
'od_type': 'Iter',
'od_wait': 100,
'random_seed': 11,
'allow_writing_files': False,
'task_type':"GPU",
'devices':'0:1'}
##iterations是迭代轮数，可以根据自己的算力与精力，选择合适的轮数，示例给出的是500
model = CatBoostRegressor(iterations=3000, **params)
model.fit(trn_x, trn_y, eval_set=(val_x, val_y),
metric_period=500,
use_best_model=True,
cat_features=[],
verbose=1)
val_pred = model.predict(val_x) # 进行预测
test_pred = model.predict(test_x)
oof[valid_index] = val_pred # 验证集的预测存储
test_predict += test_pred / kf.n_splits # 交叉验证取均值
score = mean_absolute_error(val_y, val_pred) # MAE
cv_scores.append(score)
print(cv_scores)
# 获取特征重要性打分，便于评估特征
if i == 0:
fea_ = model.feature_importances_
fea_name = model.feature_names_
fea_score = pd.DataFrame({'fea_name':fea_name, 'score':fea_})
fea_score = fea_score.sort_values('score', ascending=False)
fea_score.to_csv('feature_importances.csv', index=False)
return oof, test_predict
cat_oof, cat_test = catboost_model(traindata[cols], traindata['age'], testdata[cols])
- 打印结果
testdata['age'] = cat_test
testdata['age'] = testdata['age'].astype(float)
testdata['age'] = testdata['age'].apply(lambda x: x if x>0 else 0.0)
testdata['age'] = testdata['age'].apply(lambda x: '%.2f' % x)
testdata[['sample_id','age']].to_csv('submit_2.txt',index=False)
扫一扫
203
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
