research: molgenis_research_demoanalysis_Smith_Kruglyak_2008.R

File molgenis_research_demoanalysis_Smith_Kruglyak_2008.R, 1.1 KB (added by jvelde, 4 years ago)

Yeast demo analysis for MOLGENIS Research

Line 
1# first time only: install R/qtl
2install.packages("qtl")
3
4# load library and connect to molgenis
5library(qtl)
6source("https://molgenis83.gcc.rug.nl/molgenis.R")
7molgenis.login("admin", "mrdemo")
8
9# download and setup data
10yeastdata <- molgenis.get("base_yeast_geno_pheno_map")
11tmpfile <- tempfile("tmp_yeastdata")
12yeastdata <- yeastdata[with(yeastdata, order(chr,na.last = FALSE)), ]
13write.table(yeastdata, tmpfile, sep = ",", col.names = FALSE, row.names = FALSE, quote = FALSE, na = "")
14
15# read data and plot genotypes
16yeast <- read.cross(format="csvr", file=tmpfile, genotypes=c("1","2"))
17yeast <- calc.genoprob(yeast)
18yeast <- fill.geno(yeast)
19geno.image(yeast)
20
21# perform qtl scan on the first phenotype
22result_hk <- scanone(yeast, pheno.col=1, method="hk")
23perm_hk <- scanone(yeast, pheno.col=1, method="hk", n.perm=1000)
24plot(result_hk)
25abline(h=summary(perm_hk)[[1]],col="green",lty=2,lwd=2)
26abline(h=summary(perm_hk)[[2]],col="orange",lty=2,lwd=2)
27
28# get top marker and plot gene expression
29topmarker <- find.marker(yeast, chr=max(result_hk)$chr, pos=max(result_hk)$pos)
30plotPXG(yeast, marker=topmarker)
31
32# logout
33molgenis.logout()