test4.py 2.3 KB

123456789101112131415161718192021222324252627282930313233343536373839404142434445464748495051525354555657
  1. import pandas as pd
  2. import numpy as np
  3. def calculate_ed(y, genotype):
  4. """
  5. 计算ED值
  6. :param y: 表现型数据,一维numpy数组
  7. :param genotype: SNP基因型数据,Pandas DataFrame,行数为个体数量,列数为SNP位点数
  8. :return: ED值
  9. """
  10. # 检查输入数据是否合法
  11. if not isinstance(y, np.ndarray) or not isinstance(genotype, pd.DataFrame):
  12. raise TypeError("输入数据格式错误")
  13. # 计算每个个体的基因型得分
  14. score = genotype.sum(axis=1)
  15. # 计算每个个体的方差
  16. var_score = np.var(score)
  17. # 计算每个个体与平均表现型的差异的平方和
  18. sum_delta_y_sq = np.sum((y - np.mean(y))**2)
  19. # 计算ED值
  20. ed = sum_delta_y_sq / var_score
  21. return ed
  22. def main():
  23. # 读取MAF文件
  24. # df = pd.read_csv('example.maf', sep='\t', comment='#')
  25. # 读取MAF文件
  26. df = pd.read_csv('69dbfe7c-3efc-4126-b36d-2b59aca7a16d.wxs.aliquot_ensemble_masked.maf', sep='\t', comment='#')
  27. # 遍历每一列并将非数字值转换为 NaN
  28. for col in df.columns:
  29. # 将每列中的数据转换为数字类型
  30. df[col] = pd.to_numeric(df[col], errors='coerce')
  31. # 会使用 apply() 方法来遍历每列数据,并将大于等于1000的值设置为 NaN。如果你只想保留小于10000的数字,可以将 x < 1000 改为 x >= 10000。
  32. df[col] = df[col].apply(lambda x: x if x is np.nan or x < 10000 else np.nan)
  33. # 删除包含 NaN 的列
  34. df = df.dropna(axis=1)
  35. # 提取表现型数据和基因型数据
  36. # y = df.iloc[:, 0].values
  37. # 这里提取关键数据作为对比用户上传的SNP数据来进行对比,如果ed值越高,则意味着改用户换胃癌的风险就越高
  38. lst = [75, 420, 73, 52, 72, 62, 82, 36, 102, 56, 90, 224, 160, 66, 99, 87, 24, 93, 85, 73, 70, 38, 52, 159, 102, 37,
  39. 37, 33, 49, 65, 20, 67, 74, 84, 64, 85, 16, 12, 42, 60, 38, 49, 86, 72, 97, 176, 63, 76, 61, 54, 153, 39, 53,
  40. 70, 118, 102, 25, 49, 17, 114, 48, 52, 110, 24, 90, 61, 48, 71, 30, 48, 41, 32, 56, 59, 52, 84, 52, 53, 47,
  41. 117, 152, 48, 68, 93, 116]
  42. y = np.array(lst)
  43. genotype = df.iloc[:, 1:].to_dict()
  44. # 计算 ED 值,并保留3哥小数位
  45. ed = round(calculate_ed(y, pd.DataFrame(genotype)), 3)
  46. # print('ed', ed)
  47. return ed
  48. main()