-
Notifications
You must be signed in to change notification settings - Fork 7
/
rMVP GWAS script.R
86 lines (72 loc) · 2.63 KB
/
rMVP GWAS script.R
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
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
################################################################################
# rMVP GWAS
################################################################################
# devtools::install_github("xiaolei-lab/rMVP")
install.packages("rMVP")
library(rMVP)
################################################################################
# Step1: formatting of data
################################################################################
rm(list = ls()) #clear R environment
MVP.Data(fileHMP="mdp_genotype.hmp.txt", #Genotypic data in HapMap format
filePhe="mdp_traits.txt",
sep.hmp="\t",
sep.phe="\t",
SNP.effect="Add",
fileKin=T,
out="mvp.hmp"
)
################################################################################
# Step2: import formatted data (MVP.Data) from working directory
################################################################################
genotypic_dat <- attach.big.matrix("mvp.hmp.geno.desc")
phenotype_dat <- read.table("mvp.hmp.phe",head=TRUE)
map_info <- read.table("mvp.hmp.geno.map" , head = TRUE)
popstr <- read.table("mdp_population_structure.txt", header = T, skip = 1)
Kinship <- MVP.K.VanRaden(genotypic_dat, verbose = T)
################################################################################
##Step2: GWAS run
################################################################################
################################################################################
# PCA as covariate
################################################################################
GWAS_PCA_mvp <- MVP(
phe = phenotype_dat,
geno = genotypic_dat,
map = map_info,
method = c("GLM", "MLM", "FarmCPU"),
nPC.GLM = 5,
nPC.MLM = 5,
nPC.FarmCPU = 5
)
################################################################################
# Kinship as covariate
################################################################################
dir.create("PC_K")
setwd("./PC_K/")
getwd()
GWAS_Kin_mvp <- MVP(
phe = phenotype_dat,
geno = genotypic_dat,
map = map_info,
method = c("GLM", "MLM", "FarmCPU"),
nPC.GLM = 0,
nPC.MLM = 0,
nPC.FarmCPU = 0,
K = Kinship
)
################################################################################
# PCA and Kinship as covariates
################################################################################
dir.create("../PcaKinship" )
setwd("../PcaKinship")
GWAS_PCA_Kin_mvp <- MVP(
phe = phenotype_dat,
geno = genotypic_dat,
map = map_info,
method = c("GLM", "MLM", "FarmCPU"),
nPC.GLM = 5,
nPC.MLM = 5,
nPC.FarmCPU = 5,
K = Kinship
)