tests/testthat/old_test_pcrelate2.R

In UW-GAC/GENESIS: GENetic EStimation and Inference in Structured samples (GENESIS): Statistical methods for analyzing genetic data from samples with population structure and/or relatedness

context("pcrelate2 tests")
library(SeqVarTools)

test_that("pcrelate2", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    myrel <- suppressWarnings(pcrelate(svd, pcMat = mypcs, correct=FALSE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- SeqVarBlockIterator(svd, verbose=FALSE)
    myrel2 <- pcrelate(iterator, pcs = mypcs, verbose=FALSE)
    cols <- c("ID1", "ID2", "kin", "k0", "k2", "nsnp")
    expect_equal(myrel2$kinBtwn[,cols], kin[,cols])
    seqClose(svd)
})

test_that("pcrelate2 - variant blocks", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    myrel <- suppressWarnings(pcrelate(svd, pcMat = mypcs, correct=FALSE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- SeqVarBlockIterator(svd, variantBlock=500, verbose=FALSE)
    myrel2 <- pcrelate(iterator, pcs = mypcs, verbose=FALSE)
    cols <- c("ID1", "ID2", "kin", "k0", "k2", "nsnp")
    expect_equal(myrel2$kinBtwn[,cols], kin[,cols])
    seqClose(svd)
})

test_that("pcrelate2 - 2 sample blocks", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    myrel <- suppressWarnings(pcrelate(svd, pcMat = mypcs, correct=FALSE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- SeqVarBlockIterator(svd, verbose=FALSE)
    myrel2 <- pcrelate(iterator, pcs = mypcs, sample.block.size=50, verbose=FALSE)
    cols <- c("ID1", "ID2", "kin", "k0", "k2", "nsnp")
    expect_equal(myrel2$kinBtwn[,cols], kin[,cols])
    seqClose(svd)
})

test_that("pcrelate2 - >2 sample blocks", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    myrel <- suppressWarnings(pcrelate(svd, pcMat = mypcs, correct=FALSE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- SeqVarBlockIterator(svd, verbose=FALSE)
    myrel2 <- pcrelate(iterator, pcs = mypcs, sample.block.size=20, verbose=FALSE)
    cols <- c("ID1", "ID2", "kin", "k0", "k2", "nsnp")
    expect_equal(myrel2$kinBtwn[,cols], kin[,cols])
    seqClose(svd)
})

test_that("pcrelate2 - sample include", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    samp.incl <- sample(seqGetData(svd, "sample.id"), 50)
    myrel <- suppressWarnings(pcrelate(svd, pcMat = mypcs, scan.include=samp.incl, correct=FALSE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- SeqVarBlockIterator(svd, verbose=FALSE)
    myrel2 <- pcrelate(iterator, pcs = mypcs, sample.include=samp.incl, verbose=FALSE)
    cols <- c("ID1", "ID2", "kin", "k0", "k2", "nsnp")
    expect_equal(myrel2$kinBtwn[,cols], kin[,cols])
    seqClose(svd)
})

test_that("pcrelate2 - sample filter", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    seqSetFilter(svd, sample.sel=1:20, verbose=FALSE)
    myrel <- suppressWarnings(pcrelate(svd, pcMat = mypcs, correct=FALSE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- SeqVarBlockIterator(svd, verbose=FALSE)
    myrel2 <- pcrelate(iterator, pcs = mypcs, verbose=FALSE)
    cols <- c("ID1", "ID2", "kin", "k0", "k2", "nsnp")
    expect_equal(myrel2$kinBtwn[,cols], kin[,cols])
    seqClose(svd)
})

## this is slow
test_that("pcrelate2 - GenotypeData", {
    gd <- .testHMData()
    mypcs <- .testHMPCs(gd)
    myrel <- suppressWarnings(pcrelate(gd, pcMat = mypcs, correct=FALSE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    kin$ID1 <- as.character(kin$ID1)
    kin$ID2 <- as.character(kin$ID2)
    iterator <- GWASTools::GenotypeBlockIterator(gd)
    myrel2 <- pcrelate(iterator, pcs = mypcs, verbose=FALSE)
    tmp <- dplyr::left_join(kin, myrel2$kinBtwn, by=c("ID1", "ID2"))
    cols <- c("kin", "k0", "k2", "nsnp")
    expect_equivalent(tmp[,paste0(cols, ".x")], tmp[,paste0(cols, ".y")])
    GWASTools::close(gd)
})

test_that("pcrelate2 - GenotypeData - variant blocks", {
    gd <- .testGenoData()
    mypcs <- .testGenoPCs(gd)
    myrel <- suppressWarnings(pcrelate(gd, pcMat = mypcs, correct=FALSE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- GWASTools::GenotypeBlockIterator(gd, snpBlock=1000)
    myrel2 <- pcrelate(iterator, pcs = mypcs, verbose=FALSE)
    expect_equal(myrel2$kinBtwn$ID1, as.character(kin$ID1))
    expect_equal(myrel2$kinBtwn$ID2, as.character(kin$ID2))
    cols <- c("kin", "k0", "k2", "nsnp")
    expect_equivalent(myrel2$kinBtwn[,cols], kin[,cols])
    GWASTools::close(gd)
})

test_that("pcrelate2 - GenotypeData - sample blocks", {
    gd <- .testGenoData()
    mypcs <- .testGenoPCs(gd)
    myrel <- suppressWarnings(pcrelate(gd, pcMat = mypcs, correct=FALSE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- GWASTools::GenotypeBlockIterator(gd)
    myrel2 <- pcrelate(iterator, pcs = mypcs, sample.block.size=50, verbose=FALSE)
    expect_equal(myrel2$kinBtwn$ID1, as.character(kin$ID1))
    expect_equal(myrel2$kinBtwn$ID2, as.character(kin$ID2))
    cols <- c("kin", "k0", "k2", "nsnp")
    expect_equivalent(myrel2$kinBtwn[,cols], kin[,cols])
    GWASTools::close(gd)
})

test_that("pcrelate2 - GenotypeData - sample include", {
    gd <- .testGenoData()
    mypcs <- .testGenoPCs(gd)
    samp.incl <- sample(GWASTools::getScanID(gd), 50)
    myrel <- suppressWarnings(pcrelate(gd, pcMat = mypcs, scan.include=samp.incl, correct=FALSE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- GWASTools::GenotypeBlockIterator(gd)
    myrel2 <- pcrelate(iterator, pcs = mypcs, sample.include=samp.incl, verbose=FALSE)
    expect_equal(myrel2$kinBtwn$ID1, as.character(kin$ID1))
    expect_equal(myrel2$kinBtwn$ID2, as.character(kin$ID2))
    cols <- c("kin", "k0", "k2", "nsnp")
    expect_equivalent(myrel2$kinBtwn[,cols], kin[,cols])
    GWASTools::close(gd)
})

test_that("pcrelate2 - small sample correction", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    myrel <- suppressWarnings(pcrelate(svd, pcMat = mypcs, correct=TRUE, verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- SeqVarBlockIterator(svd, verbose=FALSE)
    myrel2 <- pcrelate(iterator, pcs = mypcs, small.samp.correct = TRUE, verbose=FALSE)
    cols <- c("ID1", "ID2", "kin", "k0", "k2", "nsnp")
    expect_equal(myrel2$kinBtwn[,cols], kin[,cols])
    seqClose(svd)
})

test_that("pcrelate2 - scale=variant", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    myrel <- suppressWarnings(pcrelate(svd, pcMat = mypcs, correct=FALSE, scale="variant", verbose=FALSE))
    kin <- suppressWarnings(pcrelateReadKinship(myrel))
    iterator <- SeqVarBlockIterator(svd, verbose=FALSE)
    myrel2 <- pcrelate(iterator, pcs = mypcs, scale="variant", verbose=FALSE)
    cols <- c("ID1", "ID2", "kin", "k0", "k2", "nsnp")
    expect_equal(myrel2$kinBtwn[,cols], kin[,cols])
    seqClose(svd)
})

test_that("pcrelate2 - scale=none", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    myrel <- suppressWarnings(pcrelate(svd, pcMat = mypcs, correct=FALSE, scale="none", ibd.probs=FALSE, verbose=FALSE))
    grm <- suppressWarnings(pcrelateMakeGRM(myrel))
    iterator <- SeqVarBlockIterator(svd, verbose=FALSE)
    myrel2 <- pcrelate(iterator, pcs = mypcs, scale="none", ibd.probs=FALSE, verbose=FALSE)
    grm2 <- pcrelateToMatrix(myrel2, verbose=FALSE)
    expect_equivalent(grm, as.matrix(grm2))
    seqClose(svd)
})

test_that("pcrelate2 - method=truncate", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    iterator <- SeqVarBlockIterator(svd, verbose=FALSE)
    myrel.f <- pcrelate(iterator, pcs = mypcs, maf.bound.method="filter", verbose=FALSE)
    myrel.t <- pcrelate(iterator, pcs = mypcs, maf.bound.method="truncate", verbose=FALSE)
    expect_true(all(myrel.t$nsnp > myrel.f$nsnp))
    expect_equal(myrel.t$kin, myrel.f$kin, tolerance=0.01)
    seqClose(svd)
})

test_that("pcrelate2 - make GRM", {
    svd <- .testData()
    mypcs <- .testPCs(svd)
    myrel <- suppressWarnings(pcrelate(svd, pcMat = mypcs, correct=FALSE, verbose=FALSE))
    grm <- suppressWarnings(pcrelateMakeGRM(myrel))
    iterator <- SeqVarBlockIterator(svd, verbose=FALSE)
    myrel2 <- pcrelate(iterator, pcs = mypcs, verbose=FALSE)
    grm2 <- pcrelateToMatrix(myrel2, verbose=FALSE)
    expect_equivalent(grm, as.matrix(grm2))
    seqClose(svd)
})