教你5分钟学会做基因单倍型分析

关键词：基因单倍型、单倍型网络图、地理分布、连锁不平衡、主效位点

怎么做单倍型分析

一、什么是单倍型？

在单倍型分析前，首先需要明白什么是单倍型、什么是基因单倍型？

一般来讲，所谓的单倍型是指同一染色体上不同变异位点的各种线性组合形式。那么基因单倍型自然就是指同一基因上（启动子、外显子、内含子、终止子）的不同变异位点的显性组合形式。

说明：基因组或染色体水平的单倍型分析不在本文范畴之内，请您参考“https://www.jianshu.com/p/4de7762aa81e”。

二、为什么要做单倍型分析？单倍型分析有什么用？

三、单倍型分析需要什么数据？数据从哪里来？

那么，进行基因单倍型分析需要那些数据呢？

1）基因型数据（必须）

既然单倍型分析需要变异位点信息，自然就需要基因型数据。基因型数据可以有很多来源，首先，可以零成本从各大公开数据库下载（推荐）；其次，可以用一代测序的方法获取基因型数据（适合于基因不太长、样本数量较少的情况）；如果你们课题组经费非常充足，可以自行设计基因芯片、高通量测序等等。

数据格式：

一代测序数据：需要拼接成序列并保存成fasta格式
二代测序数据：VCF，Hapmap，plink（map&ped）
自定义格式：表格，每行为一个变异位点，每列为一个样品，前五列固定为Chr，POS，REF，Alt和INFO

推荐一篇文章：3k水稻SNP数据集的简单利用

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_example
注释信息：Annotation_example
表型数据：Phenotype_example
其他相关数据：AccINFO_example

5.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) # 添加的注释信息

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 parameterscex = 0.8) # size of variants

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,3col.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)                # 图例中文字的大小

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)         # 主标题

6) 连锁不平衡分析

# 连锁不平衡分析
plot_LDheatmap(hap = hapResult, # 单倍型结果add.map = TRUE,  # 是否添加基因模式图gff = gff,       # 注释信息Chr = Chr,       # 染色体名称start = start,   # 基因的起始位置end = end)       # 基因的终止位置（更多参数参见帮助文档）

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, bmpfilename.prefix = geneID)            # 文件名称为: prefix + pheno_name + surfix

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

或者逐个位点进行比较分析

# 逐位点比较变异效应
hapVsPhenoPerSite(hap = hapResult,              # 单倍型分析结果pheno = pheno,                # 表型文件phenoName = names(pheno)[10], # 表型名称freq.min = 5)                 # 参与显著性计算的最小样本数
# 回车继续下一位点
# ESC退出当前命令

9) 什么？你的VCF文件太大，导不进来？？？
给你个小提示

filterLargeVCF()
filterLargeP.link()

5.2 好饭不怕晚，GUI操作界面

我对R语言了解比较少，使用上面方法比较困难，我该怎么办？
哈哈，不要担心，软件作者贴心的开发了对应的GUI界面，并且在软件中配有简单的使用说明，自己去探索吧！^_^

# 常规操作，加载geneHapR
library(geneHapR)
# 打开新世界的大门
startGUI.geneHapR()
# 剩下的交给你啦

如果有疑问可以加QQ群(255742992)询问。

六、我发现了一个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

geneHapR做基因单倍型分析相关推荐

canoco5冗余分析步骤_基因富集分析|理解
Gene Set Enrichment Analysis 基因富集分析哈罗大家好!ヾ(≧▽≦*)o 年初在和老板研究 Identifying Cell Subpopulations 有关的课题,发现 ...
基因差异表达分析——基于RSEM对比，DESeq2操作实例
以下操作全部基于R中的DESeq2函数包使用DESeq2的两点要求: DEseq2要求输入数据是由整数组成的矩阵. DESeq2要求矩阵是没有标准化的. 下面是基本流程: 一.准备工作: 确定工作路 ...
Haploview做单倍型分析
自个数据用Haploview做单倍型分析转载他人的 http://www.dxy.cn/bbs/topic/16025305 Haploview http://www.broadinstitu ...
在线作图|在线做Unifrac PCoA分析
Unifrac PCoA分析 UniFrac分析利用各样品序列间的进化信息来比较环境样品在特定的进化谱系中是否有显著的微生物群落差异.UniFrac 可用于beta 多样性的评估分析,即对样品两两之间 ...
中科院单细胞分析算法开发博士带你做单细胞转录组分析
" 福利公告:为了响应学员的学习需求,经过易生信培训团队的讨论筹备,现决定安排扩增子16S分析.宏基因组和Python课程的线上直播课.报名参加线上直播课的老师可在1年内选择参加同课程的一次 ...
R 实战 | 使用clusterProfiler进行多组基因富集分析
R 实战 | 使用clusterProfiler进行多组基因富集分析 clusterProfiler这个包我就不再介绍了,网上关于用这个包做的基础的富集分析的教程已经非常多了,今天主要介绍一下使用co ...
做生信分析平台需要什么配置的服务器？生信分析平台服务器配置建议
做生信分析平台需要什么配置的服务器? 1.CPU 2.内存 3.硬盘 4.显卡 5.不间断电源UPS 6.其它生物信息学主要研究方向:DNA/RNA/蛋白质测序,序列比对,基因发现,基因组组装,药物 ...
使用R语言包clusterProfiler做KEGG富集分析时出现的错误及解决方法
使用enrichKEGG做通路富集分析时,一直报错:显示No gene can be mapped.... k <- enrichKEGG(gene = gene, organism = &qu ...
单基因批量相关性分析 TCGA基因相关性分析单基因批量相关性分析的妙用
首先,做这个相关性分析,在这里需要安装几个R包. https://blog.csdn.net/leianuo123/article/details/102613945 https://mp.weixi ...

geneHapR做基因单倍型分析