tests/testthat/test-vcf-reader.R

In vcfppR: Rapid Manipulation of the Variant Call Format (VCF)

library(testthat)

vcffile <- system.file("extdata", "raw.gt.vcf.gz", package="vcfppR")
svfile <- system.file("extdata", "sv.vcf.gz", package="vcfppR")

test_that("vcfreader: constructor with vcfffile only", {
  br <- vcfreader$new(vcffile)
  expect_true(br$variant())
  expect_identical(br$chr(), "chr21")
  expect_identical(br$pos(), 5030082L)
  expect_identical(br$id(), "chr21:5030082:G:A")
  expect_identical(br$ref(), "G")
  expect_identical(br$alt(), "A")
})

test_that("vcfreader: constructor with vcfffile and region", {
  br <- vcfreader$new(vcffile, "chr21:5030082-")
  expect_true(br$variant())
  expect_identical(br$chr(), "chr21")
  expect_identical(br$pos(), 5030082L)
  expect_identical(br$id(), "chr21:5030082:G:A")
  expect_identical(br$ref(), "G")
  expect_identical(br$alt(), "A")
})

test_that("vcfreader: constructor with vcfffile, region and samples", {
  br <- vcfreader$new(vcffile, "chr21:5030082-", "HG00097,HG00099")
  expect_identical(br$samples(), c("HG00097", "HG00099"))
  expect_true(br$variant())
  expect_identical(br$genotypes(F), rep(0L, 4))
  expect_identical(br$chr(), "chr21")
  expect_identical(br$pos(), 5030082L)
  expect_identical(br$id(), "chr21:5030082:G:A")
  expect_identical(br$ref(), "G")
  expect_identical(br$alt(), "A")
})

test_that("vcfreader: get FORMAT item with Float type", {
  br <- vcfreader$new(vcffile, "chr21:5030347-", "HG00097,HG00099")
  expect_identical(br$samples(), c("HG00097", "HG00099"))
  expect_true(br$variant())
  expect_identical(br$pos(), 5030347L)
  expect_identical(all(is.na(br$formatFloat("AB"))), TRUE)
})

test_that("vcfreader: get FORMAT item with Integer type", {
  br <- vcfreader$new(vcffile, "chr21:5030347-", "HG00097,HG00099")
  expect_identical(br$samples(), c("HG00097", "HG00099"))
  expect_true(br$variant())
  expect_identical(br$pos(), 5030347L)
  expect_identical(br$formatInt("AD"), c(4L, 0L, 16L, 0L))
})

test_that("vcfreader: get FORMAT item with String type", {
  br <- vcfreader$new(svfile, "chr21:5114000-", "HG00096,HG00097")
  expect_identical(br$samples(), c("HG00096", "HG00097"))
  expect_true(br$variant())
  expect_identical(br$pos(), 5114000L)
  expect_identical(br$formatStr("EV"), c("RD","RD"))
})

test_that("vcfreader: get FORMAT item for the unexpected and error", {
  br <- vcfreader$new(vcffile, "chr21:5030347-", "HG00097,HG00099")
  br$variant()
  expect_identical(br$pos(), 5030347L)
  expect_error(br$formatInt("AA")) ## AA not exists
  expect_identical(br$formatInt("AD"), c(4L, 0L, 16L, 0L))
  ## use formatFloat will return unexpected values if AD in FORMAT if Integer
  expect_error(expect_identical(br$formatFloat("AD"), c(4L, 0L, 16L, 0L))) ## AA exists but of Integer type
})

test_that("vcfreader: get variants type", {
  br <- vcfreader$new(vcffile)
  i <- 0
  s <- 0
  m <- 0
  ms <- 0
  while(br$variant()) {
    if(br$isIndel()) i <- i + 1
    if(br$isSV()) s <- s + 1
    if(br$isMultiAllelics()) m <- m + 1
    if(br$isMultiAllelicSNP()) ms <- ms + 1
  }
  expect_identical(i, 1)
  expect_identical(s, 0)
  expect_identical(m, 2)
  expect_identical(ms, 2)
})

test_that("vcfreader: test variants type", {
  br <- vcfreader$new(vcffile)
  i1 <- 0
  i2 <- 0
  i3 <- 0
  i4 <- 0
  i5 <- 0
  i6 <- 0
  i7 <- 0
  i8 <- 0
  while(br$variant()) {
    if(br$hasSNP()) i1 <- i1 + 1
    if(br$hasINDEL()) i2 <- i2 + 1
    if(br$hasINS()) i3 <- i3 + 1
    if(br$hasDEL()) i4 <- i4 + 1
    if(br$hasMNP()) i5 <- i5 + 1
    if(br$hasBND()) i6 <- i6 + 1
    if(br$hasOTHER()) i7 <- i7 + 1
    if(br$hasOVERLAP()) i8 <- i8 + 1
  }
  expect_identical(i1, 14)
  expect_identical(i2, 1)
  expect_identical(i3, 0)
  expect_identical(i4, 1)
  expect_identical(i5, 0)
  expect_identical(i6, 0)
  expect_identical(i7, 0)
  expect_identical(i8, 0)
})

test_that("vcfreader: reading variant only", {
  br <- vcfreader$new(vcffile)
  br$variant()
  expect_identical(br$chr(), "chr21")
  expect_identical(br$pos(), 5030082L)
  expect_identical(br$id(), "chr21:5030082:G:A")
  expect_identical(br$ref(), "G")
  expect_identical(br$alt(), "A")
  expect_equal(br$qual(), 70.06, tolerance = 1e-6)
  expect_identical(br$filter(), "VQSRTrancheSNP99.80to100.00")
  expect_identical(br$info(), "AC=2;AF=0.000616523;AN=3244;DP=2498;FS=0;MLEAC=1;MLEAF=0.0003083;MQ=17.07;MQ0=0;QD=17.52;SOR=3.258;VQSLOD=-32.6;culprit=MQ;AN_EUR=658;AN_EAS=548;AN_AMR=520;AN_SAS=624;AN_AFR=894;AF_EUR=0;AF_EAS=0;AF_AMR=0.00384615;AF_SAS=0;AF_AFR=0;AC_EUR=0;AC_EAS=0;AC_AMR=2;AC_SAS=0;AC_AFR=0;AC_Het_EUR=0;AC_Het_EAS=0;AC_Het_AMR=0;AC_Het_SAS=0;AC_Het_AFR=0;AC_Het=0;AC_Hom_EUR=0;AC_Hom_EAS=0;AC_Hom_AMR=2;AC_Hom_SAS=0;AC_Hom_AFR=0;AC_Hom=2;HWE_EUR=1;ExcHet_EUR=1;HWE_EAS=1;ExcHet_EAS=1;HWE_AMR=0.00192678;ExcHet_AMR=1;HWE_SAS=1;ExcHet_SAS=1;HWE_AFR=1;ExcHet_AFR=1;HWE=0.000308356;ExcHet=1;ME=0;AN_EUR_unrel=508;AN_EAS_unrel=448;AN_AMR_unrel=330;AN_SAS_unrel=484;AN_AFR_unrel=666;AF_EUR_unrel=0;AF_EAS_unrel=0;AF_AMR_unrel=0;AF_SAS_unrel=0;AF_AFR_unrel=0;AC_EUR_unrel=0;AC_EAS_unrel=0;AC_AMR_unrel=0;AC_SAS_unrel=0;AC_AFR_unrel=0;AC_Het_EUR_unrel=0;AC_Het_EAS_unrel=0;AC_Het_AMR_unrel=0;AC_Het_SAS_unrel=0;AC_Het_AFR_unrel=0;AC_Hom_EUR_unrel=0;AC_Hom_EAS_unrel=0;AC_Hom_AMR_unrel=0;AC_Hom_SAS_unrel=0;AC_Hom_AFR_unrel=0")
  br$rmInfoTag("AC")
  print(br$infoInt("AC")) ## TODO: AC value still exists
  print(br$info())
  af <- br$infoFloat("AF")
  ## expect_equal(0.00061652297, 0.00061652300, tolerance = 1e-6)
  expect_equal(af, 0.000616523, tolerance = 1e-6)
  dp <- br$infoInt("DP")
  expect_identical(dp, 2498L)
  gt <- br$genotypes(TRUE)
  expect_identical(length(gt), 3202L)
  gt <- br$genotypes(FALSE)
  expect_identical(length(gt), 6404L)
  pl <- br$formatInt("PL")
  expect_identical(length(pl), 9606L)
})

test_that("vcfreader: read and output variant", {
  br <- vcfreader$new(vcffile)
  br$variant()
  br$setCHR("chr22")
  expect_identical(br$chr(), "chr22")
  br$setPOS(9999L)
  expect_identical(br$pos(), 9999L)
  br$setID("rs122334")
  expect_identical(br$id(), "rs122334")
  expect_identical(br$ref(), "G")
  expect_identical(br$alt(), "A")
  br$setRefAlt("A,GA")
  expect_identical(br$ref(), "A")
  expect_identical(br$alt(), "GA")
  expect_identical(br$infoInt("AC"), 2L)
  br$setInfoInt("AC", 3L)
  expect_identical(br$infoInt("AC"), 3L)
  br$setInfoFloat("AF", 0.3)
  expect_equal(br$infoFloat("AF"), 0.3, tolerance  = 1e-6)
  br$setInfoStr("VariantType", "indel")
  expect_identical(br$infoStr("VariantType" ), "indel")
  ## output current variant to another vcf
  s1 <- br$string()
  outvcf <- file.path(tempdir(), "test.vcf.gz")
  file.create(outvcf)
  br$output(outvcf)
  br$write()
  br$close()
  br <- vcfreader$new(outvcf)
  br$variant()
  s2 <- br$string()
  expect_identical(s1, s2)
})


test_that("vcfreader: remove tag from FORMAT", {
  br <- vcfreader$new(vcffile)
  br$variant()  ## first variant
  s <- unlist(strsplit(br$line(), "\t"))
  expect_identical(s[9], "GT:AD:DP:GQ:PL")
  br$variant()  ## second variant
  br$rmFormatTag("AD")
  br$rmFormatTag("AB")
  br$rmFormatTag("PGT")
  br$rmFormatTag("PID")
  s <- unlist(strsplit(br$line(), "\t"))
  expect_identical(s[9], "GT:DP:GQ:PL")
  expect_identical(s[10], "1/1:2:6:64,6,0")
  ## AD was removed, so we get integer(0)
  ad <- br$formatInt("AD")
  expect_identical(length(ad),0L)
  ## output current variant to another vcf
  outvcf <- file.path(tempdir(), "test.vcf.gz")
  file.create(outvcf)
  br$output(outvcf)
  br$write()
  br$close()
  ## check the output file
  br <- vcfreader$new(outvcf)
  br$variant()
  s <- unlist(strsplit(br$line(), "\t"))
  expect_identical(s[9], "GT:DP:GQ:PL")
  expect_identical(s[10], "1/1:2:6:64,6,0")
  ## AD doesn't exist. so we get error, is this a bad design?
  expect_error(br$formatInt("AD"))
})


test_that("can set genotypes for single sample", {
  
  br <- vcfreader$new(svfile, "", "HG00096")
  br$variant()
  br$genotypes(F)
  br$setGenotypes(c(1L,1L))
  outfile <- paste0(tempfile(), ".vcf.gz")
  br$output(outfile)
  br$write()
  br$close()

  vcf <- vcftable(outfile)
  expect_true(vcf$gt==2)
  
})

test_that("can change samples name and set genotypes for single sample", {
  br <- vcfreader$new(svfile, "", "HG00096")
  br$variant()
  expect_identical(br$infoStr("SVTYPE"), "DUP")
  expect_identical(br$genotypes(F), c(0L, 0L))
  br$setGenotypes(c(1L,1L))
  expect_identical(br$genotypes(F), c(1L, 1L))
  outfile <- paste0(tempfile(), ".vcf.gz")
  br$output(outfile)
  br$updateSamples("ZZZZZ")
  br$write()
  br$close()
  vcf <- vcftable(outfile)
  expect_true(vcf$gt==2)
  expect_true(vcf$samples=="ZZZZZ")
})

test_that("can getStatus of a region", {
  
  br <- vcfreader$new(svfile)
  expect_identical(br$getStatus("chr22"),-2L)
  expect_identical(br$getStatus("XXX-1-9999"),-2L)
  expect_identical(br$getStatus("chr21:5029001-5029100"), 0L)
  expect_identical(br$getStatus("chr21"), 1L)
  
})

test_that("can set region back and forth", {
  
  br <- vcfreader$new(svfile)
  br$setRegion("chr21:5227411")
  br$variant()
  expect_identical(br$infoStr("SVTYPE"), "INS")
  ## set region back
  br$setRegion("chr21:5114000")
  br$variant()
  expect_identical(br$infoStr("SVTYPE"), "DUP")
  
})