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