R语言中的geneHapR包做单倍型/单倍型网络分析（1）

这个R包对应的论文
https://bmcbioinformatics.biomedcentral.com/articles/10.1186/s12859-023-05318-9
geneHapR: an R package for gene haplotypic statistics and visualization
对应的github链接
https://github.com/ZhangRenL/geneHapR
https://gitee.com/zhangrenl/genehapr 这个链接的wiki部分有比较详细的教程
正好今天看到一篇论文
https://www.nature.com/articles/s41467-025-60436-7
Allelic variations in GA20ox3 regulate fruit length and seed germination timing for high-altitude adaptation in Arabidopsis thaliana
论文中的Figure5是单倍型的分析
论文中提供了这个部分的比对好的fasta文件
论文中的方法部分写到
The −412 to −434 bp segment of the GA20ox3 promoter region from 576 accessions was analyzed to construct a statistical parsimony network. This haplotype network was generated using Network 10.2 software
这里有一个疑问是为啥用专门的那三个碱基去做单倍型，通常某个基因的区间内或者上下游会有很多变异位点，如何选择变异位点做单倍型分析暂时还没有想明白。
论文中的FigS11，这三个碱基正好是位于一个转录因子结合基序的上面

从数据中把三个碱基的位点挑出来对应的位置是185-188，用geneHapR来分析单倍型
geneHapR 中有函数可以直接读vcf文件，也可以直接读fasta文件

library(geneHapR)

library(tidyverse)

seqs<-import_seqs(“genehapR/abc.txt”) 单倍型分析 seqs2hap(seqs) -> hapResults
hapSummary<-hap_summary(hapResults)
画图展示不同类别的单倍型
plotHapTable(hapSummary)

这个和论文中的结果符合，只不过有两个样本有gap被过滤掉了
单倍型网络分析
hapnet <- get_hapNet(hapSummary)
plotHapNet(hapnet)

给单倍型网络一个分组
序列里的id对应的是 分组数据框的行名，我这里随机构造了
group.info<-data.frame(x1=names(seqs), x2=paste0(“group”,sample(1:3,576,replace = TRUE))) %>%
column_to_rownames(“x1”)

hapnet <- get_hapNet(hapSummary,
AccINFO = group.info,
groupName = “x2”)

plotHapNet(hapnet,legend = TRUE,show_color_legend = TRUE,show_size_legend = FALSE)

这里比较有意思的是图里是需要借助鼠标点击一下作图区域，鼠标点击哪里图例就出现在哪里，这个功能不知道是怎么实现的

不同的单倍型对应的表型，对应的是论文中的fig5fghi
example.pheno<-data.frame(x1=names(seqs), x2=sample(1:1000,576,replace = TRUE)) %>%
column_to_rownames(“x1”)

example.pheno %>% dim()

hapVsPheno(hap = hapResult,
pheno = example.pheno,
hapPrefix = “H”,
filename.prefix = “A”)
这个命令通常可以直接出图，我这里遇到了报错
Error in hapVsPheno(hap = hapResult, pheno = example.pheno, hapPrefix = “H”, :
After removed NAs, observations for ‘x2’ is not enough.
此外: Warning message:
In hapVsPheno(hap = hapResult, pheno = example.pheno, hapPrefix = “H”, :
phenoName is null, will use the first pheno
暂时看不出来是什么问题，有时间再研究吧