具有遗传性疾病和性状的遗传位点分析
在正式参加研究僧数学建模比赛之前，想先做做题目练练手，因为是一名统计学的研究僧，所以选择了一道统计题，欢迎大家交流指正。实现的代码基本用的python，最后一题用了R语言~
题目数据与完整代码都在：https://github.com/hyh1608/MCM_2016
研究背景
人体的每条染色体携带一个 DNA 分子，人的遗传密码由人体中的 DNA 携带。 DNA 是由分别带有 A、T、C、G 四种碱基的脱氧核苷酸链接组成的双螺旋长链分子。在这条双螺旋的长链中，共有约 30 亿个碱基对，而基因则是 DNA 长链中有遗传效应的一些片段。在组成 DNA 的数量浩瀚的碱基对（或对应的脱氧核苷酸）中，有一些特定位置的单个核苷酸经常发生变异引起 DNA 的多态性，我们称之为位点。染色体、基因和位点的结构关系见下图
本题目针对某种遗传疾病(简称疾病 A)提供了 1000 个样本的信息，这些信息包括样本的患病与否信息（见附件 phenotype.txt）、样本的某条可能致病染色体 片段上 9445 个位点编码信息（见附件 genotype.dat），以及包含这些位点的 300 个基因的信息（见附件文件夹 gene_info），其中每个基因的信息由它所包含的 位点名称组成。
另外，人体的许多遗传疾病和性状是有关联的，比如高血压，心脏病、脂肪 肝和酒精依赖等症状具有关联性，科研人员往往把相关的性状或疾病放在一起研 究，这样能提高发现致病位点或基因的能力，因此本题还提供了这 1000 个样本 的 10 种相关性状表现与否的信息（见附件 multi_phenos.txt）。
根据题目中给出的以上数据解决以下问题：
问题一：为了方便进行数据分析，请使用适当的方法，将附件 genotype.dat 中每个位点的碱基(A,T,C,G) 编码方式转化成数值编码方式。
问题二：根据重新编码的附件 genotype.dat 和原始附件 phenotype.txt，设计或采用一个方法，找出疾病 A 最有可能的一个或几个致病位点，并给出相关的理论依据。
问题三：由于可以把基因理解为若干个位点组成的集合，遗传疾病与基因的关联性可以由基因中包含的位点的全集或其子集表现出来，根据重新编码的附件genotype.dat、原始附件 phenotype.txt，以及原始附件文件夹 gene_info，找出与疾病 A 最有可能相关的一个或几个基因，并给出相关的理论依据。
问题四：在实际的研究中，科研人员往往把相关的性状或疾病看成一个整体, 然后来探寻与它们相关的位点或基因，根据重新编码的附件 genotype.dat 和原始附件 multi_phenos.txt 找出与 multi_phenos.txt 中 10 个性状有关联的位点。
模型假设
- 假设样本数据真实，来源可靠，信息准确；
- 假设基因内的位点之间相互独立；
- 假设基因之间相互独立，不存在相互促进或者相互抑制的情况；
- 假设基因与遗传疾病之间只存在简单线性关系。
问题一：将碱基编码转换为数值编码
在后续三题的求解中，只涉及每条位点中三种编码的比例， 也就是说每条位点只与自身的三种编码方式的比例有关，与位点内的碱基类型无 关，所以六种类型统一采用 A/B 表示，而三种编码 AA、AB、BB 利用 0、1、2 重 新编码，这样有利于后文的统计分析。（原始数据中有两个位点（rs4252972 和 rs11573221）的原始 编码出现了字母 I 和 D，不属于四种碱基，我们认为这是数据出现了偏差，但是由于我们都将编码换成了 A 和 B，因此并不影响我们的重新编码过程）
import numpy as np
import pandas as pd
# 读取 数据
os.chdir('D:\\study\\研究生\\数学建模\\建模比赛')
data = pd.read_csv('./建模培训/第一轮/研究生A题-2016B遗传位点分析/genotype.csv')
# 检验数据是否有缺失值
print(data[data.isna().values==True]) # 0x9445没有缺失值
# 对每一位点不同碱基对的数量进行统计
def apply_count(arr):
dict_count = {}
for i in range(len(arr)):
if arr[i][0] not in dict_count:
dict_count[arr[i][0]] = 1
else:
dict_count[arr[i][0]] += 1
return dict_count
position = [] # 用来存放所有位点不同碱基对的数
for i in range(len(data.columns)):
column = data.columns[i]
# trans to numpy
dict_count = data[[column]].values
position.append(apply_count(dict_count))
# position[0]:{'TT':664,'TC':293,'CC':43}
# 编码时默认将AB型的碱基对默认为1，0为主纯碱基对，2为次纯碱基对，所以先删除AB型碱基对
for i in range(len(position)):
for key in position[i]:
if key[0] != key[1]:
del position[i][key]
break
def trans_code(data):
trans_arr = []
for i in range(len(data.columns)):
column = data.columns[i]
data_c = data[column].values.tolist()
for j in range(len(data_c)):
# 将每一个位点AB型的碱基对直接编码为1
if data_c[j][0] != data_c[j][1]:
data_c[j] = 1
else:
k = data_c[j]
# 如果AA和BB型的碱基对数量相等，则默认TT为2，CC为0
if len(set(position[i].values())) == 1:
if k == 'TT':
data_c[j] = 2
if k == 'CC':
data_c[j] = 0
else:
# 将主纯碱基对编码为0，次纯碱基对编码为2
if position[i][k] == max(position[i].values()):
data_c[j] = 0
else:
data_c[j] = 2
trans_arr.append(data_c)
return trans_arr
# 将每一列编码
trans_encode = trans_code(data)
trans_arr = np.array(trans_encode).T
# 把编码后的数组储存成数据框的形式
df = pd.DataFrame(trans_arr,columns=data.columns)
问题二：找出疾病 A 最有可能的致病位点
问题分析：首先采用最小等位基因频率控制法筛选具有统计意义的位点**最小等位基因频率（Minor Allele Frequency，简称 MAF）是指在给定人群中 不常见的等位基因发生的频率。**较小的 MAF 会降低统计效能，通常情况下，我们取 MAF 的阈值为 0.05，若某位点的 MAF<0.05，则认为该位点异常，将其从数据中删除。随后我们引进卡方统计量来描述某位点编码在健康和患病两组中的分布差异性，统计量的值越大说明差异性越大，进而说明该位点与患病的关联性越强。
# 导入疾病信息，0表示未患病，1表示患病
label = pd.read_table('./建模培训/第一轮/研究生A题-2016B遗传位点分析/phenotype.txt',header=None)
# 数据清洗，采用最小等位基因频率进行筛选
List1 = [] # list1每一项储存三个参数，第一个参数为位点的索引值，第二个参数为频率大的等位基因出现的概率
# 第三个参数为频率小的等位基因出现的概率
List2 = []
for i in range(len(position)):
n_0 = 0
n_1 = 0
n_2 = 0
temp = []
for j in trans_arr[:,i]:
if j == 0:
n_0 += 1
if j == 1:
n_1 += 1
if j == 2:
n_2 += 1
pro1 =( n_1 + n_0 * 2 ) / 2000
pro2 = (n_1 + n_2 * 2 ) / 2000
temp.append(i)
temp.append(pro1)
temp.append(pro2)
List1.append(temp)
# List2用来储存基因概率小于0.05的位点索引，
for i in range(len(List1)):
if List1[i][1] < 0.05 or List1[i][2] < 0.05:
List2.append(i)
print(len(List2)) # 28一共删除个位点
# 将List2中的位点给删除
for i in df.columns[List2]:
df1 = df.drop(i,1)
df = df1
# 将label添加至最后一列
df1['label'] = label
df1_0 = df1[df1.label==0] # 获取未患病的样本
df1_1 = df1[df1.label==1] # 获取患病的样本
# 获取每一位点中0，1，2编码的数量
def apply_count1(arr):
dict_count = {0:0,1:0,2:0}
for i in range(len(arr)):
if arr[i][0] not in dict_count:
dict_count[arr[i][0]] = 1
else:
dict_count[arr[i][0]] += 1
return dict_count
# 健康的人的人每种碱基对的数量
health = []
for i in range(len(df1_0.columns)):
column = df1.columns[i]
# trans to numpy
dict_count = df1_0[[column]].values
health.append(apply_count1(dict_count))
# 不健康的人每种碱基对的数量
unhealth = []
for i in range(len(df1_1.columns)):
column = df1.columns[i]
dict_count = df1_1[[column]].values
unhealth.append(apply_count1(dict_count))
# health_statics 存放列联表,每个位点都有一张列联表
health_statis = []
for i in range(len(df1.columns)):
temp = []
temp.append([health[i][0],health[i][1],health[i][2]])
temp.append([unhealth[i][0],unhealth[i][1],unhealth[i][2]])
health_statis.append(temp)
# 存放卡方检验的p值
p_value = []
for i in range(len(health_statis) - 1):
kf_data = np.array(health_statis[i])
kf = chi2_contingency(kf_data)
# print(i,kf[1])
p_value.append(kf[1])
p_value_arr = np.array(p_value)
# List3用来存放p值小于0.01的位点索引与p值,一共有73个
List3 = []
for i in range(len(p_value_arr)):
temp = []
if p_value_arr[i] < 0.01:
temp.append(i)
temp.append(p_value_arr[i])
List3.append(temp)
第二题的曼哈顿图
问题三：找出与疾病 A 最有可能相关的基因
问题三与问题二相似，都是要找出两个变量之间的相关性，在寻找显著位点的过程中已经使用一次假设检验，如果在此基础上再利用假设检验去寻找致病基因与位点的从属关系，会引入更多的不确定性，犯错误的概率会增加。因此本题不使用问题二的结果，而是把位点集合作为自变量来看待，再加上因变量是否患病为二分类变量，所以考虑本题使用逻辑回归进行求解。
# 统计的模块
import statsmodels.api as sm
# 遍历读取当前文件夹下的所有文件
def file_name(file_dir):
for root,dirs,files in os.walk(file_dir):
gene = []
for file in files:
# print(file)
with open(file_dir+file) as f:
# 取出对应的基因
# print(file_dir+file)
num = (file.split('.')[0]).split('_')[1]
# print(num)
lines = f.readlines()
SNP = []
# 存放索引值与位点
SNP.append(num)
for line in lines:
line = line.strip()
SNP.append(line)
gene.append(SNP)
return gene
gene = file_name('./建模培训/第一轮/研究生A题-2016B遗传位点分析/gene_info/')
# df_gene存放着300张表，每张表对应一个基因
df_gene = []
for i in range(len(gene)):
temp_gene = []
for j in gene[i][1:]:
if j in df1.columns:
temp_gene.append(df1[j])
# temp_gene.append(df1['label'])
df2 = pd.concat(temp_gene,axis=1)
df2 = pd.concat([df2,df1['label']],axis=1)
df_gene.append(df2)
# 采用逻辑回归
def logit(df_gene):
train_cols = df_gene.columns[:-1]
logit = sm.Logit(df_gene['label'],df_gene[train_cols])
result = logit.fit()
# p值为方程的显著性
p = result.llr_pvalue
return p
# 300个模型的p值
p_list = []
for i in range(len(df_gene)):
temp = []
p = logit(df_gene[i])
temp.append(i)
temp.append(p)
p_list.append(temp)
# 从小到大排序并返回其索引值,取出前5
print(np.argsort(np.array(p_list)[:,1])[:5]) #array([ 4, 162, 151, 164, 25], dtype=int64)
第三题的曼哈顿图
问题四：与 10 个性状有关联的位点
典型相关分析（Cononical Correlation Analysis）是研究两组变量之间相关关系的多元分析方法。它借用主成分分析降维的思想，分别对两组变量提取主成分， 且使从两组变量提取的主成分之间的相关程度达到最大，而从同一组内部提取的各主成分之间互不相关，用从两组分别提取的主成分的相关性来描述两组变量整体的线性相关关系。
# 用R语言进行典型相关分析
test1 = read.csv('df_select.csv')
test1 = test1[,-1]
ca = cancor(test1[,1:28],test1[,29:38])
corcoef.test<-function(r, n, p, q, alpha=0.1){
m<-length(r); Q<-rep(0, m); lambda <- 1
for (k in m:1){
lambda<-lambda*(1-r[k]^2);
Q[k]<- -log(lambda)
}
s<-0; i<-m
for (k in 1:m){
Q[k]<- (n-k+1-1/2*(p+q+3)+s)*Q[k]
chi<-1-pchisq(Q[k], (p-k+1)*(q-k+1))
# if (chi>alpha){
# i<-k-1; break
# }
print(chi)
s<-s+1/r[k]^2
}
i
}
corcoef.test(r=ca$cor,n=1000,p=28,q=10)
写在最后
- 代码可能看上去比较繁琐，本人小白一枚，后续如果有更好的方法，会来及时更新
- 第四小题用典型相关分析求出的特征向量a与b特别小，希望大家能有更好的方法来解决第四题~
扫一扫
621
被折叠的 条评论
为什么被折叠?
到【灌水乐园】发言
