R语言rMVP包做GWAS（全基因组关联分析）分析实例

用户7010445

发布于 2024-05-18 08:47:39

4571

发布于 2024-05-18 08:47:39

本篇推文的示例数据来源于参考书《Genome-Wide Association Studies》的第十章 A Practical Guide to Using Structural Variants for Genome-Wide Association Studies。

植物里做GWAS分析通常是选择某个群体做二代基因组测序（有的已经研究比较多的物种比如水稻、玉米可以能已经发表过很多数据，），测序数据与参考基因组进行比对鉴定变异位点，然后用变异位点和表型特征去做关联分析。变异位点包括snp indel SV SRR STR TIP 等等。表型数据通常常规的植物表型，比如株高、产量、花瓣数量等。也可以是分子层面的代谢物，基因表达量等。也可以是环境数据，鉴定哪些变异位点和环境有关。

变异位点的数据通常是用vcf文件存储（当然也有其他格式）。我这篇推文介绍用vcf文件去做GWAS。

首先介绍一下vcf文件的格式

vcf文件是文本文件，我们在自己电脑上直接用记事本打开就可以查看文件里的内容。通常样本很多的话，对应的文件也会非常大。自己电脑打开还挺费劲的。

vcf 文本里的内容按照特定的模式排列

vcf简单的可以划分为三个部分

1、两个#号开头的行

2、一个#号开头的行

3、零个#号开头的行

表型数据

两列，第一列是样本名字，第二列是表型的值

如果是用rMVP这个R包来做GWAS的话表型数据的样本顺序和vcf文件的样本顺序不一致也可以，但是其他软件有的会要求样本顺序一致

rMVP 这个R包的github主页

https://github.com/xiaolei-lab/rMVP

只使用混合线性模型的代码

setwd("D:/Jupyter/practice/rMVP_GWAS/")

library(rMVP)

MVP.Data(fileVCF="smoove_filtered.vcf",
         filePhe="phenotype.txt",
         sep.phe = "\t",
         fileKin=TRUE,
         filePC=TRUE,
         out="smoove"
)

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


imMVP <- MVP(
  phe=phenotype,
  geno=genotype,
  map=map,
  K=Kinship,
  CV.MLM=Covariates_PC,
  nPC.MLM=10,
  priority="speed",   ##for Kinship construction
  ncpus=2,
  vc.method="BRENT",  ##only works for MLM
  maxLoop=10,
  threshold=0.05,
  method=c("MLM"),
  file.output=c("pmap", "pmap.signal", "plot", "log")
)