不同单倍型和表型数据的显著性分析

前几天写了一些单倍型分析的博客：如何计算群体中的单倍型频率，单倍型的显著性分析，以及介绍了为何要做单倍型分析：GWAS分析完，要做单倍型图，还要做单倍型的显著性分析？，以及每个个体的单倍型：单倍型分析：个体所对应的单倍型是？

数据还是用之前Haploview的数据，单倍型分析全套教程参考：从入门到出家：单倍型Haploview分析（万字详解），配套数据：

图片

1，先看一下文件单倍型划分

图片

图片

发现，10号染色体的，这个区间，有几个位点位于一个block里面。

2，提取block的文件

图片

看一下map数据：

图片

没问题，就是这几个。

3，将数据变为vcf

图片

4，使用geneHapR包处理vcf

## 导入基因型vcf数据
library(geneHapR)

vcf = import_vcf("df2.vcf")

# 单倍型分型
hapResult <- vcf2hap(vcf)

write.csv(hapResult,"df2-1-hapresult.csv")

结果文件：

图片

图片

上面结果中，共有11个单倍型，每个个体都会给出具体的单倍型分型。

5，表型数据整理

表型数据准备，两列数据，第一列是ID，第二列是分析的性状，整理成txt文件，tab分割。

图片

6，读取表型数据

使用函数import_AccINFO，读取txt文件，自动将第一列ID变为行名。

图片

7，单倍型和表型数据显著性分析

图片

结果文件：

图片

结果文件中，可以看出，共有两个单倍型：H001和H002，其中H001有10个个体，H002有8个个体，观测性状为y1，两个单倍型之间没有显著性。

8，如果单倍型之间显著，是什么样子的

我们用示例数据展示一下：

data("geneHapR_test")
# plot the figs directly
hapVsPheno(hap = hapResult,
           pheno = pheno,
           phenoName = "GrainWeight.2021",
           minAcc = 3)

图片