群体进化-gwas分析

在这里插入图片描述

群体进化基础分析

PCA

• 分析原理
  • PCA(Principal Component Analysis)，即主成分分析方法，是一种使用最广泛的数据降维算法。PCA的主要思想是将n维特征映射到k维上，这k维是全新的正交特征也被称为主成分，是在原有n维特征的基础上重新构造出来的k维特征。PCA的工作就是从原始的空间中顺序地找一组相互正交的坐标轴，新的坐标轴的选择与数据本身是密切相关的。其中，第一个新坐标轴选择是原始数据中方差最大的方向，第二个新坐标轴选取是与第一个坐标轴正交的平面中使得方差最大的，第三个轴是与第1,2个轴正交的平面中方差最大的。依次类推，可以得到n个这样的坐标轴。通过这种方式获得的新的坐标轴，我们发现，大部分方差都包含在前面k个坐标轴中，后面的坐标轴所含的方差几乎为0。于是，我们可以忽略余下的坐标轴，只保留前面k个含有绝大部分方差的坐标轴。事实上，这相当于只保留包含绝大部分方差的维度特征，而忽略包含方差几乎为0的特征维度，实现对数据特征的降维处理。
  • 简洁点来讲现在有这样的数据，100个样品，2M标记，即是2000000X100的矩阵，那么就通过数学降维的方法简化到100X3甚至100X2乘（即PC1，PC2）
• 分析软件
  • GCTA
  • tassel
• EIGENSTRAT
• 结果展示
  • PCA结果矩阵（特征向量）

GWAS_1  0.0295707   0.0174155   -0.0245656
GWAS_10 0.0212291   -0.0552983  -0.0280335
GWAS_100    -0.0645872  0.00456635  0.00588907
GWAS_101    -0.0779853  -0.0317529  0.0138288
GWAS_102    -0.0790227  -0.0295285  0.0147819
GWAS_105    -0.0845384  0.000685319 0.0108059
GWAS_108    -0.0779536  -0.00380985 0.0101755
GWAS_109    -0.0789908  -0.00534946 0.012742
GWAS_11 0.0152839   0.0185823   -0.0305629
GWAS_110    -0.080786   -0.00255263 0.0131448

    *   第一列样品名称，第二列PC1的值，第三列PC2的值，第四列PC3的值（也就是平时看到的结果图的横纵坐标来源）
*   PCA解释数据结果（特征值）

54.402
32.2402
25.6809
18.0063
13.7968
9.6096
9.46086
9.00158
8.16587
7.60115

    *   这个结果每一个值对应的维度的解释情况，行数与样品数量一致，第一行代表第一维，依次类推；每一行除以所有行数的和即是其第几维解释的比例

• PCA结果图
• 在这里插入图片描述
  在这里插入图片描述
  在这里插入图片描述
  在这里插入图片描述

• 结果要点
  • 结果图中相对分群结果合理，大部分跟生产经验相一致
  • 解释度可接受，这个方面想了解的话，可以看看文献的，不是硬标准

TREE

• 分析原理

  • 系统发育分析中，最重要和最常见的内容为构建系统发育树。系统发育树也称为系统发生树（phylogenetic tree）、聚类树或者进化树（evolutionary tree）。以树状结构表示各个节点的进化关系，枝点可以是物种、同一物种的样本、基因等单元。
  • 根据SNP或者Indel 构建其系统进化树，可以展示群体中不同个体的相互关系，基因变异相似的往往会在同一个树的cluster中，一颗好的树可以给你一个群体大概的分类(你这个群体中有多少个cluster，一般同一个亚种或者有亲缘关系的个体会形成一个cluster），这是群体遗传中重要的一部分。其构建的核心原理就是把每个位点SNPs的信息提取，然后计算每个变异位点的差异得到算法中的“距离”。

• 分析软件

  • treebest
  • mega
  • taseel
  • snphylo
  • phylip

• 算法

  • nj 临近算法
  • ml 最大似然值

• 结果展示

  
  在这里插入图片描述
  在这里插入图片描述

• 结果要点
  • 分群清晰,大部分样品与已知生产经验一致

STructure

• 分析原理

  • 先预设群体由若干亚群（k=x）构成，通过模拟算法找出在k=x的情况下，最合理的样本分类方法。最后再根据每次模拟的最大似然值，找出最适用这群体的K值。

• 分析软件

  • STRUCTURE
  • ADMIXTURE
  • FRAPPE
  • fast-STRUCTURE

• 算法

  • 亚群内符合哈温平衡
    • 那么，软件在如何确定样本的最优分类方法呢？其实基于一个假设：在各个亚群内部个体应该符合哈代-温伯格平衡（哈温平衡的概念可以在百度查询），那么这个亚群内的基因频率分布应该可通过哈温平衡检验。例如，现在有40个个体的1个SNP位点的基因型，我预设亚群数k=2。我先随机将40个个体分成两份，然后检验是否符合哈温平衡。如果不符合，我继续调整分类策略，直到找到一种最优的分类方法：40个个体被分为了两份，每个亚群都由若干个体构成，每个亚群内部都最大程度地符合哈温平衡。
  • 每个位点是独立的
    • 同一个体基因组上的不同SNP可能来源不同亚群体,软件是对每个位点单独进行分群的,然后统计比例,所以要求进行分析的每个位点是独立,不然会造成比例的不准确
  • 每个样本的血统构成
    • 以k=2为例，解释一下structure是如何找到样本的最优分类。其实简单说来，就是利用了计算机超强的运行能力，一开始计算机只是随机将样本分为两份，然后在每个亚群内进行哈温平衡检验。如果不符合哈温平衡（拍脑袋的分类，一开始当然是惨不忍睹），计算机继续调整分类，然后继续检验。
    • 如此这般，在计算n次后，计算机再从这一堆结果中找到最佳的分类。这个过程称为“隐马科夫-蒙特卡罗链”的过程，计算次数n就是这个链的长度，这是structure一个重要的参数“Number of MCMC Reps”，需要预先设定。
    • 但因为这个计算的过程是从随机模拟开始的。如果一开始拍脑袋拍的不好（随机分类与真实分类差距太大），计算机一黑到底，最后把n次用完了，都没有找到一个合理的分类。所以，分析软件往往有个预实验的过程。
    • 就是在正式进行大规模运算前，计算机先尝试各种各样的随机分类，运行非常短的次数，然后评估哪种随机分类是最合理的。之后，在根据最优的随机分类，进行后续的大规模运算。这个过程就称为burn-in period，预实验的次数就称为burin-in的次数。这也是structure分析另外一个重要的参数“length of burn-in period”。
  • 选择使用那种模型
    • 主要涉及两种模型 no admixture model和admixture model。前者假设亚群间不存在杂交，后者则假设亚群间存在杂交。在绝大部分情况下，当然是选择admixture 模型更合理了。

• 结果展示

  
  在这里插入图片描述

• 结果要点
  • 最终k值选取的结果跟实际生产相符性
  • 一般会以该结果的0.8或者0.6以上血缘比例的作为标准(血缘复杂的根据情况有些可以追溯原因)

LD

• 分析原理
  • 只要两个基因不是完全独立遗传，就会表现出某种程度的连锁，这种情况就叫连锁不平衡。
  • 由于HLA不同基因座某些基因经常连锁在一起遗传，而连锁的基因并非完全的随机的组成单体型，有些基因总是较多的在一起出现，致使某些单体型在群体中呈现较高的频率，从而引起连锁不平衡。
  • 假如位于同一染色体的两个等位基因（AB）同时存在的概率大于人群中因随机分布而同时出现的概率，称这两点处于LD状态。
  • LD的基本单位是D，但是度量观察到的单倍型频率与平衡状态下期望频率的偏差。
  • 虽然D能够很好的表达LD的基本含义，但是由于其严格依赖于等位基因频率（allele frequency），故不适合应用于表述实际的LD强度。
  • 所以一般在LD的度量中最常见的是D'和r2。二者各有各的特点和用途，但都是基于D的。
    • 当D'=0，r2=0时，处于完全连锁平衡状态
    • 当D'=1，r2=1时，处于完全连锁不平衡状态。
    • 其中，从0—1之间的度量越高，LD越高，如果两个位点连锁，连锁程度也越高。
• 分析软件
  • plink
  • haploview
• 算法
  • 1、设有两个位点（A、B），等位基因分别是A、a、B、b，在群体中对应频率f(A)、f(a)、f(B)和f(b)
  • 2、两个位点共有四种单倍型AB、Ab、aB、ab，对应频率f(AB)、f(Ab)、f(aB)和f(ab)
  • 3、计算：Dab=f(AB)-f(A)*f(B)
    • 当Dab=0时，处于连锁平衡状态；
    • 当Dab≠0时，处于连锁不平衡状态。
  • LD度量：
    • 当Dab>0，|D'|=(Dab)2/min(f(AB), f(ab))；
    • 当Dab<0，|D'|=(Dab)2/min(f(Ab), f(aB))；
    • r2=(Dab)2/(f(A)f(a)f(B)*f(b))；

• 结果展示

  
  在这里插入图片描述
  在这里插入图片描述
  在这里插入图片描述

• 结果要点
  • 一般如果是GWAS项目,没有明显的分群的话,会做这个作为一个平均连锁距离的依据(文章中大多采用两种策略:1.LD的R2降到最高点的一半时的距离作为平均连锁距离;2.LD的R2降到0.2时的距离作为平均连锁距离)
  • 如果是群体进化项目,会有非常明确的群体划分,各自群体分别进行LD分析,能够从一定程度上讨论进化快慢等信息
  • 倒三角的具体区段的连锁图一般结合具体的GWAS等的位点一块展示，这个是不可能全基因组绘制的

GWAS

• 分析原理
  • 基本思想：应用于复杂性状，采用CASE-CONTROL试验设计，比较全基因组范围内所有SNP位点的等位基因或者基因型频率在case与control组中的差异,如果某个snp位点等位基因或者等位基因型case组中的频率明显高于或者低于caontrol组,则认为该位点与该疾病间存在关联
• 分析软件
  • MVP
  • GEmma
  • plink
  • taseel
  • GAPIT
  • FarmCpu
• 算法
  • LM
  • MLM
  • FarmCpu
• 方案设计要点
  • 群体大小 >300
  • 群体分层不明显(具有一致的遗传背景或者资源的群体)
  • 覆盖全基因组的高密度的标记(至少保证平均一个block内有一个标记,中玉金标记公司内现有芯片只有660k符合)
  • 表型数据记录准确性
  • 植物数据尽量是多年多点的数据
  • 表型分布比较广泛(大体成正太分布)

• 结果展示

  
  在这里插入图片描述
  在这里插入图片描述
  在这里插入图片描述
  在这里插入图片描述
  在这里插入图片描述

• 结果要点

群体进化高级分析-群体选择消除分析

pi

• 分析原理

  • π用来分析碱基多态性，多态性越低，受选择程度越高。

• 分析软件

  • vcftools

• 算法

• 结果展示

• 结果要点

fst

• 分析原理
  • 群体的固定系数F反映了群体等位基因杂合性水平。固定系数F是F统计量（Fst）的一个特例。Fst分析表示群体的分化程度，值越大，群体分化程度越高，受选择程度越高。
• 分析软件
  • vcftools
• 算法
• 结果展示

在这里插入图片描述

• 结果要点

XP-LXR

• 分析原理
• 分析软件
• 算法
• 结果展示
• 结果要点

D

w

ROH

群体进化高级分析-种群动态等

PSMC

• 分析原理
• 分析软件
• 算法

• 结果展示

  
  在这里插入图片描述

• 结果要点

Treemix

• 分析原理
• 分析软件
• 算法

• 结果展示

  
  在这里插入图片描述

• 结果要点

案例解读

群体进化

案例一:NG-熊猫群体进化-2012

• 文章:Whole-genome sequencing of giant pandas provides insights into demographic history and local adaptation
• 基础数据:34只熊猫,4.7x覆盖深度
• 分析结果:

  • 群体结构分析

    
    在这里插入图片描述

* 种群历史动态分析

在这里插入图片描述

* 选择消除分析结果注释情况
    * 该项目使用fst进行选择消除分析,分析后受选择的基因进行KEGG富集分析

在这里插入图片描述

案例二:NC-牦牛群体进化-2015

• 文章:Yak whole-genome resequencing reveals domestication signatures and prehistoric population expansions

• 基础数据:13野生牦牛和59驯化品种,6.7X测序深度,14.56M高质量SNP

• 分析结果

  • 群体结构分析

    
    在这里插入图片描述

* 选择消除分析

在这里插入图片描述

* 种群历史动态分析 && 统计分析
    * 

在这里插入图片描述

gwas

案例一:NG-韩斌2010年经典水稻14农艺性状GWAS文章

• 文章:Genome-wide association studies of 14 agronomic traits in rice landraces
• 基础数据:517水稻样品,3.6M SNP,水稻indica,japonica

• 看下文章的作者

  
  在这里插入图片描述

• 分析结果

  • tree && LD && maf

    
    在这里插入图片描述

* tree && PCA

在这里插入图片描述

* imputation 准确性评估(2010年测序成本还很高,所以进行的地深度测序,进行缺失推断)
    * 这里评估了LD,测序错误,测序覆盖情况,样品数量对推断结果的影响

在这里插入图片描述

* GWAS结果曼哈顿图和QQ图展示
    * 使用了两种算法,一般线性模型和混合线性模型

在这里插入图片描述

* gwas关联定位情况统计表格
    * 展示基本的性状,定位染色体,位点,两种基因型,最小等位基因频率,p值以及前期研究的积累的情况

在这里插入图片描述

* GWAS结果曼哈顿图局部展示与基因结构展示

在这里插入图片描述

* 定位位点及数量统计
    * 其实是有对比两种算法的结果

在这里插入图片描述

案例二:NG-韩斌2011年水稻开花期与果实性状GWAS文章

• 文章:Genome-wide association study of flowering time and grain yield traits in a worldwide collection of rice germplasm
• 基础数据:950水稻样品,来源于33个国家,4.1M snp
• 数据过滤:maf 0.05
• 为什么可以很快又发一篇NG
  • 样品来源范围及样品数量扩大
  • 性状改变

  • 重点进行了结果单体型的分析

    
    在这里插入图片描述

• 分析结果
  • tree && fst

    • 相对于2010年的文章,该次进行了群体间fst分析

      
      在这里插入图片描述

* 6个已定位基因的关键变异信息

2011-hanbin-rice-3.png

* 单体型的结果展示

在这里插入图片描述

* 单体型的数据统计

在这里插入图片描述

* GWAS结果曼哈顿图和QQ图
    * 这次可以看到没有一般线性模型了,那是因为上篇比较过,没必要再进行比较了

在这里插入图片描述

* GWAS定位基因情况
    * 本篇中比较简洁,性状,染色体,物理位置,基因,基因的描述

在这里插入图片描述

* 局部定位结果及基因结构情况
    * 注意,这里加入了表达量的情况

在这里插入图片描述

案例三:NG-日本2016年水稻开花相关基因性状GWAS分析

• 文章:Genome-wide association study using whole-genome sequencing rapidly identifies new genes influencing agronomic traits in rice
• 基础数据:176japonica,5.8X,383g,426k snp,67k indel
• 过滤maf:0.05
• 分析结果
  • 表型数据信息,以及176样品能够代表413样品(我记得该文章是从之前文章的数据中拿了176样品进行的GWAS分析)

在这里插入图片描述

* 表型处理后分布统计情况及PCA展示群体情况

在这里插入图片描述

* 不同性状GWAS分析结果情况
    * 文章证据充分:
    1. 做了转基因(转的不同的单体型)的对照实验
    2. 单体型分析到位:包含snp和indel
    3. 对gwas定位结果的分级,为后期验证的顺序有指导意义

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述

案例四:NG-棉花-纤维相关性状GWAS分析

• 文章:Resequencing a core collection of upland cotton identifies genomic variation and loci influencing fiber quality and yield

• 基础数据:419样品,6.55X数据,3.66M snp,13个性状

• 第一批棉花gwas项目,对A,D基因组差异进行了分析,环境有12个,有相应的拟南芥过表达表型验证

• 分析结果

  • 数据注释统计

    
    在这里插入图片描述

* 群体结构分析

在这里插入图片描述

* 亚群多态性和LD分析

在这里插入图片描述

* GWAS分析结果

在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述
在这里插入图片描述 在这里插入图片描述

群体进化+GWAS

案例一:NBT-田志喜-大豆-2014群体gwas文章

• 文章:resequencing 302 wild and cultivated accessions identifies genes related to domestication and improvement in soybean
• 基础数据:302个野生,栽培,地方品种大豆,11X,
• 分析结果

  • 样品来源与群体结构

    
    在这里插入图片描述

* LD分析

在这里插入图片描述

* 选择消除分析与GWAS分析结果

在这里插入图片描述
在这里插入图片描述

* fst统计分析

在这里插入图片描述
在这里插入图片描述