【生信】初探基因定位和全基因组关联分析

初探qtl和gwas

实验目的

1. 了解r语言工作环境
2. 熟悉qtl和gwas分析的一般流程
3. 掌握至少一种常用qtl和gwas分析软件

实验内容

1. 使用r包qtl2进行qtl分析
2. 使用软件tassel进行gwas分析
3. 熟悉基因型、表型等数据格式，掌握qtl和gwas分析结果的解读和可视化

实验题目

第一题：玉米magic群体的qtl分析

利用玉米magic群体数据，使用r/qtl2完成以下qtl分析：

① 本群体包含哪些数据？

② 绘制株高（ph）的lod曲线

③ 给出株高对应的qtl

第二题：tassel自带数据集的关联分析

使用tassel对自带数据集（安装目录下tutorialdata子目录中），对其它两个性状（earht和eardia）中任意一个，进行关联分析。

至少需要给出以下结果：

①与earht或eardia最显著关联的top 5位点信息（包括染色体号、位置和p值）

②两张图（曼哈顿图、qq图）及相应的解释（绘图结果说明什么）

实验过程

玉米magic群体的qtl分析

① 包含的数据

maize_magic_geno.csv: magic的基因型数据

maize_magic_foundergeno.csv: 9个founder系的基因型数据

maize_magic_gmap.csv: 遗传图谱

maize_magic_pmap.csv: 物理图谱

maize_magic_crossinfo.csv: cross信息(代数后面跟着表示九个founder相对贡献的值)

maize_magic_pheno.csv: 表型数据

maize_magic_phenocovar.csv: 表现型协变量，只描述表现型

② 绘制lod曲线

我们使用r语言的qtl2包，进行qtl定位分析。qtl2的使用方法可以参考官方文档，写的很不错。

下面就是代码实现了，注释写的也挺清楚的。

# 设置工作目录
setwd("d:/00大三上/生信原理/实践作业/qtl_gwas")

# install.packages("qtl2")
library(qtl2)
# 直接从网站上读取
file <- paste0("https://raw.githubusercontent.com/rqtl/","qtl2data/main/maizemagic/maize_magic.zip")

mm <- read_cross2(file)
# 或者把压缩包下载下来再读取
mm <- read_cross2("./maize_magic.zip")

# 观察数据 - 10条染色体、总计41324个标记
summary(mm)

# （伪）标记插入遗传图，获取假定qtl；以1cm为间隔插入伪标记
mmmap <- insert_pseudomarkers(map=mm$gmap, step=1)

# 计算基因型概率;假定基因分型误差概率0.002
# cores=4,使用多核计算
mmpr <- calc_genoprob(cross = mm, map = mmmap, error_prob=0.002, cores=4)

# 可视化查看基因型概率
# 参数依次为 基因型概率、marker图、要查看的个体编号、要查看的染色体号
# 染色体的坐标在横轴上，基因型在纵轴上。较高的基因型概率表示为暗色
png(filename = 'plot_genoprob.png',width = 3580,height = 2200,res = 300)
plot_genoprob(mmpr, mmmap, ind = 1, chr = 1)
dev.off()

# 运用 haley-knott regression 进行基因组扫描
# 可加协变量、此处不加。输出lod分数矩阵（positions × phenotypes）
mm_scan1_out <- scan1(mmpr, mm$pheno, cores=0)

# 绘制lod曲线，指定一列（表型）
png(filename = 'plot_scan1.png',width = 3580,height = 2200,res = 300)
plot_scan1(mm_scan1_out, map = mmmap, lodcolumn = "ph")
dev.off()
# 查看对ph而言哪个伪标记lod分数最高、哪个基因型标记lod得分最高
sort(mm_scan1_out[,"ph"])

# permutation test 说明scan结果的统计学意义
# 识别随机下可能出现的最大lod分数，使用1000种排列
# 这一行执行时间比较长，因为置换检验需要重新计算1000次
mm_operm <- scan1perm(genoprobs = mmpr, pheno = mm$pheno, n_perm = 1000) 

# 显著性阈值、默认5%水平
# ph为7.40 期望lod得分低于7.40是偶然事件
summary(mm_operm, alpha=0.05) 

# 寻峰，95%置信区间
thr <- summary(mm_operm)
mm_peaks <- find_peaks(scan1_output = mm_scan1_out, map = mmmap, threshold = thr, prob = 0.95, expand2markers = false)

# 查看ph表型对应多少满足阈值的峰（qtl），分别在哪
ph_qtl <- mm_peaks[mm_peaks$lodcolumn=='ph',]

# 保存结果，方便之后查看
write.csv(thr,"thr.csv")
write.csv(ph_qtl,"ph_qtl.csv")

图1 一号染色体基因型概率可视化

图2 lod曲线

由lod曲线可以看出，6、8和10号染色体很有可能有qtl位点，这也与下文的分析对应上了。

株高对应的qtl

有上面代码的输出结果（表1、表2）可以看到，ph的阈值是7.41，筛选过后可以得到三个qtl。分别位于6、8、10号染色体上，这也与上面的lod曲线对应上了。

表1 四种表型的阈值

	ps	ph	eh	gyrad
0.05	7.424502	7.408967	7.753835	7.628907

表2 株高对应的qtl

	lodcolumn	chr	pos	lod	ci_lo	ci_hi
4	ph	6	12.28468	15.91047	11.0098	12.3023
5	ph	8	90.14981	10.9792	79.44594	93.54998
6	ph	10	13.94991	7.916526	9.585541	43.14987

tassel自带数据集的关联分析

tassel简介

trait analysis by association, evolution and linkage.

tassel是一个软件包，用于评估特征的关联，进化模式和连锁不平衡。本软件的优点包括:

1. 有机会使用一些新的强大的统计方法来进行关联映射，例如通用线性模型(glm)和混合线性模型(mlm)。mlm是《nature genetics》论文—— unified mixed-model method for association mapping，该技术减少了与复杂谱系、家族、创始效应和种群结构关联映射中的i型错误。

2. 能够处理广泛的索引(插入和删除)。大多数软件忽略了这种类型的多态性;然而，在某些物种(如玉米)中，这是最常见的多态性类型。

实际操作

下载安装完tassel后，打开安装目录下tassel5\tutorialdata中的示例数据。

分别导入：

• 基因型数据mdp_genotype.hmp
• 群体结构数据 mdp_population_strucutre
• 表型性状数据 mdp_phenotype、mdp_traits

图3 导入示例文件

菜单栏 的filter可以基于每个位点的基因型统计分析，进行筛选。

我们通过maf来进行筛选。

minor allele frequency (maf)：次等位基因频率

比如：某个snp位点，包含a和g两种等位基因（allele)，100个个体中有35个为a,65个为g,则该位点的maf=0.35

过滤完成后，软件界面左侧会出现文件mdp_genotype_filter

图4 数据过滤

图5 数据过滤参数设置

数据过滤后，我们计算群体内不同个体之间的亲缘关系，可用于后续混合线性模型mlm的关联分析。

图6 计算亲缘关系

选择默认的kinship计算方法，center ibs和max alleles使用默认值。

图7 计算亲缘关系的参数设置

下表是有关参数设置的帮助。

表3 计算亲缘关系参数设置的帮助

parameter	description	values	default
kinship method	the centered_ibs (endelman - previously scaled_ibs) method produces a kinship matrix that is scaled to give a reasonable estimate of additive genetic variance. uses algorithm http://www.g3journal.org/content/2/11/1405.full.pdf equation-13. the normalized_ibs (previously gcta) uses the algorithm published here: http://www.ncbi.nlm.nih.gov/pmc/articles/pmc3014363/pdf/main.pdf.	centered _ibs, normalized _ibs, dominance _centered _ibs, dominance _normalized _ibs	centered _ibs
max alleles	max alleles	2…6	6
algorithm variation	algorithm variation	observed _allele _freq, proportion _heterozygous	observed _allele _freq

计算完成后，生成亲缘关系矩阵。

图8 亲缘关系矩阵 下面我们选择一个形状进行关联分析，在这里选择earht

图9 性状筛选

图10 选择特定性状

去掉最后一个群体结构（q3），指导手册给出的解释说，如果我们把它们全部作为协变量使用，这会产生线性相依性。

图11 去掉最后一个群体结构

摁住ctrl键同时选中上述三个文件（ mdp_genotype_filter，filtered_mdp_traits ，filtered_mdp_population_structure ）进行合并, 点击data → intersect join

图12 选择需要合并的数据

图13 合并数据

使用glm进行关联分析（默认参数）。

图14 使用glm进行关联分析

图15 glm参数设置

曼哈顿图和qq图

得到结果之后，我们就可以产看结果了，在菜单栏里的results里面有画曼哈顿图和qq图的选项。

曼哈顿图中每个点代表一个snp，纵轴为每个snp计算出来的pvalue取-log10，横轴为snp所在的染色体。基因位点的pvalue越小即-log10(pvalue)越大，其与表型性状或疾病等关联程度越强。

可以发现1、8和10号染色体与earth的性状可能有关联，从qq图上也可以看到，偏离了对角线，有一个凸起，表明有些标记可能是与earth性状相关。

图16 glm的曼哈顿图

图17 glm的qq图

接下来我们再试试混合线性模型，操作和glm差不多。

图18 使用mlm进行关联分析

从qq图来看，和性状关联的位点不多，故后面分析用glm模型的结果。

图19 mlm的曼哈顿图

图20 mlm的qq图

top5位点

我们将p值从小到大排序则可以看到最显著关联的5个位点。

图21 筛选p值

导出之后的表格如下

表4 最显著关联的top 5位点信息

chr	pos	p
10	33619582	0.000016724
1	245136244	0.000021302
4	194749287	0.000028892
4	194749230	0.000048693
8	134813437	0.000052242

讨论

1.qtl分析和gwas分析的一般流程是什么？

qtl分析的一般流程包括以下步骤：

1. 构建分离群体（作图群体）
2. 获取每个个体的目标性状的表型数据
3. 对作图群体进行基因分型：获取个体基因型数据、筛选分子标记
4. 使用统计方法检测qtl

qtl分析的一般流程包括以下步骤：

1. 自然群体资源收集和鉴定
2. 获取目标性状的表型数据
3. 群体重测序（全基因组、转录组、外显子组、gbs等），获取基因型数据（测序读段比对、遗传变异如snp的识别、gwas变异标记（如snp）构建）
4. 基因型补缺和过滤，表型数据分析
5. 遗传多样性、群体结构、亲缘关系、ld等分析
6. 关联分析
7. 候选基因筛选与挖掘

2.qtl分析中，如何确定目标性状的候选遗传位点？

在qtl分析中，目标性状的候选遗传位点通常使用统计方法来识别，该方法分析个体样本中遗传变异和性状之间的关系。

有几种方法可用于识别候选qtl，包括:

1. 连锁分析:该方法利用家族中遗传标记(如snp)的遗传模式来识别可能与性状相关的位点。
2. 关联分析:这种方法在不相关的个体样本中测试遗传变异和性状之间的统计关联。
3. 混合模型分析:该方法结合了连锁和关联分析的元素，以解释可能影响性状的已知和未知因素的影响。

一旦确定了候选qtl，就可以使用其他数据和方法进一步验证它们，例如在不同的样本中重复研究或使用功能测定来确定关联背后的生物学机制。考虑潜在的混杂因素也很重要，比如环境或生活方式因素，这些因素可能会影响这一特征。

3.gwas分析中，如何确定目标性状的候选遗传位点？如何缩小候选区间？

在gwas分析中，通常使用统计方法分析遗传变异与样本中个体的目标特征之间的关系，从而确定目标特征的候选遗传位点（）。

可以使用以下几种方法来确定候选gwas位点：

1. 单标记分析：对于每个单个的位点（例如snp）进行统计分析，以识别与特征关联的位点。
2. 多标记分析：对于一组位点进行统计分析，以识别与特征关联的位点。
3. 基于统计模型的分析：使用统计模型来考虑可能影响特征的各种因素（包括环境和生活方式因素）。

确定候选gwas位点后，可以使用附加数据和方法进一步验证它们，例如在不同样本中复制研究或使用功能性分析来确定关联的生物学机制。

要缩小候选区间，可以使用多种方法，包括：

1. 增加样本大小：使用更多的个体来提高统计力，使候选区间更加明确。
2. 使用高分辨率基因分型技术：使用高分辨率基因分型技术（如下一代测序）可以提供更多的遗传变异数据，从而帮助缩小候选区间。
3. 引入附加信息：使用其他数据（如生物学或功能性分析结果）来帮助缩小候选区间。
4. 复制研究：在不同的样本中复制研究可以帮助确定候选区间的稳定性。
5. 使用多种统计方法：使用多种不同的统计方法可以帮助缩小候选区间，因为不同的方法可能会得出不同的结论。

4.如何降低qtl、gwas分析结果的假阳性？

选择重组事件多、变异丰富、群体分化不大的群体。并且选择合适的模型。
候选gwas位点后，可以使用附加数据和方法进一步验证它们，例如在不同样本中复制研究或使用功能性分析来确定关联的生物学机制。

要缩小候选区间，可以使用多种方法，包括：

1. 增加样本大小：使用更多的个体来提高统计力，使候选区间更加明确。
2. 使用高分辨率基因分型技术：使用高分辨率基因分型技术（如下一代测序）可以提供更多的遗传变异数据，从而帮助缩小候选区间。
3. 引入附加信息：使用其他数据（如生物学或功能性分析结果）来帮助缩小候选区间。
4. 复制研究：在不同的样本中复制研究可以帮助确定候选区间的稳定性。
5. 使用多种统计方法：使用多种不同的统计方法可以帮助缩小候选区间，因为不同的方法可能会得出不同的结论。

4.如何降低qtl、gwas分析结果的假阳性？

选择重组事件多、变异丰富、群体分化不大的群体。并且选择合适的模型。