图就用的rMVP画就好了，挺好看的
注意先根据群体质控一下基因型文件
注意先检查文件格式，检查基因型和表型文件物种名是否一致

基因型和表型文件转rMVP格式

1
2
3
4
5
6
7
8
9
10
11
## 切换 Kernel 为 R
library('rMVP')
### plink 和 phenotype 转 MVP
MVP.Data(fileBed="./data/filteredSNP-2/filter_v2",
         filePhe="./data/Phenotypes/phenotype240313.txt",
         fileKin=TRUE,
         filePC =TRUE,       
         #priority="speed",
         #maxLine=10000,
         out    ="./GWAS/manhattan_plot/results/mvp.plink"
)

分析

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
### 读 Pheno 文件
library(bigmemory)
library(rMVP)

setwd('./GWAS/04.manhattan_plot/results/')

phenotype <- read.table("./GWAS/04.manhattan_plot/results/mvp.plink.phe", header = TRUE)
names(phenotype) <- c("taxa", "root_length:compacted")

### 读 MVP.Geno 文件
genotype<- attach.big.matrix("mvp.plink.geno.desc")
map <- read.table("mvp.plink.geno.map", header = TRUE)

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

### 执行rMVP
imMVP <- MVP(
      phe =phenotype,
      geno=genotype,
      map =map,
      K   =Kinship,
      CV.GLM=Covariates_PC,  ##if you have additional covariates, please keep there open.
      CV.MLM=Covariates_PC,
      CV.FarmCPU=Covariates_PC,
      nPC.GLM=3,   ##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"),
      # file.output=c("pmap", "pmap.signal", "plot", "log")
      file.output=c("pmap", "log")
)

计算阈值

其实上一步有个阈值，非常严格的阈值。
这里用GEC计算阈值。

1
2
3
4
5
!java -jar -Xmx8g ./bin/gec/gec.jar \
--effect-number \
--plink-binary ./data/filteredSNP-2/filter_v2 \
--genome \
--out ./GWAS/manhattan_plot/results/gec/gec.result

Manhatton可视化

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
### 读 rMVP 计算的 GWAS 结果
data <- read.csv('log2CdU.FarmCPU.csv')
MVP.Report(data[, -4:-7],
	plot.type    = "m",
	col          = c("dodgerblue4","deepskyblue"),
	LOG10        = TRUE,
	ylim         = NULL,
	threshold    = c(1.01E-7,2.03E-6), # GEC阈值
	threshold.lty= c(1,2),
	threshold.lwd= c(1,1),
	threshold.col= c("black","grey"),
	amplify      = TRUE,
	chr.den.col  = c("darkgreen", "yellow", "red"),
	signal.col   = c("red","green"),
	signal.cex   = c(1,1),
	signal.pch   = c(19,19),
	file.type    = "jpg",
	memo         = "",
	dpi          = 300
	)

Q-Q Plot 可视化

1
2
3
4
5
6
7
8
9
10
data <- read.csv('UminersC.FarmCPU.csv', header=TRUE)
CMplot(data[,c(1:3,8)],
	   plot.type    = "q",
	   conf.int.col = 'grey',
	   box          = TRUE,
	   file         = "jpg",
	   dpi          = 300,
	   file.output  = TRUE,
	   verbose      = TRUE
	   )