rMVP 速查表 v1.4.6 · GWAS

# 稳定版（CRAN）
install.packages("rMVP")

# 开发版（GitHub；国内可用 Gitee 镜像）
devtools::install_github("xiaolei-lab/rMVP")

library(rMVP)
?rMVP           # 全部参数
?MVP.Report     # 绘图参数

# 表型：ID + 每列一个性状，缺失写 NA，带表头
Taxa   trait1  trait2  trait3
33-16  101.5   0.25    0
4722   105.5   NA      0

# Map（Numeric 格式需要）：带表头三列
SNP        Chr  Pos
rs3683945  1    3197400

# Numeric：m×n 或 n×m 均可，值为 0/1/2

genotype  <- attach.big.matrix("mvp.geno.desc")
phenotype <- read.table("mvp.phe", head=TRUE)
map       <- read.table("mvp.geno.map", head=TRUE)

imMVP <- MVP(
  phe = phenotype,          # 两列，可含 NA
  geno = genotype,
  map = map,
  #K = Kinship,             # 有预算好的 K 才打开
  #CV.GLM = Covariates,     # 有协变量时打开全部 CV.*
  #CV.MLM = Covariates,
  #CV.FarmCPU = Covariates,
  nPC.GLM = 5,              # 不加 PC：注释掉，别写 0
  nPC.MLM = 3,
  nPC.FarmCPU = 3,
  maxLine = 10000,          # 调小可降内存
  #ncpus = 10,              # 默认用全部 CPU
  vc.method = "BRENT",      # 仅 MLM
  method.bin = "static",    # 仅 FarmCPU
  threshold = 0.05,
  method = c("GLM", "MLM", "FarmCPU"),
  file.output = c("pmap", "pmap.signal", "plot", "log")
)

# 亲缘矩阵（n×n）
MVP.Data.Kin("mvp.kin.txt", out="mvp", maxLine=1e4, sep="\t")  # 读文件
MVP.Data.Kin(TRUE, mvp_prefix="mvp.vcf", out="mvp")              # 计算

# 主成分（n×nPC）
MVP.Data.PC("mvp.pc.txt", out="mvp", sep="\t")                   # 读文件
MVP.Data.PC(TRUE, mvp_prefix="mvp.vcf", pcs.keep=5)             # 计算

# 加载
Kinship       <- attach.big.matrix("mvp.kin.desc")
Covariates_PC <- bigmemory::as.matrix(attach.big.matrix("mvp.pc.desc"))

# 品种、性别作因子，体重作数值协变量
Covariates <- model.matrix.lm(
  ~ as.factor(breed) + as.factor(sex) + as.numeric(weight),
  data = yourdata, na.action = "na.pass")   # 允许 NA

# 把 PC 并入协变量（此时 MVP() 中不要设 nPC.*）
Covariates <- cbind(Covariates, Covariates_PC)

# 双阈值 + 显著点放大 + 密度条
MVP.Report(pig60K, plot.type="m", LOG10=TRUE,
  threshold=c(1e-6,1e-4), threshold.lty=c(1,2),
  threshold.col=c("black","grey"), amplify=TRUE,
  signal.col=c("red","green"),
  chr.den.col=c("darkgreen","yellow","red"), bin.size=1e6,
  file.type="jpg", dpi=300)

# 多性状 / 多方法
MVP.Report(pig60K, plot.type="m", multracks=TRUE,
  threshold=c(1e-6,1e-4), amplify=TRUE, bin.size=1e6,
  chr.den.col=c("darkgreen","yellow","red"), file.type="jpg", dpi=300)

MVP.Report(pig60K, plot.type="c", r=0.4,
  col=c("grey30","grey60"),
  chr.labels=paste("Chr", c(1:18,"X"), sep=""),
  threshold=c(1e-6,1e-4), cir.chr.h=1.5, amplify=TRUE,
  signal.line=1, signal.col=c("red","green"),
  chr.den.col=c("darkgreen","yellow","red"), bin.size=1e6,
  outward=FALSE, file.type="jpg", dpi=300)

# 单个 QQ
MVP.Report(pig60K, plot.type="q", conf.int.col=NULL, box=TRUE,
  file.type="jpg", dpi=300)

# 多个 QQ
MVP.Report(imMVP, plot.type="q",
  col=c("dodgerblue1","olivedrab3","darkgoldenrod1"),
  conf.int.col="grey", box=FALSE, multracks=TRUE,
  file.type="jpg", dpi=300)

# SNP 密度
MVP.Report(pig60K[, 1:3], plot.type="d",
  col=c("darkgreen","yellow","red"), file.type="jpg", dpi=300)

# 画 SNP 效应（非 p 值）
MVP.Report(cattle50K, plot.type="m", band=0, LOG10=FALSE,
  ylab="Abs(SNP effect)", threshold=0.015, threshold.col="red",
  amplify=TRUE, signal.col=NULL,
  col=c("dodgerblue4","deepskyblue"), chr.den.col=NULL,
  file.type="jpg", dpi=300)

# 表型分布直方图
MVP.Hist(phe=phenotype, file.type="jpg", breakNum=18, dpi=300)

# PCA（前 3 个 PC；class 可传品种/地点，pch 设点形）
pca <- attach.big.matrix("mvp.pc.desc")[, 1:3]
MVP.PCAplot(PCA=pca, Ncluster=3, class=NULL,
  col=c("red","green","yellow"), file.type="jpg")

names(imMVP)              # map, glm.results, mlm.results, farmcpu.results
head(imMVP$mlm.results)   # Effect, SE, P

# 拼出完整表并按 p 值排序
res <- cbind(imMVP$map, imMVP$mlm.results)
res <- res[order(res[, ncol(res)]), ]
head(res, 20)

# 用返回值重画
MVP.Report(imMVP, plot.type="m", threshold=0.05/nrow(map))

for (i in 2:ncol(phenotype)) {
  imMVP <- MVP(phe = phenotype[, c(1, i)],
               geno = genotype, map = map,
               nPC.GLM = 5, nPC.MLM = 3, nPC.FarmCPU = 3,
               maxLine = 10000, vc.method = "BRENT",
               method.bin = "static", threshold = 0.05,
               method = c("GLM", "MLM", "FarmCPU"),
               file.output = c("pmap", "pmap.signal", "plot", "log"))
  gc()   # 每轮释放内存
}