UKBiobank Olink蛋白数据预处理-python

最新推荐文章于 2024-09-13 19:41:31 发布

MedicineAI

最新推荐文章于 2024-09-13 19:41:31 发布

阅读量468

点赞数 6

文章标签： python 数据挖掘开发语言

本文链接：https://blog.csdn.net/MedicineAI/article/details/138719552

版权

项目来源：GitHub - jasonHKU0907/FutureHealthProteomicPrediction

一、读取原始数据

读取原始olink数据，并只使用基线时间的数据，随访数据不考虑

mydf = pd.read_csv(dpath + 'olink_data.csv')
mydf = mydf.loc[mydf.ins_index == 0]   # 基线时间数据

提取患者id和蛋白id的列表

eid_lst = list(np.sort(list(set(mydf.eid))))   # 提取患者id的list
protein_lst = list(np.sort(list(set(mydf.protein_id))))    # 提取蛋白id的list

由于Olink蛋白数据的保存为下图格式，每一个患者都有多个蛋白数据，即同一个eid对应了多个protein_id行，因此，为了转化为一个患者id一个数据行的格式，需要对每个eid依次读取。

由于患者数量5w+，将他们以步长5000分为11个epoch读取。依次读取每个患者的所有行，融合此数据与完整protein_id。每个患者一列依次拼接组合成protein_id行 * eid+1列的矩阵。转置后成为我们需要的格式。

myout_df_final = pd.DataFrame()
step = 5000

for i in range(11):
    protein_tmp_df = pd.DataFrame({'protein_id': protein_lst})
    myout_df = pd.DataFrame({'protein_id': protein_lst})
    j = 0
    for eid_item in eid_lst[i*step:(i+1)*step]:   # 依次选取5000个患者id中的每一个
        tmpdf0 = mydf.loc[mydf.eid == eid_item]   # 此患者数据所有数据行
        tmpdf1 = tmpdf0[['protein_id', 'result']]    # 提取'protein_id'和'result'列
        merged_df = pd.merge(protein_tmp_df, tmpdf1, how='left', on=['protein_id'])    # 按protein_id列，将此患者result数据与完整protein_id融合
        myout_df = pd.concat([myout_df, merged_df.iloc[:, 1:2]], axis=1)   # 按列拼接每一个患者的'result'列
        myout_df.rename(columns={'result': eid_item}, inplace=True)   # 用患者的id命名他的数据列
        print(j)
        j+=1
    myout_df = myout_df.transpose()   # 转置后每一行为一个患者蛋白数据行
    myout_df.columns = myout_df.iloc[0].astype(int).tolist()
    myout_df.drop('protein_id', axis=0, inplace=True)
    myout_df['eid'] = myout_df.index.astype(int).tolist()
    myout_df_final = pd.concat([myout_df_final, myout_df], axis = 0)
    myout_df_final.to_csv(dpath + 'Protemics_ins0.csv', index = False)
    print(i)

二、缺失值处理

1、蛋白的缺失值

pro_df = myout_df_final
my_eid_lst = pro_df.eid.tolist()   # 患者id list
nb_eids = len(my_eid_lst)    # 患者个数
my_pros_lst = pro_df.columns.tolist()[:-1]   # 蛋白id list

# 每一个蛋白中的NA个数列表
na_prop_lst = []
for my_pros in my_pros_lst:  # 对每一个蛋白列，找出其中的na个数，计算na比列
    tmpdf = pro_df[my_pros]   
    nb_na = len(np.where(tmpdf.isna() == True)[0])
    na_prop_lst.append(np.round(nb_na/nb_eids, 3))   # 比例=na个数/患者个数，保留3位小数

na_pros_df = pd.DataFrame({'Pros': my_pros_lst, 'NA_prop':na_prop_lst})   # 构造一个2列的表，一列是蛋白id，一列是对应NA比例
na_pros_df.sort_values(ascending=False, by = ['NA_prop'], inplace = True)    # 按NA比例降序重排
pros_beyond_na30_lst = na_pros_df.loc[na_pros_df.NA_prop>0.3].Pros.tolist()   # 找出NA比例大于30%的蛋白id
pro_df.drop(pros_beyond_na30_lst, axis = 1, inplace=True)   # 删除这些蛋白id列

2、患者的缺失值

nb_pros = pro_df.shape[1] -1   # 蛋白个数

# 每一个患者中的NA个数列表
na_prop_pros_lst = []
for my_eid in my_eid_lst:   # 对每一个患者行，找出其中的NA个数，计算NA比列
    tmpdf = pro_df.loc[pro_df.eid == my_eid]
    nb_na = len(np.where(tmpdf.isna() == True)[0])
    na_prop_pros_lst.append(np.round(nb_na/nb_pros, 3))   # 比例=NA个数/蛋白个数，保留3位小数

na_eid_df = pd.DataFrame({'eid': my_eid_lst, 'NA_prop':na_prop_pros_lst})  # 构造一个2列的表，一列是患者id，一列是对应NA比例
na_eid_df.sort_values(ascending=False, by = ['NA_prop'], inplace = True)   # 按NA比例降序重排
eid_beyond_na30_lst = na_eid_df.loc[na_eid_df.NA_prop>0.3].eid.tolist()    # 找出NA比例大于30%的患者id
pro_df = pro_df[~pro_df.eid.isin(eid_beyond_na30_lst)]   # 不属于这些患者id的行组成新的总表

3、插补

把患者所在UKB评估中心的地理位置合并到英国的10个地区作为划分标准，对这10个地区插补时做留一交叉验证。

pro_f_lst = pro_df.columns.tolist()[1:]   # protein id
# 每个患者及其所对应的region代码
region_code_df = pd.read_csv(dpath + 'Eid_info_data.csv', usecols = ['eid', 'Region_code'])
pro_df = pd.merge(region_code_df, pro_df, how = 'right', on = 'eid')
region_code_lst = [0, 1, 2, 3, 4, 5, 6, 7, 8, 9]

首先对蛋白数据进行标准化，以降低离群值的影响。

def get_standardization(mydf):
    tmp_df = mydf.copy()
    for col in tmp_df.columns:
        ubd = tmp_df[col].quantile(0.995)
        lbd = tmp_df[col].quantile(0.005)
        tmp_df[col].iloc[tmp_df[col]>ubd] = ubd
        tmp_df[col].iloc[tmp_df[col]<lbd] = lbd
        tmp_df[col] = (tmp_df[col] - lbd) / (ubd - lbd)
    return tmp_df

mydf[pro_f_lst] = get_standardization(mydf[pro_f_lst])

其次将所有患者按照性别分为两组，做插补时加上人工选出的7个人口统计学数据，分别为Age，Ethnicity_White, Ethnicity_Asian, Ethnicity_Black, Ethnicity_Others, TDI (Townsend deprivation index at recruiment), Education。

# Male组
mydf_male = mydf.loc[mydf.DM_GENDER == 1]
mydf_out_male = pd.DataFrame()
time0 = time.time()
# 人工选出的7个因素
demo_f = ['DM_AGE', 'DM_ETH-W', 'DM_ETH-A', 'DM_ETH-B', 'DM_ETH-O', 'DM_TDI', 'DM_EDUC']

KNN插补操作（以Male组为例）：首先将所有男性按照10个地区的编码进行划分，其中一个区域作为测试集，其他作为训练集，通过调用sklearn.impute.KNNImputer，并计算与每一个蛋白pearson相关系数最大的50个蛋白，利用这50个相关蛋白与对缺失值进行插补。

for region_id in region_code_lst:    # 0~9
    tmpdf_train = mydf_male.loc[mydf_male.Region_code != region_id]
    tmpdf_test = mydf_male.loc[mydf_male.Region_code == region_id]   # 留一交叉验证
    tmpdf_test.reset_index(inplace = True)
    corr_df = tmpdf_train[pro_f_lst].corr().abs()    # 计算每列protein之间的pearson相关系数的绝对值
    i = 0
    for pro_f in pro_f_lst:   # 对于每一个protein，在tmpdf_test中更新插补完的protein列
        # 与pro_f相关性最强的前50个protein ID
        f2select = get_correlated_pros(pro_f = pro_f, nb_f2select = 50, corr_df = corr_df)
        tmpdf_train_sub = tmpdf_train[f2select + demo_f]
        tmpdf_test_sub = tmpdf_test[f2select + demo_f]
        # 训练集中的50个权重相同的邻居的平均值进行插补
        knn_imputer = KNNImputer(n_neighbors=50, weights="uniform")
        knn_imputer.fit(tmpdf_train_sub)
        tmpdf_test_sub_imputed = pd.DataFrame(knn_imputer.transform(tmpdf_test_sub))
        tmpdf_test_sub_imputed.columns = tmpdf_test_sub.columns.tolist() 
        tmpdf_test[pro_f] = tmpdf_test_sub_imputed[pro_f]   # 把插补好了的某一个蛋白列填回原test数据
        i+=1
        if i%500 == 0:    # 每500个protein print一次
            print((region_id, i, time.time() - time0))
    mydf_out_male = pd.concat([mydf_out_male, tmpdf_test], axis=0)