mGWAS/htu_gwas_vesto_v3.R at main · mabro7766/mGWAS · GitHub

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
#R version 3.1.0 (2014-04-10) -- "Spring Dance"
#Copyright (C) 2014 The R Foundation for Statistical Computing
#Platform: x86_64-w64-mingw32/x64 (64-bit)


### how to use GWAS R script.

t<-as.integer(commandArgs(trailingOnly=TRUE))

cat(t)

### set workingdir: dir where phenotype data is

setwd("/scratch/ngs-bioenergy/25okt2017")

### load library(qqman)
library(qqman)

### step 1:  you need to source the following scripts
#####################################################

source('/scratch/ngs-bioenergy/25okt2017/emma.r')
## this the emma package from Hyun Min Kang
## http://mouse.cs.ucla.edu/emma/install.html

source('/scratch/ngs-bioenergy/25okt2017/amm_gwas_vesto_v3.R')
## this runs the GWAS

# source('/scratch/ngs-bioenergy/25okt2017/plots_gwas.r')
# use library qqman instead

## step 2: load the genotype and kinship data
##############################################################

load('/scratch/ngs-bioenergy/25okt2017/snps2016.RData')

load('/scratch/ngs-bioenergy/25okt2017/kinship.RData')


## step 3: load your phenotype data
##############################################################

ph <- read.delim("tRg_Slope_Sign.txt")

## step 4: run the GWAS analysis
##############################################################

# function call: amm_gwas_vesto_v3<-function(Y,X,K,p=0.001,n=2,run=T,calculate.effect.size=FALSE,
#                                            include.lm=FALSE,use.SNP_INFO=FALSE,report=T)
# p proportion minor allele freq
# Y phenotype
# X genotype
# K kinship
# n trait column (1st column is ecotype_id)

amm_gwas_vesto_v3(ph,snps2016,kinship,p=0.05,n=t, report=F)

# Sys.time()-a
## the function will output a data.frame with n=number of SNPS rows and 9 columns

## test data (165 ecotypes, 210K SNPs) run 1.4 mins on vesto's dell latitude , 64-bit version.

outputFilename=paste("output_trait",t-1,".txt",sep="")
write.table(output.sort,file=outputFilename,quote=FALSE,col.names=TRUE, row.names=FALSE,sep="\t",na=".")

gwasResults = output.sort[,c(1,2,3,8)]
colnames(gwasResults) = c("SNP","CHR","BP","P")
jpeg(paste("manhattan_trait",t-1,".jpg",sep=""), width=6, height=5, units="in", res=300, quality=100)
manhattan(gwasResults)
dev.off()