单倍型分析
一、什么是单倍型？
在单倍型分析前，首先需要明白什么是单倍型、什么是基因单倍型？
一般来讲，所谓的单倍型是指同一染色体上不同变异位点的各种线性组合形式。那么基因单倍型自然就是指同一基因上（启动子、外显子、内含子、终止子）的不同变异位点的显性组合形式。
说明：基因组或染色体水平的单倍型分析不在本文范畴之内，请您参考“https://www.jianshu.com/p/4de7762aa81e”。
二、单倍型分析需要什么数据？数据从哪里来？
那么，进行基因单倍型分析需要那些数据呢？
1）基因型数据（必须）
既然单倍型分析需要变异位点信息，自然就需要基因型数据。基因型数据可以有很多来源，首先，可以零成本从各大公开数据库下载（推荐）；其次，可以用一代测序的方法获取基因型数据（适合于基因不太长、样本数量较少的情况）；如果你们课题组经费非常充足，可以自行做制作基因芯片、高通量测序等等。
数据格式：
- 一代测序数据：需要拼接成序列并保存成fasta格式
- 二代测序数据：VCF，Hapmap，plink（map&ped）
- 自定义格式：表格，每行为一个变异位点，每列为一个样品，前五列固定为Chr，POS，REF，Alt和INFO
2）基因注释信息（可选）
如果需要查看单倍型中的变异位点在基因上面的分布的话，则会需要用到基因注释信息。如果没有这方面需求，则可以忽略。一般来说，各个已测序的物种都有对应的注释信息，这个做基因功能研究的肯定知道从哪里搞到适合你们所研究的物种的注释信息的。比如已测序的作物可以从phytozome数据库中下载，或者物种对应的数据库下载玉米（maizeGDB）。当然，这种注释信息不一定百分百正确。如果非常不行的你研究的基因恰巧注释是错误的，那么就需要你自行创建一个BED格式的注释文件了。
数据格式
- GFF/GFF3：一般从数据库下载的注释文件属于这种格式，可以手动将对应基因的注释信息提取出来，也可以直接使用源文件
- BED6：一般自定义为BED6格式；第一列为染色体名称，第二列为起始位置，第三列为终止位置，第四列为片段名称（基因ID+空格+CDS/UTR），第五列留空，第六列为方向（+/-）
3）样本信息（表型数据、地理坐标、样本分类等等，可选）
在分析完单倍型之后需要评估优势单倍型（产量高、高耐逆等）、变异位点的效应分析，则需要用到表型数据。单倍型的地理分布需要用到地理坐标。如果提供样本分类数据的话，单倍型网络分析能提供更多的线索。
数据格式
- 文件中：每行为一个样本，每一列作为一种数据，
- R中： 数据框（data.frame），行名为样本名，列名为数据信息（type，subgroup，location，等）
三、单倍型需要什么工具？
什么，都2302年了，您还在用Excel分析单倍型呢？？！
在这里向推荐大家使用geneHapR。这款软件完美地支持从fasta序列以及vcf、p.link、HapMap等格式的二代测序结果进行单倍型鉴定。
软件安装：由于这款软件是基于R语言开发的，所以需要首先安装R和RStudio两款软件（注意一定要先安装R再安装RStudio）。
最新版R下载链接：https://cloud.r-project.org/bin/windows/base/。
最新版RStudio下载链接：https://posit.co/download/rstudio-desktop/。
都安装完成后，打开RStudio，输入如下命令。
# 先安装一些依赖的BiocManager的包，否则后续安装geneHapR容易出错
install.packages("BiocManager")
BiocManager::install(c("Biostrings", "GenomicRanges", "muscle", "IRanges", "rtracklayer", "trackViewer"))
# 安装geneHapR
install.packages("geneHapR")
如果安装过程遇到错误了。请不要灰心，不要放弃，命运总是坎坷的，相信阳光总在风雨后！请先安装缺失的软件包，等完成后再尝试运行install.packages("geneHapR")安装geneHapR。
如果看到程序包‘geneHapR’打开成功，MD5和检查也通过，恭喜你可以使用geneHapR进行单倍型分析了！！！
四、工具装好了，该怎么操作呢？
4.0 测试数据准备
将下列4个文件下载下来并存放到工作目录中，后续分析会用到
基因型数据：Genotype
注释信息：Annotation
表型数据：Phenotype
其他相关信息：AccINFO
4.1 RStudio代码
如果有一定的R语言基础，推荐您使用这种方式。
1） 从数据导入到单倍型鉴定
# 首先把软件加载进来
library(geneHapR)
# 设置工作目录（windows的同学注意"\"和"/"的问题）
setwd("D:/Haplotype")
# 导入各种数据
gff <- import_gff("gff/OsGHD7.gff3") # 导入GFF格式的注释数据
gff <- import_bed("12859_2023_5318_MOESM3_ESM.bed6") # 导入BED格式的注释数据
pheno <- import_AccINFO("12859_2023_5318_MOESM4_ESM.tsv") # 导入表型数据
AccINFO <- import_AccINFO("12859_2023_5318_MOESM5_ESM.csv",
sep = ",", # 分隔符号，默认为制表符"\t"
na.strings = "NA") # 导入其他样本信息
# 导入VCF格式的基因型数据及变异信息的筛选
vcf <- import_vcf("OsGHD7.vcf.gz")
vcf <- filter_vcf(vcf, gff,
mode = "both", # both, POS, Type的其中一个
start = start, end = end, Chr = Chr, # 保留哪些位置的变异信息
type = "CDS", cusTyp = c("CDS")) # 保留哪些类型的变异信息
# 导入表格形式的基因型数据
tbl <- read.csv("12859_2023_5318_MOESM2_ESM.geno")
# 导入其他格式的基因型数据请查看帮助
# 这里就不再一一展示了
help("geneHapR")
library(help = 'geneHapR')
browseVignettes('geneHapR')
# 一些基本参数的设置
geneID <- "OsGHD7" # 基因ID
Chr <- "Chr7" # 基因所处的染色体名称
start <- 9152403 # 基因的起始位置（染色体坐标）
end <- 9155185 # 基因的终止位置（染色体坐标）
hapPrefix <- "H" # 单倍型名称的前缀
# 常用格式的单倍型鉴定，具体参数请查看帮助文档
# VCF: vcf2hap()
# P.link(ped&map): plink.pedmap2hap()
# Fasta: seqs2hap()
# HapMap: hmp2hap()
# Table: table2hap()
# 从VCF开始单倍型鉴定
hapResult <- vcf2hap(vcf, hapPrefix = hapPrefix,
hetero_remove = TRUE, # 移除包含杂合位点的样本
na_drop = TRUE) # 移除包含基因型缺失的样本
# 从表格形式的基因型数据开始单倍型分析
hapResult <- table2hap(vcf, hapPrefix = hapPrefix,
hetero_remove = TRUE, # 移除包含杂合位点的样本
na_drop = TRUE) # 移除包含基因型缺失的样本
# 对单倍型结果进行汇总整理
hapSummary <- hap_summary(hapResult, hapPrefix = hapPrefix)
# 将单倍型鉴定结果保存到硬盘
write.hap(hapResult, file = "GeneID.hapResult")
write.hap(hapSummary, file = "GeneID.hapSummary")
# 导入之前的单倍型分析结果
hapResult <- import_hap(file = "GeneID.hapResult")
hapSummary <- import_hap(file = "GeneID.hapSummary")
2）单倍型结果展示
# 单倍型结果可视化分析
# 以表格形式展示各单倍型的基因型
plotHapTable(hapSummary, # 单倍型结果
hapPrefix = hapPrefix, # 单倍型名称前缀
angle = 45, # 物理位置的角度
displayIndelSize = 0, # 图中展示最大的Indel大小
title = geneID) # 图片标题
# 在表格中添加注释信息
plotHapTable(hapSummary, # 单倍型结果
hapPrefix = hapPrefix, # 单倍型名称前缀
angle = 45, # 物理位置的角度
displayIndelSize = 4, # 图中展示最大的Indel大小
title = geneID, # 图片标题
INFO_tag = "ANN", tag_field = 11, geneName = geneID) # 添加的注释信息
Table1.png
3) 变异位点都在哪里
# 在基因模式图上展示变异位点的信息
displayVarOnGeneModel(gff = gff, hapSummary = hapSummary,
startPOS = start-10,
endPOS = end+10,
CDS_h = 0.05, fiveUTR_h = 0.25, threeUTR_h = 0.25, # gene model parameters
cex = 0.8) # size of variants
model.png
4) 单倍型网络分析
# 单倍型网络分析
hapSummary[hapSummary == "DEL"] = "N"
hapnet <- get_hapNet(hapSummary, # 单倍型结果
AccINFO = AccINFO, # 包含样本分类信息的数据框(data.frame)
groupName = "Subpopulation", # 含有样本分类信息的列名称
na.label = "Unknown") # 未知分类样本的类别
plotHapNet(hapnet, # 单倍型网络
scale = "log2", # 标准化方法"log10"或"log2"或"none"
show.mutation = 2, # 是否展示变异位点数量, 0,1,2,3
col.link = 2, link.width = 2, # 单倍型之间连线的颜色和宽度
main = geneID, # 主标题
pie.lim = c(0.5, 2), # 圆圈的大小
legend_version = 1, # 图例形式（0或1）
labels = T, # 是否在单倍型上添加label
# legend = FALSE) # 不添加图例
# legend = TRUE) # 添加图例,但需要单击添加的位置
legend = c(12,0), # 图例的坐标
cex.legend = 0.6) # 图例中文字的大小
nework.png
5) 单倍型的地理分布
# 地理分布
AccINFO$Longitude <- as.numeric(AccINFO$Longitude)
AccINFO$Latitude <- as.numeric(AccINFO$Latitude)
hapDistribution(hapResult, # 单倍型结果
AccINFO = AccINFO, # 含有地理坐标的数据框（data.frame）
hapNames = c("H001",
"H002",
"H003"), # 展示的单倍型名称建议不超过3个
symbol.lim = c(3, 6), # 圆圈的大小
LON.col = "Longitude", # 经纬度所处的列名称
LAT.col = "Latitude", # 经纬度所处的列名称
legend = "bottomleft", # 图例所处的位置
cex.legend = 1, # 图例大小
lwd.pie = 0.2, # 圆圈线条的粗细
lwd = 1.5, # 地图线条的粗细
main = geneID) # 主标题
maps.png
6) 连锁不平衡分析
# 连锁不平衡分析
plot_LDheatmap(hap = hapResult, # 单倍型结果
add.map = TRUE, # 是否添加基因模式图
gff = gff, # 注释信息
Chr = Chr, # 染色体名称
start = start, # 基因的起始位置
end = end) # 基因的终止位置（更多参数参见帮助文档）
LDheatmap.png
7）表型关联分析
# 表型关联分析
# 单个表型的分析
hapVsPheno(hap = hapResult, # 单倍型分析结果
pheno = pheno, # 表型
hapPrefix = hapPrefix, # 单倍型名称的前缀
title = geneID, # 主标题
minAcc = 4, # 参与p值计算所需的最小样本数
symnum.args = list( # 定义显著性标注方式
cutpoints = c(0, 0.001, 0.01, 0.05, 1),
symbols = c("***", "**", "*", "ns")),
mergeFigs = TRUE) # 结果包括两个图，是否融合成一张图
# 多个表型的分析
hapVsPhenos(hap = hapResult,
pheno = pheno,
hapPrefix = hapPrefix,
title = geneID,
compression = "lzw", # tiff文件的压缩方式
res = 300, width = 12, height = 12, # 图片大小的单位"inch"
outPutSingleFile = TRUE, # 只有pdf格式支持输出单个文件
filename.surfix = "pdf", # 文件格式: pdf, png, jpg, tiff, bmp
filename.prefix = geneID) # 文件名称为: prefix + pheno_name + surfix
pheno.png
8）变异位点的效应估算
# 位点效应估算（结果仅供参考）
# EFF <- siteEFF(hapResult, pheno)
# plotEFF(EFF, gff = gff,
# Chr = Chr, start = start, end = end,
# showType = c("five_prime_UTR", "CDS", "three_prime_UTR"), # see help(plotEFF)
# y = "effect", # the means of y axis, one of effect or pvalue
# ylab = "effect", # label of y axis
# cex = 0.5, # Cex
# legend.cex = 0.8, # legend size
# main = geneID, # main title
# CDS.height = 1, # controls the height of CDS, heights of others will be half of that
# markMutants = TRUE, # mark mutants by short lines
# mutants.col = 1, mutants.type = 1, # parameters for appearance of mutants
# pch = 20) # points type
effect.png
或者逐个位点进行比较分析
# 逐位点比较变异效应
hapVsPhenoPerSite(hap = hapResult, # 单倍型分析结果
pheno = pheno, # 表型文件
phenoName = names(pheno)[10], # 表型名称
freq.min = 5) # 参与显著性计算的最小样本数
# 回车继续下一位点
# ESC退出当前命令
EFF.png
9) 什么？你的VCF文件太大，导不进来？？？
给你个小提示
filterLargeVCF()
filterLargeP.link()
4.2 好饭不怕晚，GUI操作界面
我对R语言了解比较少，学习上面方法比较困难，我该怎么办？
哈哈，不要担心，软件作者贴心的开发了对应的GUI界面，并且在软件中配有简单的使用说明，自己去探索吧！^_^
# 常规操作，加载geneHapR
library(geneHapR)
# 打开新世界的大门
startGUI.geneHapR()
# 剩下的交给你啦
五、我发现了一个bug怎么办？
在软件的gitee仓库给作者留言。或者直接联系作者，作者联系方式还是在软件的gitee仓库首页最底部^_^。
geneHapR的gitee仓库地址：(https://gitee.com/zhangrenl/genehapr)。
一定认清是gitee.com不是github.com。
六、最后的最后
如果您在您的研究中的使用了geneHapR，欢迎您引用如下文章：
Zhang, R., Jia, G. & Diao, X. geneHapR: an R package for gene haplotypic statistics and visualization. BMC Bioinformatics 24, 199 (2023). https://doi.org/10.1186/s12859-023-05318-9
