Sunyoung Shin
| Snapshot Date: 2020-10-16 14:40:19 -0400 (Fri, 16 Oct 2020) |
| URL: https://git.bioconductor.org/packages/atSNP |
| Branch: RELEASE_3_11 |
| Last Commit: 1a804fa |
| Last Changed Date: 2020-04-27 15:23:26 -0400 (Mon, 27 Apr 2020) |
| Back to Multiple platform build/check report for BioC 3.11 |
|
This page was generated on 2020-10-17 11:58:04 -0400 (Sat, 17 Oct 2020).
| TO THE DEVELOPERS/MAINTAINERS OF THE atSNP PACKAGE: Please make sure to use the following settings in order to reproduce any error or warning you see on this page. |
| Package 97/1905 | Hostname | OS / Arch | INSTALL | BUILD | CHECK | BUILD BIN | ||||||
| atSNP 1.4.0 Sunyoung Shin
| malbec2 | Linux (Ubuntu 18.04.4 LTS) / x86_64 | OK | ERROR | skipped | |||||||
| tokay2 | Windows Server 2012 R2 Standard / x64 | OK | OK | OK | OK | | ||||||
| machv2 | macOS 10.14.6 Mojave / x86_64 | OK | OK | [ OK ] | OK | |
| Package: atSNP |
| Version: 1.4.0 |
| Command: /Library/Frameworks/R.framework/Versions/Current/Resources/bin/R CMD check --install=check:atSNP.install-out.txt --library=/Library/Frameworks/R.framework/Versions/Current/Resources/library --no-vignettes --timings atSNP_1.4.0.tar.gz |
| StartedAt: 2020-10-16 22:47:11 -0400 (Fri, 16 Oct 2020) |
| EndedAt: 2020-10-16 22:52:54 -0400 (Fri, 16 Oct 2020) |
| EllapsedTime: 342.1 seconds |
| RetCode: 0 |
| Status: OK |
| CheckDir: atSNP.Rcheck |
| Warnings: 0 |
############################################################################## ############################################################################## ### ### Running command: ### ### /Library/Frameworks/R.framework/Versions/Current/Resources/bin/R CMD check --install=check:atSNP.install-out.txt --library=/Library/Frameworks/R.framework/Versions/Current/Resources/library --no-vignettes --timings atSNP_1.4.0.tar.gz ### ############################################################################## ############################################################################## * using log directory ‘/Users/biocbuild/bbs-3.11-bioc/meat/atSNP.Rcheck’ * using R version 4.0.3 (2020-10-10) * using platform: x86_64-apple-darwin17.0 (64-bit) * using session charset: UTF-8 * using option ‘--no-vignettes’ * checking for file ‘atSNP/DESCRIPTION’ ... OK * checking extension type ... Package * this is package ‘atSNP’ version ‘1.4.0’ * checking package namespace information ... OK * checking package dependencies ... OK * checking if this is a source package ... OK * checking if there is a namespace ... OK * checking for hidden files and directories ... OK * checking for portable file names ... OK * checking for sufficient/correct file permissions ... OK * checking whether package ‘atSNP’ can be installed ... OK * checking installed package size ... OK * checking package directory ... OK * checking ‘build’ directory ... OK * checking DESCRIPTION meta-information ... OK * checking top-level files ... OK * checking for left-over files ... OK * checking index information ... OK * checking package subdirectories ... OK * checking R files for non-ASCII characters ... OK * checking R files for syntax errors ... OK * checking whether the package can be loaded ... OK * checking whether the package can be loaded with stated dependencies ... OK * checking whether the package can be unloaded cleanly ... OK * checking whether the namespace can be loaded with stated dependencies ... OK * checking whether the namespace can be unloaded cleanly ... OK * checking dependencies in R code ... NOTE Namespaces in Imports field not imported from: ‘graphics’ ‘testit’ All declared Imports should be used. * checking S3 generic/method consistency ... OK * checking replacement functions ... OK * checking foreign function calls ... OK * checking R code for possible problems ... NOTE ComputePValues: no visible binding for global variable ‘motif’ ComputePValues: no visible binding for global variable ‘snpid’ ComputePValues: no visible binding for global variable ‘snpbase’ ComputePValues: no visible binding for global variable ‘pval_ref’ ComputePValues: no visible binding for global variable ‘pval_snp’ ComputePValues: no visible binding for global variable ‘pval_cond_ref’ ComputePValues: no visible binding for global variable ‘pval_cond_snp’ ComputePValues: no visible binding for global variable ‘pval_diff’ ComputePValues: no visible binding for global variable ‘pval_rank’ LoadSNPData: no visible binding for global variable ‘IUPAC_CODE_MAP’ MatchSubsequence: no visible binding for global variable ‘motif’ MatchSubsequence: no visible binding for global variable ‘snpid’ MatchSubsequence: no visible binding for global variable ‘snpbase’ MatchSubsequence: no visible binding for global variable ‘len_seq’ MatchSubsequence: no visible binding for global variable ‘ref_seq’ MatchSubsequence : <anonymous>: no visible binding for global variable ‘motif’ dtMotifMatch: no visible binding for global variable ‘ref_seq’ dtMotifMatch: no visible binding for global variable ‘len_seq’ dtMotifMatch: no visible binding for global variable ‘snp_ref_start’ dtMotifMatch: no visible binding for global variable ‘ref_start’ dtMotifMatch: no visible binding for global variable ‘snp_start’ dtMotifMatch: no visible binding for global variable ‘snp_ref_end’ dtMotifMatch: no visible binding for global variable ‘ref_end’ dtMotifMatch: no visible binding for global variable ‘snp_end’ dtMotifMatch: no visible binding for global variable ‘snp_ref_length’ dtMotifMatch: no visible binding for global variable ‘ref_aug_match_seq_forward’ dtMotifMatch: no visible binding for global variable ‘ref_aug_match_seq_reverse’ dtMotifMatch: no visible binding for global variable ‘snp_aug_match_seq_forward’ dtMotifMatch: no visible binding for global variable ‘snp_seq’ dtMotifMatch: no visible binding for global variable ‘snp_aug_match_seq_reverse’ dtMotifMatch: no visible binding for global variable ‘ref_strand’ dtMotifMatch: no visible binding for global variable ‘ref_location’ dtMotifMatch: no visible binding for global variable ‘snp_strand’ dtMotifMatch: no visible binding for global variable ‘snp_location’ dtMotifMatch: no visible binding for global variable ‘ref_extra_pwm_left’ dtMotifMatch: no visible binding for global variable ‘ref_extra_pwm_right’ dtMotifMatch: no visible binding for global variable ‘snp_extra_pwm_left’ dtMotifMatch: no visible binding for global variable ‘snp_extra_pwm_right’ dtMotifMatch: no visible binding for global variable ‘snpid’ match_subseq_par: no visible binding for global variable ‘snpid’ match_subseq_par: no visible binding for global variable ‘motif’ match_subseq_par: no visible binding for global variable ‘snpbase’ match_subseq_par: no visible binding for global variable ‘ref_strand’ match_subseq_par: no visible binding for global variable ‘ref_match_seq’ match_subseq_par: no visible binding for global variable ‘ref_seq’ match_subseq_par: no visible binding for global variable ‘ref_start’ match_subseq_par: no visible binding for global variable ‘ref_end’ match_subseq_par: no visible binding for global variable ‘ref_seq_rev’ match_subseq_par: no visible binding for global variable ‘len_seq’ match_subseq_par: no visible binding for global variable ‘snp_strand’ match_subseq_par: no visible binding for global variable ‘snp_match_seq’ match_subseq_par: no visible binding for global variable ‘snp_seq’ match_subseq_par: no visible binding for global variable ‘snp_start’ match_subseq_par: no visible binding for global variable ‘snp_end’ match_subseq_par: no visible binding for global variable ‘snp_seq_rev’ match_subseq_par: no visible binding for global variable ‘snp_seq_ref_match’ match_subseq_par: no visible binding for global variable ‘ref_seq_snp_match’ match_subseq_par: no visible binding for global variable ‘motif_len’ match_subseq_par: no visible binding for global variable ‘log_lik_ref’ match_subseq_par: no visible binding for global variable ‘log_lik_snp’ match_subseq_par: no visible binding for global variable ‘log_lik_ratio’ match_subseq_par: no visible binding for global variable ‘log_enhance_odds’ match_subseq_par: no visible binding for global variable ‘log_reduce_odds’ match_subseq_par: no visible binding for global variable ‘IUPAC’ motif_score_par: no visible binding for global variable ‘motif’ motif_score_par: no visible binding for global variable ‘snpbase’ results_motif_par: no visible binding for global variable ‘p.value’ Undefined global functions or variables: IUPAC IUPAC_CODE_MAP len_seq log_enhance_odds log_lik_ratio log_lik_ref log_lik_snp log_reduce_odds motif motif_len p.value pval_cond_ref pval_cond_snp pval_diff pval_rank pval_ref pval_snp ref_aug_match_seq_forward ref_aug_match_seq_reverse ref_end ref_extra_pwm_left ref_extra_pwm_right ref_location ref_match_seq ref_seq ref_seq_rev ref_seq_snp_match ref_start ref_strand snp_aug_match_seq_forward snp_aug_match_seq_reverse snp_end snp_extra_pwm_left snp_extra_pwm_right snp_location snp_match_seq snp_ref_end snp_ref_length snp_ref_start snp_seq snp_seq_ref_match snp_seq_rev snp_start snp_strand snpbase snpid * checking Rd files ... OK * checking Rd metadata ... OK * checking Rd cross-references ... OK * checking for missing documentation entries ... OK * checking for code/documentation mismatches ... OK * checking Rd \usage sections ... OK * checking Rd contents ... OK * checking for unstated dependencies in examples ... OK * checking contents of ‘data’ directory ... OK * checking data for non-ASCII characters ... OK * checking data for ASCII and uncompressed saves ... OK * checking line endings in C/C++/Fortran sources/headers ... OK * checking compiled code ... NOTE Note: information on .o files is not available File ‘/Library/Frameworks/R.framework/Versions/4.0/Resources/library/atSNP/libs/atSNP.so’: Found ‘_printf’, possibly from ‘printf’ (C) Compiled code should not call entry points which might terminate R nor write to stdout/stderr instead of to the console, nor use Fortran I/O nor system RNGs. The detected symbols are linked into the code but might come from libraries and not actually be called. See ‘Writing portable packages’ in the ‘Writing R Extensions’ manual. * checking files in ‘vignettes’ ... OK * checking examples ... OK * checking for unstated dependencies in ‘tests’ ... OK * checking tests ... Running ‘test.R’ Running ‘test_change.R’ Running ‘test_diff.R’ Running ‘test_is.R’ OK * checking for unstated dependencies in vignettes ... OK * checking package vignettes in ‘inst/doc’ ... OK * checking running R code from vignettes ... SKIPPED * checking re-building of vignette outputs ... SKIPPED * checking PDF version of manual ... OK * DONE Status: 3 NOTEs See ‘/Users/biocbuild/bbs-3.11-bioc/meat/atSNP.Rcheck/00check.log’ for details.
atSNP.Rcheck/00install.out
############################################################################## ############################################################################## ### ### Running command: ### ### /Library/Frameworks/R.framework/Versions/Current/Resources/bin/R CMD INSTALL atSNP ### ############################################################################## ############################################################################## * installing to library ‘/Library/Frameworks/R.framework/Versions/4.0/Resources/library’ * installing *source* package ‘atSNP’ ... ** using staged installation ** libs clang++ -mmacosx-version-min=10.13 -std=gnu++11 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I'/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include' -I/usr/local/include -fPIC -Wall -g -O2 -c ImportanceSample.cpp -o ImportanceSample.o clang++ -mmacosx-version-min=10.13 -std=gnu++11 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I'/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include' -I/usr/local/include -fPIC -Wall -g -O2 -c ImportanceSampleChange.cpp -o ImportanceSampleChange.o clang++ -mmacosx-version-min=10.13 -std=gnu++11 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I'/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include' -I/usr/local/include -fPIC -Wall -g -O2 -c ImportanceSampleDiff.cpp -o ImportanceSampleDiff.o clang++ -mmacosx-version-min=10.13 -std=gnu++11 -I"/Library/Frameworks/R.framework/Resources/include" -DNDEBUG -I'/Library/Frameworks/R.framework/Versions/4.0/Resources/library/Rcpp/include' -I/usr/local/include -fPIC -Wall -g -O2 -c MotifScore.cpp -o MotifScore.o clang++ -mmacosx-version-min=10.13 -std=gnu++11 -dynamiclib -Wl,-headerpad_max_install_names -undefined dynamic_lookup -single_module -multiply_defined suppress -L/Library/Frameworks/R.framework/Resources/lib -L/usr/local/lib -o atSNP.so ImportanceSample.o ImportanceSampleChange.o ImportanceSampleDiff.o MotifScore.o -F/Library/Frameworks/R.framework/.. -framework R -Wl,-framework -Wl,CoreFoundation installing to /Library/Frameworks/R.framework/Versions/4.0/Resources/library/00LOCK-atSNP/00new/atSNP/libs ** R ** data ** inst ** byte-compile and prepare package for lazy loading ** help *** installing help indices ** building package indices ** installing vignettes ** testing if installed package can be loaded from temporary location ** checking absolute paths in shared objects and dynamic libraries ** testing if installed package can be loaded from final location ** testing if installed package keeps a record of temporary installation path * DONE (atSNP)
atSNP.Rcheck/tests/test.Rout
R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(atSNP)
> library(BiocParallel)
> library(testthat)
>
> ## process the data
> data(example)
>
> motif_scores <- ComputeMotifScore(motif_library, snpInfo, ncores = 1)
>
> motif_scores <- MatchSubsequence(motif_scores$snp.tbl, motif_scores$motif.scores, ncores = 1, motif.lib = motif_library)
>
> motif_scores[which(motif_scores$snpid == "rs7412" & motif_scores$motif == "SIX5_disc1"), ]
snpid motif
4 rs7412 SIX5_disc1
ref_seq
4 CTCCTCCGCGATGCCGATGACCTGCAGAAGCGCCTGGCAGTGTACCAGGCCGGGGCCCGCG
snp_seq motif_len
4 CTCCTCCGCGATGCCGATGACCTGCAGAAGTGCCTGGCAGTGTACCAGGCCGGGGCCCGCG 10
ref_start ref_end ref_strand snp_start snp_end snp_strand log_lik_ref
4 29 38 - 22 31 + -42.60672
log_lik_snp log_lik_ratio log_enhance_odds log_reduce_odds IUPAC
4 -38.4083 -4.198418 23.013 -2.917768 GARWTGTAGT
ref_match_seq snp_match_seq ref_seq_snp_match snp_seq_ref_match snpbase
4 GCCAGGCGCT CTGCAGAAGT CTGCAGAAGC GCCAGGCACT T
>
> len_seq <- sapply(motif_scores$ref_seq, nchar)
> snp_pos <- as.integer(len_seq / 2) + 1
>
> i <- which(motif_scores$snpid == "rs7412" & motif_scores$motif == "SIX5_disc1")
>
> test_that("Error: reference bases are not the same as the sequence matrix.", {
+ expect_equal(sum(snpInfo$sequence_matrix[31, ] != snpInfo$ref_base), 0)
+ expect_equal(sum(snpInfo$sequence_matrix[31, ] == snpInfo$snp_base), 0)
+ })
>
> test_that("Error: log_lik_ratio is not correct.", {
+ expect_equal(motif_scores$log_lik_ref - motif_scores$log_lik_snp, motif_scores$log_lik_ratio)
+ })
>
> test_that("Error: log likelihoods are not correct.", {
+
+ log_lik <- sapply(seq(nrow(motif_scores)),
+ function(i) {
+ motif_mat <- motif_library[[motif_scores$motif[i]]]
+ colind<-which(snpInfo$snpids==motif_scores$snpid[i])
+ bases <- snpInfo$sequence_matrix[motif_scores$ref_start[i]:motif_scores$ref_end[i], colind]
+ if(motif_scores$ref_strand[i] == "-")
+ bases <- 5 - rev(bases)
+ log(prod(
+ motif_mat[cbind(seq(nrow(motif_mat)),
+ bases)]))
+ })
+
+ expect_equal(log_lik, motif_scores$log_lik_ref)
+
+ snp_mat <- snpInfo$sequence_matrix
+ snp_mat[cbind(snp_pos, seq(ncol(snp_mat)))] <- snpInfo$snp_base
+ log_lik <- sapply(seq(nrow(motif_scores)),
+ function(i) {
+ motif_mat <- motif_library[[motif_scores$motif[i]]]
+ colind<-which(snpInfo$snpids==motif_scores$snpid[i])
+ bases <- snp_mat[motif_scores$snp_start[i]:motif_scores$snp_end[i], colind]
+ if(motif_scores$snp_strand[i] == "-")
+ bases <- 5 - rev(bases)
+ log(prod(
+ motif_mat[cbind(seq(nrow(motif_mat)),
+ bases)]))
+ })
+
+ expect_equal(log_lik, motif_scores$log_lik_snp)
+ })
>
> test_that("Error: log_enhance_odds not correct.", {
+
+ len_seq <- sapply(motif_scores$ref_seq, nchar)
+ snp_pos <- as.integer(len_seq / 2) + 1
+
+ ## log odds for reduction in binding affinity
+
+ pos_in_pwm <- snp_pos - motif_scores$ref_start + 1
+ neg_ids <- which(motif_scores$ref_strand == "-")
+ pos_in_pwm[neg_ids] <- motif_scores$ref_end[neg_ids]- snp_pos[neg_ids] + 1
+ snp_base <- sapply(substr(motif_scores$snp_seq, snp_pos, snp_pos), function(x) which(c("A", "C", "G", "T") == x))
+ ref_base <- sapply(substr(motif_scores$ref_seq, snp_pos, snp_pos), function(x) which(c("A", "C", "G", "T") == x))
+ snp_base[neg_ids] <- 5 - snp_base[neg_ids]
+ ref_base[neg_ids] <- 5 - ref_base[neg_ids]
+ my_log_reduce_odds <- sapply(seq(nrow(motif_scores)),
+ function(i)
+ log(motif_library[[motif_scores$motif[i]]][pos_in_pwm[i], ref_base[i]]) -
+ log(motif_library[[motif_scores$motif[i]]][pos_in_pwm[i], snp_base[i]])
+ )
+
+ expect_equal(my_log_reduce_odds, motif_scores$log_reduce_odds)
+
+ ## log odds in enhancing binding affinity
+
+ pos_in_pwm <- snp_pos - motif_scores$snp_start + 1
+ neg_ids <- which(motif_scores$snp_strand == "-")
+ pos_in_pwm[neg_ids] <- motif_scores$snp_end[neg_ids]- snp_pos[neg_ids] + 1
+ snp_base <- sapply(substr(motif_scores$snp_seq, snp_pos, snp_pos), function(x) which(c("A", "C", "G", "T") == x))
+ ref_base <- sapply(substr(motif_scores$ref_seq, snp_pos, snp_pos), function(x) which(c("A", "C", "G", "T") == x))
+ snp_base[neg_ids] <- 5 - snp_base[neg_ids]
+ ref_base[neg_ids] <- 5 - ref_base[neg_ids]
+ my_log_enhance_odds <- sapply(seq(nrow(motif_scores)),
+ function(i)
+ log(motif_library[[motif_scores$motif[i]]][pos_in_pwm[i], snp_base[i]]) -
+ log(motif_library[[motif_scores$motif[i]]][pos_in_pwm[i], ref_base[i]])
+ )
+
+ expect_equal(my_log_enhance_odds, motif_scores$log_enhance_odds)
+
+
+ })
>
> test_that("Error: the maximum log likelihood computation is not correct.", {
+
+ snp_mat <- snpInfo$sequence_matrix
+ snp_mat[cbind(snp_pos, seq(ncol(snp_mat)))] <- snpInfo$snp_base
+
+ .findMaxLog <- function(seq_vec, pwm) {
+ snp_pos <- as.integer(length(seq_vec) / 2) + 1
+ start_pos <- snp_pos - nrow(pwm) + 1
+ end_pos <- snp_pos
+ rev_seq <- 5 - rev(seq_vec)
+
+ maxLogProb <- -Inf
+ for(i in start_pos : end_pos) {
+ LogProb <- log(prod(pwm[cbind(seq(nrow(pwm)),
+ seq_vec[i - 1 + seq(nrow(pwm))])]))
+ if(LogProb > maxLogProb)
+ maxLogProb <- LogProb
+ }
+ for(i in start_pos : end_pos) {
+ LogProb <- log(prod(pwm[cbind(seq(nrow(pwm)),
+ rev_seq[i - 1 + seq(nrow(pwm))])]))
+ if(LogProb > maxLogProb)
+ maxLogProb <- LogProb
+ }
+ return(maxLogProb)
+ }
+
+ ## find the maximum log likelihood on the reference sequence
+ my_log_lik_ref <- sapply(seq(nrow(motif_scores)),
+ function(x) {
+ colind<-which(snpInfo$snpids==motif_scores$snpid[x])
+ seq_vec<- snpInfo$sequence_matrix[, colind]
+ pwm <- motif_library[[motif_scores$motif[x]]]
+ return(.findMaxLog(seq_vec, pwm))
+ })
+
+ ## find the maximum log likelihood on the SNP sequence
+
+ my_log_lik_snp <- sapply(seq(nrow(motif_scores)),
+ function(x) {
+ colind<-which(snpInfo$snpids==motif_scores$snpid[x]) #ADDED
+ seq_vec<- snp_mat[, colind]
+ pwm <- motif_library[[motif_scores$motif[x]]]
+ return(.findMaxLog(seq_vec, pwm))
+ })
+
+ expect_equal(my_log_lik_ref, motif_scores$log_lik_ref)
+ expect_equal(my_log_lik_snp, motif_scores$log_lik_snp)
+
+ })
>
> proc.time()
user system elapsed
12.522 0.785 13.474
atSNP.Rcheck/tests/test_change.Rout
R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(atSNP)
> library(BiocParallel)
> library(testthat)
> data(example)
>
> trans_mat <- matrix(rep(snpInfo$prior, each = 4), nrow = 4)
> test_pwm <- motif_library$SIX5_disc1
> scores <- as.matrix(motif_scores$motif.scores[3:4, 4:5])
> score_diff <- abs(scores[,2]-scores[,1])
>
> pval_a <- .Call("test_p_value", test_pwm, snpInfo$prior, snpInfo$transition, scores, 0.15, 100)
> pval_ratio <- abs(log(pval_a[seq(nrow(scores)),1]) - log(pval_a[seq(nrow(scores)) + nrow(scores), 1]))
>
> test_score <- test_pwm
> for(i in seq(nrow(test_score))) {
+ for(j in seq(ncol(test_score))) {
+ test_score[i, j] <- exp(mean(log(test_pwm[i, j] / test_pwm[i, -j])))
+ }
+ }
>
> adj_mat <- test_pwm + 0.25
> motif_len <- nrow(test_pwm)
>
> ## these are functions for this test only
> drawonesample <- function(theta) {
+ prob_start <- rev(rowSums(test_score ^ theta) / rowSums(adj_mat))
+ id <- sample(seq(motif_len), 1, prob = prob_start)
+ sample <- sample(1:4, 2 * motif_len - 1, replace = TRUE, prob = snpInfo$prior)
+ delta <- adj_mat
+ delta[motif_len - id + 1, ] <- test_score[motif_len - id + 1, ] ^ theta
+ sample[id - 1 + seq(motif_len)] <- apply(delta, 1, function(x) sample(seq(4), 1, prob = x))
+ ## compute weight
+ sc <- 0
+ for(s in seq(motif_len)) {
+ delta <- adj_mat
+ delta[motif_len + 1 - s, ] <- test_score[motif_len + 1 - s, ] ^ theta
+ sc <- sc + prod(delta[cbind(seq(motif_len), sample[s - 1 + seq(motif_len)])]) /
+ prod(snpInfo$prior[sample[s - 1 + seq(motif_len)]])
+ }
+ sample <- c(sample, id, sc)
+ return(sample)
+ }
> jointprob <- function(x) prod(test_pwm[cbind(seq(motif_len), x)])
> maxjointprob <- function(x) {
+ maxp <- -Inf
+ p <- -Inf
+ for(i in 1:motif_len) {
+ p <- jointprob(x[i:(i+motif_len - 1)])
+ if(p > maxp)
+ maxp <- p
+ }
+ for(i in 1:motif_len) {
+ p <- jointprob(5 - x[(i+motif_len - 1):i])
+ if(p > maxp)
+ maxp <- p
+ }
+ return(maxp)
+ }
> get_freq <- function(sample) {
+ emp_freq <- matrix(0, nrow = 2 * motif_len - 1, ncol = 4)
+ for(i in seq(2 * motif_len - 1)) {
+ for(j in seq(4)) {
+ emp_freq[i, j] <- sum(sample[i, ] == j - 1)
+ }
+ }
+ emp_freq <- emp_freq / rowSums(emp_freq)
+ return(emp_freq)
+ }
>
> test_that("Error: quantile function computing are not equivalent.", {
+ for(p in c(0.01, 0.1, 0.5, 0.9, 0.99) ) {
+ delta <- .Call("test_find_percentile_change", score_diff, p, package = "atSNP")
+ delta.r <- as.double(sort(abs(scores[,2]-scores[,1]))[ceiling((1 - p) * (nrow(scores)))])
+ expect_equal(delta, delta.r)
+ }
+ })
>
> test_that("Error: the scores for samples are not equivalent.", {
+ p <- 0.1
+ delta <- .Call("test_find_percentile_change", score_diff, p, package = "atSNP")
+ theta <- .Call("test_find_theta_change", test_score, adj_mat, delta, package = "atSNP")
+ ## Use R code to generate a random sample
+ for(i in seq(10)) {
+ sample <- drawonesample(theta)
+ sample_score <- .Call("test_compute_sample_score_change", test_pwm, test_score, adj_mat, sample[seq(2 * motif_len - 1)] - 1, snpInfo$prior, trans_mat, sample[2 * motif_len] - 1, theta, package = "atSNP")
+ expect_equal(sample[2 * motif_len + 1], sample_score[1])
+ sample1 <- sample2 <- sample3 <- sample
+ sample1[motif_len] <- seq(4)[-sample[motif_len]][1]
+ sample2[motif_len] <- seq(4)[-sample[motif_len]][2]
+ sample3[motif_len] <- seq(4)[-sample[motif_len]][3]
+ sample_score_r <- log(maxjointprob(sample[seq(2 * motif_len - 1)])) -
+ log(c(maxjointprob(sample1[seq(2 * motif_len - 1)]),
+ maxjointprob(sample2[seq(2 * motif_len - 1)]),
+ maxjointprob(sample3[seq(2 * motif_len - 1)])))
+ expect_equal(sample_score_r, sample_score[2:4])
+ }
+
+ ## Use C code to generate a random sample
+ for(i in seq(10)) {
+ sample <- .Call("test_importance_sample_change", test_score, snpInfo$prior, trans_mat, test_pwm, theta, package = "atSNP")
+ start_pos <- sample[2 * motif_len] + 1
+ adj_score <- 0
+ for(s in seq_len(motif_len)) {
+ adj_s <- sum(log(adj_mat[cbind(seq(motif_len), sample[s - 1 + seq(motif_len)] + 1)]) -
+ log(snpInfo$prior[sample[s - 1 + seq(motif_len)] + 1]))
+ adj_s <- adj_s + theta * log(test_score[motif_len + 1 - s, sample[motif_len] + 1]) -
+ log(adj_mat[motif_len + 1 - s, sample[motif_len] + 1])
+ adj_score <- adj_score + exp(adj_s)
+ }
+ sample_score <- .Call("test_compute_sample_score_change", test_pwm, test_score, adj_mat, sample[seq(2 * motif_len - 1)], snpInfo$prior, trans_mat, sample[2 * motif_len], theta, package = "atSNP")
+ expect_equal(adj_score, sample_score[1])
+ }
+ })
>
> test_that("Error: compute the normalizing constant.", {
+ ## parameters
+ for(p in seq(9) / 10) {
+ delta <- .Call("test_find_percentile_change", score_diff, p, package = "atSNP")
+ theta <- .Call("test_find_theta_change", test_score, adj_mat, delta, package = "atSNP")
+ const <- .Call("test_func_delta_change", test_score, adj_mat, theta, package = "atSNP")
+ ## in R
+ adj_sum <- rowSums(adj_mat)
+ wei_sum <- rowSums(test_score ^ theta)
+ const.r <- prod(adj_sum) * sum(wei_sum / adj_sum)
+ expect_equal(const, const.r)
+ }
+ })
>
> test_that("Error: sample distributions are not expected.", {
+ ## parameters
+ p <- 0.1
+ delta <- .Call("test_find_percentile_change", score_diff, p, package = "atSNP")
+ theta <- .Call("test_find_theta_change", test_score, adj_mat, delta, package = "atSNP")
+ prob_start <- rev(rowSums(test_score ^ theta) / rowSums(adj_mat))
+ ## construct the delta matrix
+ delta <- matrix(1, nrow = 4 * motif_len, ncol = 2 * motif_len - 1)
+ for(pos in seq(motif_len)) {
+ delta[seq(4) + 4 * (pos - 1), ] <- snpInfo$prior
+ delta[seq(4) + 4 * (pos - 1), pos - 1 + seq(motif_len)] <- t(test_pwm)
+ delta[seq(4) + 4 * (pos - 1), motif_len] <- test_score[motif_len + 1 - pos, ] ^ theta
+ delta[seq(4) + 4 * (pos - 1), ] <- delta[seq(4) + 4 * (pos - 1),] / rep(colSums(delta[seq(4) + 4 * (pos - 1), ]), each = 4)
+ }
+ target_freq <- matrix(0, nrow = 4, ncol = 2 * motif_len - 1)
+ for(pos in seq(motif_len)) {
+ target_freq <- target_freq + delta[seq(4) + 4 * (pos - 1), ] * prob_start[pos]
+ }
+ target_freq <- t(target_freq)
+ target_freq <- target_freq / rowSums(target_freq)
+
+ results_i <- function(i) {
+ ## generate 100 samples
+ sample1 <- sapply(seq(100), function(x)
+ .Call("test_importance_sample_change",
+ adj_mat, snpInfo$prior, trans_mat, test_score, theta, package = "atSNP"))
+ emp_freq1 <- get_freq(sample1)
+ sample2 <- sapply(rep(theta, 100), drawonesample)
+ emp_freq2 <- get_freq(sample2 - 1)
+ ## print(rbind(emp_freq1[10, ], emp_freq2[10, ], target_freq[10, ]))
+ max(abs(emp_freq1 - target_freq)) > max(abs(emp_freq2 - target_freq))
+ }
+
+ if(Sys.info()[["sysname"]] == "Windows"){
+ snow <- SnowParam(workers = 1, type = "SOCK")
+ results<-bpmapply(results_i, seq(20), BPPARAM = snow,SIMPLIFY = FALSE)
+ }else{
+ results<-bpmapply(results_i, seq(20), BPPARAM = MulticoreParam(workers = 1),
+ SIMPLIFY = FALSE)
+ }
+
+ print(sum(unlist(results)))
+ print(pbinom(sum(unlist(results)), size = 20, prob = 0.5))
+ })
[1] 11
[1] 0.7482777
>
> test_that("Error: the chosen pvalues should have the smaller variance.", {
+ .structure_diff <- function(pval_mat) {
+ id <- apply(pval_mat[, c(2, 4)], 1, which.min)
+ return(cbind(pval_mat[, c(1, 3)][cbind(seq_along(id), id)],
+ pval_mat[, c(2, 4)][cbind(seq_along(id), id)]))
+ }
+ for(p in c(0.05, 0.1, 0.2, 0.5)) {
+ p_values <- .Call("test_p_value_change", test_pwm, test_score, adj_mat, snpInfo$prior, snpInfo$transition, score_diff, pval_ratio, quantile(score_diff, 1 - p), 100, package = "atSNP")$score
+ p_values_s <- .structure_diff(p_values)
+ expect_equal(p_values_s[, 2], apply(p_values[, c(2, 4)], 1, min))
+ }
+ })
>
> proc.time()
user system elapsed
16.044 0.761 16.857
atSNP.Rcheck/tests/test_diff.Rout
R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(atSNP)
> library(BiocParallel)
> library(testthat)
> data(example)
>
> trans_mat <- matrix(rep(snpInfo$prior, each = 4), nrow = 4)
> test_pwm <- motif_library$SIX5_disc1
> scores <- as.matrix(motif_scores$motif.scores[3:4, 4:5])
> score_diff <- abs(scores[,2]-scores[,1])
>
> test_score <- test_pwm
> for(i in seq(nrow(test_score))) {
+ for(j in seq(ncol(test_score))) {
+ test_score[i, j] <- exp(mean(log(test_pwm[i, j] / test_pwm[i, -j])))
+ }
+ }
>
> adj_mat <- test_pwm + rowMeans(test_pwm)
> motif_len <- nrow(test_pwm)
>
> ## these are functions for this test only
> drawonesample <- function(theta) {
+ prob_start <- sapply(seq(motif_len),
+ function(j)
+ sum(snpInfo$prior * test_score[motif_len + 1 - j, ] ^ theta *
+ adj_mat[motif_len + 1 - j, ]) /
+ sum(snpInfo$prior * adj_mat[motif_len + 1 - j, ])
+ )
+ id <- sample(seq(motif_len), 1, prob = prob_start)
+ sample <- sample(1:4, 2 * motif_len - 1, replace = TRUE, prob = snpInfo$prior)
+ delta <- adj_mat
+ delta[motif_len + 1 - id, ] <- delta[motif_len + 1 - id, ] * test_score[motif_len + 1 - id, ] ^ theta
+ sample[id - 1 + seq(motif_len)] <- apply(delta, 1, function(x)
+ sample(seq(4), 1, prob = x * snpInfo$prior))
+ sc <- 0
+ for(s in seq(motif_len)) {
+ delta <- adj_mat
+ delta[motif_len + 1 - s, ] <- delta[motif_len + 1 - s, ] * test_score[motif_len + 1 - s, ] ^ theta
+ sc <- sc + prod(delta[cbind(seq(motif_len), sample[s - 1 + seq(motif_len)])])
+ }
+ sample <- c(sample, id, sc)
+ return(sample)
+ }
> jointprob <- function(x) prod(test_pwm[cbind(seq(motif_len), x)])
> maxjointprob <- function(x) {
+ maxp <- -Inf
+ p <- -Inf
+ for(i in 1:motif_len) {
+ p <- jointprob(x[i:(i+motif_len - 1)])
+ if(p > maxp)
+ maxp <- p
+ }
+ for(i in 1:motif_len) {
+ p <- jointprob(5 - x[(i+motif_len - 1):i])
+ if(p > maxp)
+ maxp <- p
+ }
+ return(maxp)
+ }
> get_freq <- function(sample) {
+ emp_freq <- matrix(0, nrow = 2 * motif_len - 1, ncol = 4)
+ for(i in seq(2 * motif_len - 1)) {
+ for(j in seq(4)) {
+ emp_freq[i, j] <- sum(sample[i, ] == j - 1)
+ }
+ }
+ emp_freq <- emp_freq / rowSums(emp_freq)
+ return(emp_freq)
+ }
>
> test_that("Error: quantile function computing are not equivalent.", {
+ for(p in c(0.01, 0.1, 0.5, 0.9, 0.99)) {
+ delta <- .Call("test_find_percentile_diff", score_diff, p, package = "atSNP")
+ delta.r <- as.double(sort(abs(scores[,2]-scores[,1]))[ceiling((1 - p) * (nrow(scores)))])
+ expect_equal(delta, delta.r)
+ }
+ })
>
> test_that("Error: the scores for samples are not equivalent.", {
+ p <- 0.1
+ delta <- .Call("test_find_percentile_diff", score_diff, p, package = "atSNP")
+ theta <- .Call("test_find_theta_diff", test_score, adj_mat, snpInfo$prior, snpInfo$transition, delta, package = "atSNP")
+ ## Use R code to generate a random sample
+ for(i in seq(10)) {
+ sample <- drawonesample(theta)
+ sample_score <- .Call("test_compute_sample_score_diff", test_pwm, test_score, adj_mat, sample[seq(2 * motif_len - 1)] - 1, sample[2 * motif_len] - 1, theta, package = "atSNP")
+ expect_equal(sample[2 * motif_len + 1], sample_score[1])
+ sample1 <- sample2 <- sample3 <- sample
+ sample1[motif_len] <- seq(4)[-sample[motif_len]][1]
+ sample2[motif_len] <- seq(4)[-sample[motif_len]][2]
+ sample3[motif_len] <- seq(4)[-sample[motif_len]][3]
+ sample_score_r <- log(maxjointprob(sample[seq(2 * motif_len - 1)])) -
+ log(c(maxjointprob(sample1[seq(2 * motif_len - 1)]),
+ maxjointprob(sample2[seq(2 * motif_len - 1)]),
+ maxjointprob(sample3[seq(2 * motif_len - 1)])))
+ expect_equal(sample_score_r, sample_score[-1])
+ }
+
+ ## Use C code to generate a random sample
+ delta <- matrix(1, nrow = 4 * motif_len, ncol = 2 * motif_len - 1)
+ for(pos in seq(motif_len)) {
+ for(j in (pos + motif_len - 1) : 1) {
+ if(j < pos + motif_len - 1) {
+ delta[4 * (pos - 1) + seq(4), j] <- sum(snpInfo$prior * delta[4 * (pos - 1) + seq(4), j + 1])
+ }
+ if(j >= pos) {
+ delta[4 * (pos - 1) + seq(4), j] <- delta[4 * (pos - 1) + seq(4), j] * adj_mat[j - pos + 1, ]
+ }
+ if(j == motif_len) {
+ delta[4 * (pos - 1) + seq(4), j] <- delta[4 * (pos - 1) + seq(4), j] * test_score[j - pos + 1, ] ^ theta
+ }
+ }
+ }
+ for(i in seq(10)) {
+ sample <- .Call("test_importance_sample_diff", delta, snpInfo$prior, trans_mat, test_pwm, theta, package = "atSNP")
+ start_pos <- sample[2 * motif_len] + 1
+ adj_score <- 0
+ for(s in seq_len(motif_len)) {
+ adj_s <- sum(log(adj_mat[cbind(seq(motif_len), sample[s - 1 + seq(motif_len)] + 1)]))
+ adj_s <- adj_s + theta * log(test_score[motif_len + 1 - s, sample[motif_len] + 1])
+ adj_score <- adj_score + exp(adj_s)
+ }
+ sample_score <- .Call("test_compute_sample_score_diff", test_pwm, test_score, adj_mat, sample[seq(2 * motif_len - 1)], sample[2 * motif_len], theta, package = "atSNP")
+ expect_equal(adj_score, sample_score[1])
+ }
+ })
>
> test_that("Error: compute the normalizing constant.", {
+
+ ## parameters
+ p <- 0.1
+ delta <- .Call("test_find_percentile_diff", score_diff, p, package = "atSNP")
+ theta <- .Call("test_find_theta_diff", test_score, adj_mat, snpInfo$prior, snpInfo$transition, delta, package = "atSNP")
+
+ ##
+ const <- .Call("test_func_delta_diff", test_score, adj_mat, snpInfo$prior, trans_mat, theta, package = "atSNP")
+
+ prob_start <- sapply(seq(motif_len),
+ function(j)
+ sum(snpInfo$prior * test_score[motif_len + 1 - j, ] ^ theta *
+ adj_mat[motif_len + 1 - j, ]) /
+ sum(snpInfo$prior * adj_mat[motif_len + 1 - j, ])
+ )
+
+ const.r <- prod(colSums(snpInfo$prior * t(adj_mat))) * sum(prob_start)
+ expect_equal(const, const.r)
+ })
>
> test_that("Error: sample distributions are not expected.", {
+
+ ## parameters
+ p <- 0.1
+ delta <- .Call("test_find_percentile_diff", score_diff, p, package = "atSNP")
+ theta <- .Call("test_find_theta_diff", test_score, adj_mat, snpInfo$prior, snpInfo$transition, delta, package = "atSNP")
+
+ ## construct the delta matrix
+ delta <- matrix(1, nrow = 4 * motif_len, ncol = 2 * motif_len - 1)
+ for(pos in seq(motif_len)) {
+ for(j in (pos + motif_len - 1) : 1) {
+ if(j < pos + motif_len - 1) {
+ delta[4 * (pos - 1) + seq(4), j] <- sum(snpInfo$prior * delta[4 * (pos - 1) + seq(4), j + 1])
+ }
+ if(j >= pos) {
+ delta[4 * (pos - 1) + seq(4), j] <- delta[4 * (pos - 1) + seq(4), j] * adj_mat[j - pos + 1, ]
+ }
+ if(j == motif_len) {
+ delta[4 * (pos - 1) + seq(4), j] <- delta[4 * (pos - 1) + seq(4), j] * test_score[j - pos + 1, ] ^ theta
+ }
+ }
+ }
+
+ target_freq <- matrix(0, nrow = 4, ncol = 2 * motif_len - 1)
+
+ mat <- snpInfo$prior * matrix(delta[, 1], nrow = 4)
+ wei <- colSums(mat)
+ for(j in seq(2 * motif_len - 1)) {
+ for(pos in seq(motif_len)) {
+ tmp <- delta[seq(4) + 4 * (pos - 1), j] * snpInfo$prior
+ target_freq[, j] <- target_freq[, j] + tmp / sum(tmp) * wei[pos]
+ }
+ }
+ target_freq <- t(target_freq)
+ target_freq <- target_freq / rowSums(target_freq)
+
+ results_i <- function(i) {
+ ## generate 100 samples
+ sample1 <- sapply(seq(100), function(x)
+ .Call("test_importance_sample_diff",
+ delta, snpInfo$prior, trans_mat, test_score, theta, package = "atSNP"))
+ emp_freq1 <- get_freq(sample1)
+
+ sample2 <- sapply(rep(theta, 100), drawonesample)
+ emp_freq2 <- get_freq(sample2 - 1)
+
+ ## print(rbind(emp_freq1[10, ], emp_freq2[10, ], target_freq[10, ]))
+ max(abs(emp_freq1 - target_freq)) > max(abs(emp_freq2 - target_freq))
+ }
+
+ if(Sys.info()[["sysname"]] == "Windows"){
+ snow <- SnowParam(workers = 1, type = "SOCK")
+ results<-bpmapply(results_i, seq(20), BPPARAM = snow,SIMPLIFY = FALSE)
+ }else{
+ results<-bpmapply(results_i, seq(20), BPPARAM = MulticoreParam(workers = 1),
+ SIMPLIFY = FALSE)
+ }
+
+ print(sum(unlist(results)))
+
+ print(pbinom(sum(unlist(results)), size = 20, prob = 0.5))
+
+ })
[1] 9
[1] 0.4119015
>
> test_that("Error: the chosen pvalues should have the smaller variance.", {
+
+ .structure_diff <- function(pval_mat) {
+ id <- apply(pval_mat[, c(2, 4)], 1, which.min)
+ return(cbind(pval_mat[, c(1, 3)][cbind(seq_along(id), id)],
+ pval_mat[, c(2, 4)][cbind(seq_along(id), id)]))
+ }
+
+ for(p in c(0.05, 0.1, 0.2, 0.5)) {
+ p_values <- .Call("test_p_value_diff", test_pwm, test_score, adj_mat, snpInfo$prior, snpInfo$transition, score_diff, quantile(score_diff, 1 - p), 100, package = "atSNP")
+ p_values_s <- .structure_diff(p_values)
+ expect_equal(p_values_s[, 2], apply(p_values[, c(2, 4)], 1, min))
+ }
+ })
>
> proc.time()
user system elapsed
17.149 0.718 17.869
atSNP.Rcheck/tests/test_is.Rout
R version 4.0.3 (2020-10-10) -- "Bunny-Wunnies Freak Out"
Copyright (C) 2020 The R Foundation for Statistical Computing
Platform: x86_64-apple-darwin17.0 (64-bit)
R is free software and comes with ABSOLUTELY NO WARRANTY.
You are welcome to redistribute it under certain conditions.
Type 'license()' or 'licence()' for distribution details.
R is a collaborative project with many contributors.
Type 'contributors()' for more information and
'citation()' on how to cite R or R packages in publications.
Type 'demo()' for some demos, 'help()' for on-line help, or
'help.start()' for an HTML browser interface to help.
Type 'q()' to quit R.
> library(atSNP)
> library(BiocParallel)
> library(testthat)
> data(example)
>
> trans_mat <- matrix(rep(snpInfo$prior, each = 4), nrow = 4)
> test_pwm <- motif_library$SIX5_disc1
> scores <- as.matrix(motif_scores$motif.scores[3:4, 4:5])
>
> motif_len <- nrow(test_pwm)
>
> ## these are functions for this test only
> drawonesample <- function(theta) {
+ delta <- snpInfo$prior * t(test_pwm ^ theta)
+ delta <- delta / rep(colSums(delta), each = 4)
+ sample <- sample(1:4, 2 * motif_len - 1, replace = TRUE, prob = snpInfo$prior)
+ id <- sample(seq(motif_len), 1)
+ sample[id : (id + motif_len - 1)] <- apply(delta, 2, function(x) sample(1:4, 1, prob = x))
+ sc <- s_cond <- 0
+ for(s in seq(motif_len)) {
+ sc <- sc + prod(test_pwm[cbind(seq(motif_len),
+ sample[s : (s + motif_len - 1)])]) ^ theta
+ }
+ s_cond <- prod(test_pwm[cbind(seq(motif_len),
+ sample[id : (id + motif_len - 1)])]) ^ theta
+ sample <- c(sample, id, sc, s_cond)
+ return(sample)
+ }
> jointprob <- function(x) prod(test_pwm[cbind(seq(motif_len), x)])
> maxjointprob <- function(x) {
+ maxp <- -Inf
+ p <- -Inf
+ for(i in 1:motif_len) {
+ p <- jointprob(x[i:(i+motif_len - 1)])
+ if(p > maxp)
+ maxp <- p
+ }
+ for(i in 1:motif_len) {
+ p <- jointprob(5 - x[(i+motif_len - 1):i])
+ if(p > maxp)
+ maxp <- p
+ }
+ return(maxp)
+ }
> get_freq <- function(sample) {
+ ids <- cbind(
+ rep(sample[motif_len * 2, ], each = motif_len) + seq(motif_len),
+ rep(seq(100), each = motif_len))
+ sample_motif <- matrix(sample[ids], nrow = motif_len) + 1
+ emp_freq <- matrix(0, nrow = motif_len, ncol = 4)
+ for(i in seq(motif_len)) {
+ for(j in seq(4)) {
+ emp_freq[i, j] <- sum(sample_motif[i, ] == j)
+ }
+ }
+ emp_freq <- emp_freq / rowSums(emp_freq)
+ return(emp_freq)
+ }
>
> test_that("Error: quantile function computing are not equivalent.", {
+ for(p in c(0.01, 0.1, 0.5, 0.9, 0.99)) {
+ delta <- .Call("test_find_percentile", c(scores), p, package = "atSNP")
+ delta.r <- -sort(-c(scores))[as.integer(p * length(scores)) + 1]
+ expect_equal(delta, delta.r)
+ }
+ })
>
> test_that("Error: the scores for samples are not equivalent.", {
+ p <- 0.01
+ delta <- .Call("test_find_percentile", scores, p, package = "atSNP")
+ theta <- .Call("test_find_theta", test_pwm, snpInfo$prior, snpInfo$transition, delta, package = "atSNP")
+ ## Use R code to generate a random sample
+ for(i in seq(10)) {
+ sample <- drawonesample(theta)
+ sample_score <- .Call("test_compute_sample_score", test_pwm, sample[seq(2 * motif_len - 1)] - 1, sample[motif_len * 2] - 1, theta, package = "atSNP")
+ expect_equal(sample[2 * motif_len + 1], sample_score[2])
+ expect_equal(sample[2 * motif_len + 2], sample_score[3])
+ }
+ ## Use C code to generate a random sample
+ for(i in seq(10)) {
+ delta <- t(test_pwm ^ theta)
+ delta <- cbind(matrix(
+ sum(snpInfo$prior * delta[, 1]),
+ nrow = 4, ncol = motif_len - 1), delta)
+ sample <- .Call("test_importance_sample", delta, snpInfo$prior, trans_mat, test_pwm, theta, package = "atSNP")
+ start_pos <- sample[motif_len * 2]
+ adj_score <- 0
+ for(s in seq(motif_len) - 1) {
+ adj_score <- adj_score + prod(test_pwm[cbind(seq(motif_len),
+ sample[s + seq(motif_len)] + 1)]) ^ theta
+ }
+ adj_score_cond <- prod(test_pwm[cbind(seq(motif_len), sample[start_pos + seq(motif_len)] + 1)]) ^ theta
+ sample_score <- .Call("test_compute_sample_score", test_pwm, sample[seq(2 * motif_len - 1)], sample[motif_len * 2], theta, package = "atSNP")
+ expect_equal(adj_score, sample_score[2])
+ expect_equal(adj_score_cond, sample_score[3])
+ }
+ })
>
> test_that("Error: compute the normalizing constant.", {
+ ## parameters
+ p <- 0.01
+ delta <- .Call("test_find_percentile", scores, p, package = "atSNP")
+ theta <- .Call("test_find_theta", test_pwm, snpInfo$prior, snpInfo$transition, delta, package = "atSNP")
+ ##
+ const <- .Call("test_func_delta", test_pwm, snpInfo$prior, trans_mat, theta, package = "atSNP")
+ const.r <- prod(colSums(snpInfo$prior * t(test_pwm) ^ theta)) * motif_len
+ expect_equal(abs(const - const.r) / const < 1e-5, TRUE)
+ })
>
> test_that("Error: sample distributions are not expected.", {
+ ## parameters
+ p <- 0.1
+ delta <- .Call("test_find_percentile", scores, p, package = "atSNP")
+ theta <- .Call("test_find_theta", test_pwm, snpInfo$prior, trans_mat, delta, package = "atSNP")
+ delta <- t(test_pwm ^ theta)
+ delta <- cbind(matrix(
+ sum(snpInfo$prior * delta[, 1]),
+ nrow = 4, ncol = motif_len - 1), delta)
+
+ results_i <- function(i) {
+ ## generate 100 samples
+ sample <- sapply(seq(100), function(x)
+ .Call("test_importance_sample",
+ delta, snpInfo$prior, trans_mat, test_pwm, theta, package = "atSNP"))
+ emp_freq1 <- get_freq(sample)
+ target_freq <- test_pwm ^ theta * snpInfo$prior
+ target_freq <- target_freq / rowSums(target_freq)
+ ## generate samples in R
+ sample <- sapply(rep(theta, 100), drawonesample)
+ emp_freq2 <- get_freq(sample[seq(2 * motif_len), ] - 1)
+ max(abs(emp_freq1 - target_freq)) > max(abs(emp_freq2 - target_freq))
+ }
+
+ if(Sys.info()[["sysname"]] == "Windows"){
+ snow <- SnowParam(workers = 1, type = "SOCK")
+ results<-bpmapply(results_i, seq(20), BPPARAM = snow,SIMPLIFY = FALSE)
+ }else{
+ results<-bpmapply(results_i, seq(20), BPPARAM = MulticoreParam(workers = 1),
+ SIMPLIFY = FALSE)
+ }
+
+ print(sum(unlist(results)))
+ print(pbinom(sum(unlist(results)), size = 20, prob = 0.5))
+ })
[1] 18
[1] 0.99998
>
> test_that("Error: the chosen pvalues should have the smaller variance.", {
+ .structure <- function(pval_mat) {
+ id1 <- apply(pval_mat[, c(2, 4)], 1, which.min)
+ return(cbind(
+ pval_mat[, c(1, 3)][cbind(seq_along(id1), id1)],
+ pval_mat[, c(2, 4)][cbind(seq_along(id1), id1)])
+ )
+ }
+ for(p in c(0.01, 0.05, 0.1)) {
+ theta <- .Call("test_find_theta", test_pwm, snpInfo$prior, trans_mat, quantile(c(scores), 1 - p), package = "atSNP")
+ p_values <- .Call("test_p_value", test_pwm, snpInfo$prior, snpInfo$transition, c(scores), theta, 100, package = "atSNP")
+ p_values_s <- .structure(p_values)
+ expect_equal(p_values_s[, 2], apply(p_values[, c(2, 4)], 1, min))
+ }
+ })
>
> proc.time()
user system elapsed
18.388 0.756 19.137
atSNP.Rcheck/atSNP-Ex.timings
| name | user | system | elapsed | |
| ComputeMotifScore | 0.134 | 0.127 | 0.186 | |
| ComputePValues | 0.354 | 0.133 | 0.314 | |
| GetIUPACSequence | 0.003 | 0.002 | 0.004 | |
| LoadFastaData | 0.592 | 0.077 | 0.684 | |
| LoadMotifLibrary | 0.433 | 0.048 | 0.485 | |
| LoadSNPData | 0.000 | 0.000 | 0.001 | |
| MatchSubsequence | 0.110 | 0.077 | 0.125 | |
| dtMotifMatch | 0.093 | 0.044 | 0.152 | |
| plotMotifMatch | 1.486 | 0.409 | 1.858 | |