In byandell/qtl2feather: Feather Database for R/qtl2 Genotype Probabilities

knitr::opts_chunk$set(echo = TRUE, fig.width = 7, fig.height = 5)

The goal is to use feather_genoprob object with scan1 and other tools in R/qtl2.

suppressPackageStartupMessages({
  library(qtl2)
  library(qtl2feather)
})

DOex <- 
  read_cross2(
    file.path("https://raw.githubusercontent.com/rqtl",
              "qtl2data/master/DOex",
              "DOex.zip"))

Calculate genotype probabilities and convert to allele probabilities

pr <- calc_genoprob(DOex, error_prob=0.002)
apr <- genoprob_to_alleleprob(pr)

Write feather database of pr.

tmpdir <- tempdir()

fpr <- feather_genoprob(pr, "pr", tmpdir)
fapr <- feather_genoprob(apr, "apr", tmpdir)

Genome scan

scan_apr <- scan1(fapr, DOex$pheno)

find_peaks(scan_apr, DOex$pmap)

plot(scan_apr, DOex$pmap)

coefs <- scan1coef(apr, DOex$pheno)

plot(coefs, DOex$pmap)

Plot allele effects over LOD scan.

plot(coefs, DOex$pmap, scan1_output = scan_apr)

SNP probabilities

Since SNPs are in RDS, need to first download file, then read it. This file only has snpinfo for chr 2.

file <- file.path("https://raw.githubusercontent.com/rqtl",
                  "qtl2data/master/DOex", 
                  "c2_snpinfo.rds")
tmpfile <- tempfile()
download.file(file, tmpfile, quiet=TRUE)
snpinfo <- readRDS(tmpfile)
unlink(tmpfile)

Create index to non-equivalent set of SNPs.

snpinfo <- index_snps(DOex$pmap, snpinfo)

Convert to snp probabilities using feather database. Need to pick a chromosome, or first one is chosen.

snppr <- genoprob_to_snpprob(subset(fapr, chr = "2"), snpinfo)

object.size(snppr)

Kinship calculation

kinship <- calc_kinship(apr, "loco")

fkinship <- calc_kinship(fapr, "loco")

summary(c(kinship[["2"]] - fkinship[["2"]]))

Perform SNP association analysis.

scan_snppr <- scan1(snppr, DOex$pheno, fkinship["2"])

Plot results

plot(scan_snppr, snpinfo, drop_hilit = 1.5)

byandell/qtl2feather documentation built on May 13, 2019, 9:30 a.m.