全基因组关联分析数据分析流程
in 生信工具 with 0 comment

全基因组关联分析数据分析流程

in 生信工具 with 0 comment

全基因组关联分析(genome-wide association study, GWAS)的数据分析流程的一般步骤为:数据比对、call SNP鉴定基因型、表型统计和基因型表型关联分析。由于重测序数据质控和比对、变异位点鉴定和分型两部分过于简单,在此略去。下面主要对后续的基因型过滤、群体结构、亲缘关系、LD衰减分析及GWAS具体步骤进行说明。

基因型过滤

为了保证关联分析的计算效率和统计学效力,一般使用VCFtools或PLINK过滤缺失率高和次要等位频率较低的SNP。而对于是否还需要根据LD进行过滤,我认为没有必要,不过滤对于最后结果的影响不大,无非就是加快运算速度。因此,下面对基因型将不进行LD过滤。

使用VCFtools过滤

vcftools --vcf [vcf] [--plink] --max-missing 0.8 --maf 0.05 [--remove-indels] --out [outfile] 

使用PLINK过滤

plink --vcf [vcf] --maf 0.05 --geno 0.1 --recode vcf-iid --out [outfile]

然后可使用snpEff对SNP进行注释,将SNP按其在基因组上的相对位置进行分类,包括基因上游、5‘-UTR区、外显子区、内含子区、3’-UTR区和基因下游等。 同时,注释SNP对蛋白产物的影响,如同义突变、非同义突变、移码突变和终止密码子提前等。

java -jar snpEff.jar ann [database] [vcf] -v -ud 5000 -t >[outfile]

群体结构、亲缘关系和LD衰减分析

为了降低群体结构和家系亲缘关系对全基因组关联分析的影响,需利用SNP信息计算出代表群体结构的Q矩阵和家系亲缘关系的K矩阵。

PCA分析确定主成分来控制群体结构,对群体结构进行检验和校正,主成分得分信息还可用于关联分析的混合线性模型中,从而减少群体结构带来的假阳性关联。

PCA分析

这里使用plink和GCTA进行分析,首先需要将vcf转为ped和bed格式。

vcftools --vcf [vcf] --plink --out [outfile]
plink --noweb --file [vcf] --make-bed --out [outfile]

plink计算PCA

这里只计算前三个PC,因一般也只看前三个。

plink --bfile [file_prefix] --pca 3 --out [outfile]

该步骤会生成2个文件,一个是eigenval结尾的文件,记录特征值,用来计算每个PC所占的比值。另一个是eigenvec结尾的文件,记录特征向量。

GCTA计算PCA

首先计算近交系数,生成grm矩阵:

gcta64 --bfile [file_prefix] --make-grm --autosome --out [outfile]

然后计算PCA

gcta64 --grm [grm_file] --pca 3 --out [outfile]

同样生成以.eigenval和.eigenvec结尾的文件.

PCA分析一般保留1-10个主成分来完成GWAS关联分析中混杂因素的矫正,一般选取3个主成分来做后续分析。

群体结构推断

为了了解群体遗传构成,一般使用admixture进行群体结构推断,可以添加参数j(多线程)来加快运算速度。

for i in {1..15}
do
admixture --cv [ped_file] ${i} -j50 >> log.txt
done

一般选取CV值最小时的K值作为最佳分群,其Q矩阵用于后续分析。

亲缘关系分析

亲缘关系K矩阵可以通过EMMAX计算获得:

emmax-kin-inter64 [tped_file] -v -d 10

计算完成后会输出一个aBN.kinf结尾的文件,即为后续分析使用的K矩阵。

LD衰减分析

LD衰减分析常用的软件有plink、Haploview、PopLDdecay.

plink计算LD

plink --file [name] --r2 --ld-window 99999 --ld-window-r2 0 --ld-window-kb 1000 --out [fileout]

PopLDdecay计算LD

PopLDdecay -InVCF [name].vcf.gz -OutStat [name].LD

LD半衰距离为连锁不平衡参数r2衰减至最大值的一半时对应的距离。实践中常用该值来评估群体中遗传标记连锁与重组情况,确定关联分析所需标记密度及GWAS结果中的显著信号在基因组候选基因的选取范围。

GWAS具体步骤

GWAS一般以群体结构和亲缘关系矩阵作为协变量,通过某种模型将SNP与表型关联起来。

现以EMMAX、GAPIT和rMVP为例进行说明。

利用EMMAX进行GWAS

EMMAX的运行需要准备的文件包括基因型文件、亲缘关系矩阵、群体结构协变量文件和表型文件。

基因型文件

基因型文件一般直接使用VCF过滤后得到,具体如下(这里略去了将染色体名改为数值和添加SNP ID两步)。

按次要等位频率(0.05)和缺失率(0.1)进行过滤

plink --vcf 299samples.SNP.recode.onlyChr.renameChr.vcf --maf 0.05 --geno 0.1 --recode vcf-iid --out 299samples.SNP.recode.onlyChr.renameChr.maf0.05.int0.9

亲缘关系矩阵

emmax-kin-intel64 299samples.SNP.recode.onlyChr.renameChr.maf0.05.int0.9 -v -d 10

群体结构协变量文件

  1. 将上述过滤的基因型文件转为admixture需要的格式
plink --vcf 299samples.SNP.recode.onlyChr.renameChr.maf0.05.int0.9.vcf --recode 12 --out 299samples.SNP.recode.onlyChr.renameChr.maf0.05.int0.9
  1. 使用admixture进行群体结构分析
for i in {1..15}
do
admixture --cv 299samples.SNP.recode.onlyChr.renameChr.maf0.05.int0.9.ped ${i} -j50 >> log.txt
done
  1. 查看最佳分群
grep CV log.txt
CV error (K=1): 0.78595
CV error (K=2): 0.73205
CV error (K=3): 0.69132
CV error (K=4): 0.67968
CV error (K=5): 0.66617
CV error (K=6): 0.66016
CV error (K=7): 0.64714
CV error (K=8): 0.64350
CV error (K=9): 0.64214
CV error (K=10): 0.63923
CV error (K=11): 0.64538
CV error (K=12): 0.63823
CV error (K=13): 0.64402
CV error (K=14): 0.64683
CV error (K=15): 0.64720

K=12时CV值最小,所以最佳分群为12 K=12时候的Q矩阵用于后续分析。admixture产生的Q矩阵格式如下(前5行)

0.000010 0.000010 0.016867 0.000010 0.227751 0.072513 0.065477 0.239721 0.014277 0.014286 0.216792 0.132285
0.000010 0.123528 0.096458 0.000010 0.112097 0.066272 0.333813 0.124592 0.073127 0.000010 0.030636 0.039448
0.000010 0.000010 0.041969 0.009214 0.076072 0.105772 0.164234 0.169847 0.017128 0.000010 0.415723 0.000010
0.009623 0.207993 0.000010 0.000017 0.000010 0.125707 0.049496 0.034851 0.201822 0.000012 0.000010 0.370449
0.113456 0.000010 0.000010 0.000010 0.000010 0.000010 0.000010 0.711426 0.162180 0.000010 0.000010 0.012858

因文件较小,可以手动在excel中进行格式调整(删除最后一列),调整后格式如下:

TH001   TH001   1       0.00001 0.00001 0.016867        0.00001 0.227751        0.072513        0.065477        0.239721        0.014277        0.014286        0.216792
TH002   TH002   1       0.00001 0.123528        0.096458        0.00001 0.112097        0.066272        0.333813        0.124592        0.073127        0.00001 0.030636
TH003   TH003   1       0.00001 0.00001 0.041969        0.009214        0.076072        0.105772        0.164234        0.169847        0.017128        0.00001 0.415723
TH004   TH004   1       0.009623        0.207993        0.00001 0.000017        0.00001 0.125707        0.049496        0.034851        0.201822        0.000012        0.00001
TH005   TH005   1       0.113456        0.00001 0.00001 0.00001 0.00001 0.00001 0.00001 0.711426        0.16218 0.00001 0.00001

将此文件保存为emmax_cov.txt作为协变量文件。

表型文件

EMMAX每次只能运行一个性状,需要的表型文件格式与tfam文件格式类似。

开始时我的表型数据格式如下:

Taxa    ts
TH001    289.1420384
TH002    278.9890168
TH003    293.6657575
TH004    285.1359311
TH005    358.0060475

转换后的格式如下:

TH001    TH001    289.1420384
TH002    TH002    278.9890168
TH003    TH003    293.6657575
TH004    TH004    285.1359311
TH005    TH005    358.0060475

GWAS运行

所有所需文件准备好后,可以开始进行GWAS分析:

emmax-intel64 -v -d 10 -t 299samples.SNP.recode.onlyChr.renameChr.maf0.05.int0.9 -o 299samples.SNP.emmax -p 2018-2019-total_sterols.overlap.average.format.txt -k 299samples.SNP.recode.onlyChr.renameChr.maf0.05.int0.9.aBN.kinf -c  emmax_cov.txt

最后生成的299samples.SNP.emmax.ps就是后续画图所需的文件。

该文件格式为(第1-4列分别为ID、回归系数、回归系数的标准差、P值):

1_124600        -4.712679096    3.657835178     0.198655268
1_124713        6.921010685     5.754786514     0.2301053679
1_124751        -0.8592932091   3.401319244     0.8007315474
1_124787        -3.745950601    3.81551729      0.3270442088
1_125023        -1.72110757     3.516697732     0.624926493

可视化

可视化需要将上述生成的ps文件转换下格式,抓换后的格式为(第1-4列分别为ID、chr、position、P值):

1_124600        1       124600  0.198655268
1_124713        1       124713  0.2301053679
1_124751        1       124751  0.8007315474
1_124787        1       124787  0.3270442088
1_125023        1       125023  0.624926493

转换过程可由单行命令完成:

perl -alne '@a=split(/_/,$F[0]);print "$F[0]\t$a[0]\t$a[1]\t$F[3]";' 299samples.SNP.emmax.ps >299samples.SNP.emmax.for_vis.ps

然后可根据如下代码进行可视化:

if(!require(CMplot)){
  install.packages("CMplot")
  library(CMplot)
}
gwas <- read.table("299samples.SNP.emmax.for_vis.ps", header = T)
head(gwas)
threshold <- 1/nrow(gwas[!is.na(gwas$pos),])
CMplot(gwas, threshold = threshold, amplify = F, memo = "", file = "tiff", plot.type=c("m","q"), ylim = c(0,12))

其中,

plot.type参数可选择"d", "c", "m", "q" or "b",

plot.type="d", SNP密度图
plot.type="c", 环形曼哈顿图
plot.type="m", 曼哈顿图
plot.type="q", QQ图
plot.type="b", 同时绘制环形曼哈顿图、曼哈顿图和QQ图
plot.type=c("m","q"), 绘制曼哈顿图和QQ图

利用rMVP进行GWAS

rMVP的运行需要准备以下文件:表型文件、基因型文件、亲缘关系矩阵、PCA矩阵。

表型文件

表型文件可以包含一个性状或多个性状,而且材料名的顺序可以不与基因型文件中的顺序一致,MVP可以自动自动调整基因型文件中材料的顺序使得其与表型文件里的一致。

一个性状:

Taxa    ts
TH001    289.1420384
TH002    278.9890168
TH003    293.6657575
TH004    285.1359311
TH005    358.0060475

多个性状:

Taxa    trait1    trait2    trait3
33-16    101.5    0.25    0
38-11    102.7    0.23    1
4226    101.2    -0.17    1
4722    105.5    NA    0
A188    108.1    0.57    1
A214N    95.13    0.87    0
A239    100.2    -0.16    1

基因型文件

rMVP支持的基因型文件格式比较多,包括plink bed、VCF、hapmap、numeric.

如果你的基因型文件是plink binary(bed)格式,则可以按照如下方式将格式转换为MVP格式:

# Full-featured function (Recommended)
MVP.Data(fileBed="plink",
         filePhe=NULL,
         fileKin=FALSE,
         filePC=FALSE,       
         #priority="speed",
         #maxLine=10000,
         out="mvp.plink"
         )

# Only convert genotypes
MVP.Data.Bfile2MVP(bfile="plink", out='mvp', maxLine=1e4, priority='speed') # the genotype data should be fully imputed before using this function

其中,

fileBed:plink binary格式的基因型文件名

fileKin:TRUE 或 FALSE,如果为TRUE,亲缘关系矩阵将自动计算

filePC:TRUE 或 FALSE,如果为TRUE,PCA矩阵将会自动计算

out:输出文件名的前缀

priority:"speed" 或 "memory", speed模式将消耗更多内存,memory模式使用内存较少,但也更慢

maxLine:如果上一个参数为memory模式,该值就是读进内存的数据行数

VCF

如果基因型文件格式为VCF格式,则可以按照如下方式将格式转换为MVP格式:

# Full-featured function (Recommended)
MVP.Data(fileVCF="myVCF.vcf",
         #filePhe="Phenotype.txt",
         fileKin=FALSE,
         filePC=FALSE,
         out="mvp.vcf"
         )

# Only convert genotypes
MVP.Data.VCF2MVP("myVCF.vcf", out='mvp') # the genotype data should be fully imputed before using this function

其中:

fileVCF, VCF格式的基因型文件
filePhe, 表型文件
fileKin, TRUE 或 FALSE,如果为TRUE,亲缘关系矩阵将自动计算
filePC, TRUE 或 FALSE,如果为TRUE,PCA矩阵将会自动计算
out, 输出文件名的前缀

Hapmap

如果基因型文件格式为Hapmap格式,则可以按照如下方式将格式转换为MVP格式:

# hapmap format example
rs#    alleles    chrom    pos    strand    assembly#    center    protLSID    assayLSID    panelLSID    QCcode    33-16    38-11    4226    4722    A188    ...    A239
rs3683945    G/A    1    3197400    +    NA    NA    NA    NA    NA    NA    AG    AG    GG    AG    GG    ...    AA
rs3707673    A/G    1    3407393    +    NA    NA    NA    NA    NA    NA    GA    GA    AA    GA    AA    ...    GG
rs6269442    G/A    1    3492195    +    NA    NA    NA    NA    NA    NA    AG    GG    GG    AG    GG    ...    AA
rs6336442    G/A    1    3580634    +    NA    NA    NA    NA    NA    NA    AG    AG    GG    AG    GG    ...    AA
rs13475699    G    1    3860406    +    NA    NA    NA    NA    NA    NA    GG    GG    GG    GG    GG    ...    GG
# Full-featured function (Recommended)
MVP.Data(fileHMP="hapmap.txt",
         filePhe="Phenotype.txt",
         sep.hmp="\t",
         sep.phe="\t",
         SNP.effect="Add",
         fileKin=FALSE,
         filePC=FALSE,
         #priority="memory",
         #maxLine=10000,
         out="mvp.hmp"
         )

# Only convert genotypes
MVP.Data.Hapmap2MVP("hapmap.txt", out='mvp') # the genotype data should be fully imputed before using this function

其中:

fileHMP, hapmap格式的基因型文件
filePhe, 表型文件
sep.phe, 表型文件分隔符
fileKin, TRUE 或 FALSE,如果为TRUE,亲缘关系矩阵将自动计算
filePC, TRUE 或 FALSE,如果为TRUE,PCA矩阵将会自动计算
priority, "speed" 或 "memory", speed模式将消耗更多内存,memory模式使用内存较少,但也更慢
maxLine, 如果上一个参数为memory模式,该值就是读进内存的数据行数

Numeric

如果基因型文件格式为Numeric格式,则可以按照如下方式将格式转换为MVP格式:

# Numeric.txt
1    1    2    1    2    …    0
1    1    0    1    0    …    2
1    2    2    1    2    …    0
1    1    2    1    2    …    0
0    0    0    0    0    …    0
# Map.txt
SNP    Chr    Pos
rs3683945    1    3197400
rs3707673    1    3407393
rs6269442    1    3492195
rs6336442    1    3580634
rs13475699    1    3860406
# Full-featured function (Recommended)
MVP.Data(fileNum="Numeric.txt",
         filePhe="Phenotype.txt",
         fileMap="Map.txt",
         sep.num="\t",
         sep.map="\t", 
         sep.phe="\t",
         fileKin=FALSE,
         filePC=FALSE,
         #priority="memory",
         #maxLine=10000,
         out="mvp.num"
         )

# Only convert genotypes
MVP.Data.Numeric2MVP("Numeric.txt", out='mvp', maxLine=1e4, priority='speed', auto_transpose=T) # the genotype data should be fully imputed before using this function

亲缘关系矩阵

如果你已经有了亲缘关系矩阵,并且格式类似下面的(n*n矩阵,n为材料数),可以直接读入:

0.3032    -0.0193    0.0094    0.0024    0.0381    ...    -0.0072
-0.0193    0.274    -0.0243    0.0032    -0.0081    ...    0.0056
0.0094    -0.0243    0.3207    -0.0071    -0.0045    ...    -0.0407
0.0024    0.0032    -0.0071    0.321    -0.008    ...    -0.0093
0.0381    -0.0081    -0.0045    -0.008    0.3498    ...    -0.0238
...    ...    ...    ...    ...    ...    ...
-0.0072    0.0056    -0.0407    -0.0093    -0.0238    ...    0.3436
# read from file
MVP.Data.Kin("mvp.kin.txt", out="mvp", maxLine=1e4, priority='memory', sep='\t')

# calculate from mvp_geno_file
MVP.Data.Kin(TRUE, mvp_prefix='mvp', out='mvp')

PCA矩阵

如果你已经有了如下格式的PCA矩阵:

0.010175524    -0.037989071    0.009588312
-0.009138673    -0.036763080    -0.006396714
-0.004723734    -0.047837625    0.021687731
0.012887843    -0.048418352    0.054298850
0.003871951    -0.038070387    0.008020508
-0.079505846    0.005818163    -0.206364549

可以直接读入:

# read from file
MVP.Data.PC("mvp.pc.txt", out='mvp', out=NULL, sep='\t')

或者根据基因型文件直接计算:

# calculate from mvp_geno_file
MVP.Data.PC(TRUE, mvp_prefix='mvp.plink', out='mvp', pcs.keep=5)

开始运行GWAS

先读入MVP生成的基因型、表型、map、亲缘关系矩阵和协变量(PCA)文件。

genotype <- attach.big.matrix("mvp.geno.desc")
phenotype <- read.table("mvp.phe",head=TRUE)
map <- read.table("mvp.geno.map" , head = TRUE)
Kinship <- attach.big.matrix("mvp.kin.desc")
Covariates_PC <- bigmemory::as.matrix(attach.big.matrix("mvp.pc.desc"))

如果是表型文件只含有一个性状,则可以使用如下代码运行:

imMVP <- MVP(
    phe=phenotype,
    geno=genotype,
    map=map,
    #K=Kinship,
    #CV.GLM=Covariates,  ##if you have additional covariates, please keep there open.
    #CV.MLM=Covariates,
    #CV.FarmCPU=Covariates,
    nPC.GLM=5,   ##if you have added PCs into covariates, please keep there closed.
    nPC.MLM=3,  ##if you don't want to add PCs as covariates, please comment out the parameters instead of setting the nPC to 0.
    nPC.FarmCPU=3,
    priority="speed",   ##for Kinship construction
    #ncpus=10,
    vc.method="BRENT",  ##only works for MLM
    maxLoop=10,
    method.bin="static",   ## "FaST-LMM", "static" (#only works for FarmCPU)
    #permutation.threshold=TRUE,
    #permutation.rep=100,
    threshold=0.05,
    method=c("GLM", "MLM", "FarmCPU")
)

如果是表型文件中包含多个性状,则就可以使用如下代码:

for(i in 2:ncol(phenotype)){
  imMVP <- MVP(
    phe=phenotype[, c(1, i)],
    geno=genotype,
    map=map,
    #K=Kinship,
    #CV.GLM=Covariates,
    #CV.MLM=Covariates,
    #CV.FarmCPU=Covariates,
    nPC.GLM=5,
    nPC.MLM=3,
    nPC.FarmCPU=3,
    priority="speed",
    #ncpus=10,
    vc.method="BRENT",
    maxLoop=10,
    method.bin="static",
    #permutation.threshold=TRUE,
    #permutation.rep=100,
    threshold=0.05,
    method=c("GLM", "MLM", "FarmCPU")
  )
  gc()
}

运行完之后会输出多个文件,包括表型分布图、SNP密度分布图、PCA图、曼哈顿圈图、常规曼哈顿图、QQ图以及候选的显著SNP。如果还需要其他格式的图片,可以参考rMVP的说明文档

利用GAPIT进行GWAS

运行GAPIT的相关代码如下,需要输入的文件包括基因型文件(hapmap或numeric格式)、表型文件:

library("gplots")
library("LDheatmap")
library("genetics")
library("ape")
library("compiler")

library("EMMREML")
library("scatterplot3d")
library("multtest")
library("snpStats")
source("https://zzlab.net/GAPIT/gapit_functions.txt")https://zzlab.net/GAPIT/gapit_functions.txt
# import data
myY=read.table(file="../2018-2019-ts.overlap.average.txt", head = TRUE)
myG=read.table(file="../299samples.SNP.onlyChr.hmp.txt", head = FALSE)
# run GAPIT with multiple models
myGAPIT <- GAPIT(
  #Y=myY[,c(1,3)],
  Y=myY, # 表型文件
  #GD=myGD,
  #GM=myGM,
  G=myG, # 基因型文件
  SNP.MAF=0.05, # Minor Allele Frequency to Filter SNPs
  model=c("GLM","MLM","MLMM","FarmCPU","Blink"),
  PCA.total=3, # 前三个主成分进行群体结构校正
  Model.selection = TRUE, 
  kinship.cluster=c("average", "complete", "ward"), # change lustering algorithms
  kinship.group=c("Mean", "Max"), # change kinship  summary  statistic  
  NJtree.group=4,     # set the number of clusting group in Njtree plot
  #QTN.position=mysimulation$QTN.position,
  Inter.Plot=TRUE,  # perform interactive plot
  PCA.3d=TRUE,      # plot 3d interactive PCA
  file.output=T
  ) 

其他详细说明可以参考GAPIT说明文档

GWAS结果筛选

GWAS结果中的曼哈顿图显示的是每个SNP在关联分析中的显著性水平;QQ图则反映的是关联分析的效果。

基因位点在y轴的高度对应该位点与表型的关联程度,关联程度越强,y值越大。受LD影响,基因组上强关联位点周围的SNP也会呈现出关联性由高到低连续变化的信号强度,从而在P值小的地方出现峰。峰值点附近的这种信号变化符合群体遗传重组模式,可能是一个可靠位点。

QQ图通过比较每个SNP期望P值与观测P值的差异来对GWAS结果进行质控。GWAS假设只有一小部分SNP与表型相关,因此大部分SNP期望P值应该与观测P值重合。

基于LD衰减距离和显著关联SNP,通常有2种方式确定候选区间。

  1. 将显著关联的SNP在N kb以内的位置确定为相关区间;

  2. N kb以内的位置相近的SNP定义为一个cluster。

其中,N为LD衰减距离。

GWAS鉴定出候选区间后整合多方面信息来精选候选基因。一般符合以下条件的基因值得进行验证和深入研究。

  1. 信号模式:关联性从高到低连续变化;

  2. 峰值区间内的基因功能注释与表型相关;

  3. 其他实验功能研究或组学数据支持GWAS结果。

参考资料