#' Fisher test
#'
#' @param obj an object
#'
#' @param ... other variables
#' @return the output will be a test_snps (data.table) object including 3 columns:
#' \dQuote{snp_pos}, \dQuote{case_control}, and \dQuote{value} which is
#' a p-values.
#' @export
fisher_test <- function(obj, ...) {
UseMethod("fisher_test")
}
#' Fisher's Exact Test for gwas count data
#'
#' @param obj object of class gwas
#'
#' @param reference reference group of subjects in which we want to perform fisher test
#' test
#' @param snp_start select starting position of snps.
#' @param snp_end select ending position of snp.
#' @param alternative indicates the alternative hypothesis and must be one of "two.sided", "greater" or "less".
#' You can specify just the initial letter. Only used in the 2 by 2 case
#' @param ... optional arguments to fisher.test
#'
#' @return the output will be a test_snps (data.table) object including 3 columns:
#' \dQuote{snp_pos}, \dQuote{case_control}, and \dQuote{value} which is
#' a p-values.
#'
#'@examples
#'\donttest{
#'piggyback::pb_download(repo = "soroushmdg/gwid",tag = "v0.0.1",dest = tempdir())
#'ibd_data_file <- paste0(tempdir(),"//chr3.ibd")
#'genome_data_file <- paste0(tempdir(),"//chr3.gds")
#'phase_data_file <- paste0(tempdir(),"//chr3.vcf")
#'case_control_data_file <- paste0(tempdir(),"//case-cont-RA.withmap.Rda")
#'# case-control data
#'case_control <- gwid::case_control(case_control_rda = case_control_data_file)
#'names(case_control) #cases and controls group
#'summary(case_control) # in here, we only consider cases,cont1,cont2,cont3 groups in the study
#'case_control$cases[1:3] # first three subject names of cases group
#'# read SNP data (use SNPRelate to convert it to gds) and count number of minor alleles
#'snp_data_gds <- gwid::build_gwas(gds_data = genome_data_file,
#'caco = case_control,gwas_generator = TRUE)
#'class(snp_data_gds)
#'names(snp_data_gds)
#'head(snp_data_gds$snps) # it has information about counts of minor alleles in each location.
#'# read haplotype data (output of beagle)
#'haplotype_data <- gwid::build_phase(phased_vcf = phase_data_file,caco = case_control)
#'class(haplotype_data)
#'names(haplotype_data)
#'dim(haplotype_data$Hap.1) #22302 SNP and 1911 subjects
#'# read IBD data (output of Refined-IBD)
#'ibd_data <- gwid::build_gwid(ibd_data = ibd_data_file,gwas = snp_data_gds)
#'class(ibd_data)
#'ibd_data$ibd # refined IBD output
#'ibd_data$res # count number of IBD for each SNP location
#'# plot count of IBD in chromosome 3
#'plot(ibd_data,y = c("cases","cont1"),ly = FALSE)
#'# Further investigate location between 117M and 122M
#'# significant number of IBD's in group cases, compare to cont1, cont2 and cont3.
#'plot(ibd_data,y = c("cases","cont1"),snp_start = 119026294,snp_end = 120613594,ly = FALSE)
#'model_fisher <- gwid::fisher_test(ibd_data,case_control,reference = "cases",
#'snp_start = 119026294,snp_end = 120613594)
#'class(model_fisher)
#'plot(model_fisher, y = c("cases","cont1"),ly = FALSE)
#'hap_str <- gwid::haplotype_structure(ibd_data,phase = haplotype_data,w = 10,
#'snp_start = 119026294,snp_end = 120613594)
#'haplo_freq <- gwid::haplotype_frequency(hap_str)
#'plot(haplo_freq,y = c("cases", "cont1"),plot_type = "haplotype_structure_frequency",
#'nwin = 1, type = "version1",ly = FALSE)
#'}
#' @export
fisher_test.gwas <- function(obj, reference, snp_start , snp_end, alternative = c("two.sided","greater","less"), ...) {
if (missing(obj)){
stop("please provide gwas object (output of function build_gwas)")
}
if (missing(reference)){
stop("please provide reference names: e.x. 'case1' " )
}
alternative <- match.arg(alternative)
if (sum(unlist(lapply(obj, inherits,"result_snps")))==0){
obj$snps <- extract(obj)
}
snp_pos <- case_control <- NULL
list_ind_result_snps <- which(unlist(lapply(obj, inherits,"result_snps")))
snp_pos_total <- unique(obj[[list_ind_result_snps]][["snp_pos"]])
if (missing(snp_start)) {
snp_start <- snp_pos_total[1]
}
if (missing(snp_end)) {
snp_end <- snp_pos_total[length(snp_pos_total)]
}
snp_indx <- which(snp_pos_total >= snp_start & snp_pos_total <= snp_end)
leni <- diff(range(snp_indx)) + 1
snp_pos_plot <- snp_pos_total[snp_indx][1:leni]
snps_mat <- data.table::dcast(obj[[list_ind_result_snps]][snp_pos %in% snp_pos_plot,], snp_pos ~ case_control , value.var = "value")[,-1,with=FALSE]
snps_mat <- as.matrix(snps_mat)
nas <- extract(obj,type = "nas")
nas_mat <- data.table::dcast(nas[snp_pos %in% snp_pos_plot,], snp_pos ~ case_control , value.var = "value")[,-1,with=FALSE]
nas_mat <- as.matrix(nas_mat)
idx <- which(colnames(snps_mat) == reference)
cc <- array(0, dim = c(2, nrow(snps_mat), ncol(snps_mat)))
pval <- matrix(0, nrow = nrow(snps_mat), ncol = ncol(snps_mat))
cc[1, , ] <- snps_mat # fill the first matrix with snps_mat
nums <- as.integer(summary(obj[["caco"]])[1:length(obj[["caco"]])])
cc[2, , ] <- rep(nums * 2, each = nrow(snps_mat)) - snps_mat - nas_mat * 2
cc[2,,][cc[2,,] <0] <- 0
for (i in 1:nrow(snps_mat)) {
for (j in (1:ncol(snps_mat))[-idx]) pval[i, j] <- stats::fisher.test(cc[, i, c(idx, j)], alternative = alternative, ...)$p.value
}
pval[,idx] <- 1
colnames(pval) <- colnames(snps_mat)
pval_table <- data.table::as.data.table(pval)[,snp_pos:=snp_pos_total[snp_indx]]
pval_table <- data.table::melt(pval_table,id.vars="snp_pos",value.name = "value",
measure.vars=colnames(pval),variable.name="case_control")
output <- pval_table
output <- output[case_control!=reference]
class(output) <- append("test_snps",class(output))
return(output)
}
#' Fisher's Exact Test for gwid count data
#'
#' @param obj An object of class gwid. Output of \code{build_gwid} function
#'
#' @param caco An object of class caco. Output of \code{case_control} function.
#' @param snp_start select starting position of snps.
#' @param snp_end select ending position of snp.
#' @param reference reference group of subjects in which we want to perform fisher test
#' @param alternative indicates the alternative hypothesis and must be one of "two.sided", "greater" or "less".
#' You can specify just the initial letter. Only used in the 2 by 2 case
#' @param ... optional arguments to fisher.test
#'
#' @return the output will be a test_snps (data.table) object including 3 columns:
#' \dQuote{snp_pos}, \dQuote{case_control}, and \dQuote{value} which is
#' a p-values.
#'
#'@examples
#'\donttest{
#'piggyback::pb_download(repo = "soroushmdg/gwid",tag = "v0.0.1",dest = tempdir())
#'ibd_data_file <- paste0(tempdir(),"//chr3.ibd")
#'genome_data_file <- paste0(tempdir(),"//chr3.gds")
#'phase_data_file <- paste0(tempdir(),"//chr3.vcf")
#'case_control_data_file <- paste0(tempdir(),"//case-cont-RA.withmap.Rda")
#'# case-control data
#'case_control <- gwid::case_control(case_control_rda = case_control_data_file)
#'names(case_control) #cases and controls group
#'summary(case_control) # in here, we only consider cases,cont1,cont2,cont3 groups in the study
#'case_control$cases[1:3] # first three subject names of cases group
#'# read SNP data (use SNPRelate to convert it to gds) and count number of minor alleles
#'snp_data_gds <- gwid::build_gwas(gds_data = genome_data_file,
#'caco = case_control,gwas_generator = TRUE)
#'class(snp_data_gds)
#'names(snp_data_gds)
#'head(snp_data_gds$snps) # it has information about counts of minor alleles in each location.
#'# read haplotype data (output of beagle)
#'haplotype_data <- gwid::build_phase(phased_vcf = phase_data_file,caco = case_control)
#'class(haplotype_data)
#'names(haplotype_data)
#'dim(haplotype_data$Hap.1) #22302 SNP and 1911 subjects
#'# read IBD data (output of Refined-IBD)
#'ibd_data <- gwid::build_gwid(ibd_data = ibd_data_file,gwas = snp_data_gds)
#'class(ibd_data)
#'ibd_data$ibd # refined IBD output
#'ibd_data$res # count number of IBD for each SNP location
#'# plot count of IBD in chromosome 3
#'plot(ibd_data,y = c("cases","cont1"),ly = FALSE)
#'# Further investigate location between 117M and 122M
#'# significant number of IBD's in group cases, compare to cont1, cont2 and cont3.
#'plot(ibd_data,y = c("cases","cont1"),snp_start = 119026294,snp_end = 120613594,ly = FALSE)
#'model_fisher <- gwid::fisher_test(ibd_data,case_control,reference = "cases",
#'snp_start = 119026294,snp_end = 120613594)
#'class(model_fisher)
#'plot(model_fisher, y = c("cases","cont1"),ly = FALSE)
#'hap_str <- gwid::haplotype_structure(ibd_data,phase = haplotype_data,w = 10,
#'snp_start = 119026294,snp_end = 120613594)
#'haplo_freq <- gwid::haplotype_frequency(hap_str)
#'plot(haplo_freq,y = c("cases", "cont1"),plot_type = "haplotype_structure_frequency",
#'nwin = 1, type = "version1",ly = FALSE)
#'}
#' @export
fisher_test.gwid <- function(obj, caco,snp_start ,snp_end ,reference , alternative = c("two.sided","greater","less"), ...) {
if (missing(obj)){
stop("please provide gwid object (output of function build_gwid)")
}
if (missing(caco)){
stop("please provide caco object(case_control)")
}
if (missing(reference)){
stop("please provide a reference e.x. 'case1' ")
}
snp_pos <- NULL
alternative <- match.arg(alternative)
#type <- match.arg(type)
if (sum(unlist(lapply(obj, inherits,"result_snps")))==0){
obj$res <- obj(obj)
}
list_ind_result_snps <- which(unlist(lapply(obj, inherits,"result_snps")))
snp_pos_total <- unique(obj[[list_ind_result_snps]]$snp_pos)
if (missing(snp_start)) {
snp_start <- snp_pos_total[1]
}
if (missing(snp_end)) {
snp_end <- snp_pos_total[length(snp_pos_total)]
}
snp_indx <- which(snp_pos_total >= snp_start & snp_pos_total <= snp_end)
leni <- diff(range(snp_indx)) + 1
snp_pos_plot <- snp_pos_total[snp_indx][1:leni]
temp <- obj[[list_ind_result_snps]][snp_pos %in% snp_pos_plot,]
output <- fisher_test(temp,caco=caco,reference=reference,alternative=alternative, ...)
return(output)
}
#' permutation test
#'
#' @param obj object
#'
#' @param ... other variables
#'
#' @return the output will be a test_snps (data.table) object including 3 columns:
#' \dQuote{snp_pos}, \dQuote{case_control}, and \dQuote{value} which is
#' a p-values.
#' @export
permutation_test <- function(obj, ...) {
UseMethod("permutation_test",obj)
}
#' permutation test for gwid count data
#'
#' @param obj An object of class gwid. Output of \code{build_gwid} function
#'
#' @param gwas object of class gwas
#' @param snp_start select starting position of snps.
#' @param snp_end select ending position of snp.
#' @param nperm Number of permutations.
#' @param reference reference group
#' @param ... other variables
#' @return the output will be a test_snps (data.table) object including 3 columns:
#' \dQuote{snp_pos}, \dQuote{case_control}, and \dQuote{value} which is
#' a p-values.
#'@examples
#'\donttest{
#'piggyback::pb_download(repo = "soroushmdg/gwid",tag = "v0.0.1",dest = tempdir())
#'ibd_data_file <- paste0(tempdir(),"//chr3.ibd")
#'genome_data_file <- paste0(tempdir(),"//chr3.gds")
#'phase_data_file <- paste0(tempdir(),"//chr3.vcf")
#'case_control_data_file <- paste0(tempdir(),"//case-cont-RA.withmap.Rda")
#'# case-control data
#'case_control <- gwid::case_control(case_control_rda = case_control_data_file)
#'names(case_control) #cases and controls group
#'summary(case_control) # in here, we only consider cases,cont1,cont2,cont3 groups in the study
#'case_control$cases[1:3] # first three subject names of cases group
#'# read SNP data (use SNPRelate to convert it to gds) and count number of minor alleles
#'snp_data_gds <- gwid::build_gwas(gds_data = genome_data_file,
#'caco = case_control,gwas_generator = TRUE)
#'class(snp_data_gds)
#'names(snp_data_gds)
#'head(snp_data_gds$snps) # it has information about counts of minor alleles in each location.
#'# read haplotype data (output of beagle)
#'haplotype_data <- gwid::build_phase(phased_vcf = phase_data_file,caco = case_control)
#'class(haplotype_data)
#'names(haplotype_data)
#'dim(haplotype_data$Hap.1) #22302 SNP and 1911 subjects
#'# read IBD data (output of Refined-IBD)
#'ibd_data <- gwid::build_gwid(ibd_data = ibd_data_file,gwas = snp_data_gds)
#'class(ibd_data)
#'ibd_data$ibd # refined IBD output
#'ibd_data$res # count number of IBD for each SNP location
#'# plot count of IBD in chromosome 3
#'plot(ibd_data,y = c("cases","cont1"),ly = FALSE)
#'# Further investigate location between 117M and 122M
#'# significant number of IBD's in group cases, compare to cont1, cont2 and cont3.
#'plot(ibd_data,y = c("cases","cont1"),snp_start = 119026294,snp_end = 120613594,ly = FALSE)
#'model_fisher <- gwid::fisher_test(ibd_data,case_control,reference = "cases",
#'snp_start = 119026294,snp_end = 120613594)
#'model_permutation <- permutation_test(ibd_data,snp_data_gds,
#'snp_start = 119026294,snp_end = 120613594,nperm=20,reference = "cases")
#'class(model_fisher)
#'plot(model_fisher, y = c("cases","cont1"),ly = FALSE)
#'hap_str <- gwid::haplotype_structure(ibd_data,phase = haplotype_data,w = 10,
#'snp_start = 119026294,snp_end = 120613594)
#'haplo_freq <- gwid::haplotype_frequency(hap_str)
#'plot(haplo_freq,y = c("cases", "cont1"),plot_type = "haplotype_structure_frequency",
#'nwin = 1, type = "version1",ly = FALSE)
#'}
#' @export
permutation_test.gwid <- function(obj, gwas, snp_start, snp_end,
nperm = 100, reference = "cases", ...) {
if (missing(obj)) {
stop("please provide gwid object (output of function build_gwid)")
}
if (missing(gwas)) {
stop("please provide gwas object")
}
if (missing(reference)) {
stop("please provide reference group")
}
if (sum(unlist(lapply(obj, inherits, "result_snps"))) == 0) {
obj$res <- obj(obj)
}
if (missing(snp_start)) {
snp_start <- snp_pos_total[1]
}
if (missing(snp_end)) {
snp_end <- snp_pos_total[length(snp_pos_total)]
}
list_ind_result_snps <- which(unlist(lapply(obj, inherits, "result_snps")))
value <- snp_pos <- case_control <- V1 <- V3 <- NULL
snp_pos_total <- unique(obj[[list_ind_result_snps]]$snp_pos)
snp_indx <- which(snp_pos_total >= snp_start & snp_pos_total <= snp_end)
leni <- diff(range(snp_indx)) + 1
snp_pos_plot <- snp_pos_total[snp_indx][1:leni]
temp <- obj[[list_ind_result_snps]][snp_pos %in% snp_pos_plot, ]
myres <- as.matrix(data.table::dcast(obj$res[snp_pos %in% snp_pos_plot], snp_pos ~ case_control, value.var = "value")[, -1, with = FALSE])
myres_gap <- myres[,reference] - myres
perm_matrix <- matrix(0,nrow = nrow(myres_gap),ncol = ncol(myres_gap))
gwas_temp <- gwas
gwas_temp$caco <- list(all_subj = unique(unlist(gwas_temp$caco, use.names = FALSE)))
class(gwas_temp$caco) <- class(gwas[["caco"]])
build_gwid_modify <- function(gwid = "object of class gwid", gwas = "object of class gwas", gwid_generator = TRUE) {
ibd <- gwid$ibd
V1 <- V2 <- V3 <- V4 <- V5 <- V6 <- V7 <- V8 <- V9 <- NULL
ibd <- ibd[V1 %in% unlist(unique(gwas[["caco"]])) & V3 %in% unlist(unique(gwas[["caco"]]))]
class(ibd) <- append("IBD", class(ibd))
seq2 <- Vectorize(seq.default, vectorize.args = c("from", "to"))
profile <- ind <- vector(mode = "list", length = length(gwas[["caco"]])) # list length 6
for (j in seq_along(gwas[["caco"]])) {
ind[[j]] <- which(ibd$V1 %in% gwas[["caco"]][[j]] & ibd$V3 %in% gwas[["caco"]][[j]])
a1 <- ibd[V1 %in% gwas[["caco"]][[j]] & V3 %in% gwas[["caco"]][[j]]]
a2 <- seq2(match(a1$V6, gwas$snp.pos), match(a1$V7, gwas$snp.pos))
Un1 <- unlist(a2)
profile[[j]] <- Matrix::sparseMatrix(
i = rep(seq_along(a2), lengths(a2)),
j = Un1,
x = 1
)
if (ncol(profile[[j]]) < length(gwas$snp.pos)) {
if (min(Un1) > 1) {
mytemp <- Matrix(0, nrow = nrow(profile[[j]]), ncol = (min(Un1) - 1))
profile[[j]] <- cbind(mytemp, profile[[j]])
}
if (max(Un1) < length(gwas$snp.pos)) {
mytemp <- Matrix(0, nrow = nrow(profile[[j]]), ncol = (length(gwas$snp.pos) - max(Un1)))
profile[[j]] <- cbind(profile[[j]], mytemp)
}
profile[[j]] <- methods::as(profile[[j]], "sparseMatrix")
}
}
names(profile) <- names(ind) <- names(gwas[["caco"]])
class(profile) <- "profile"
output <- list(profile = profile, IND = ind)
output$snp_pos <- gwas[["snp.pos"]]
output$ibd <- ibd
class(output) <- "gwid"
if (gwid_generator) {
output$res <- extract(output)
}
return(output)
}
myregion2_temp <- build_gwid_modify(gwid = obj, gwas = gwas_temp)
mres <- myregion2_temp$profile$all_subj
ibd <- myregion2_temp$ibd
subj.ind <- list()
subj.ind[[1]] <- (ibd[,V1])
subj.ind[[2]] <- (ibd[,V3])
caco_temp <- gwas[["caco"]]
mylength <- length(caco_temp[[reference]])
caco_temp_unlist <- unique(unlist(caco_temp, use.names = FALSE))
for (i in 1:nperm) {
cacoi <- matrix(sample(caco_temp_unlist,2*mylength),ncol=2)
indi <- list()
mresi <- list()
indi[[1]] <- which(subj.ind[[1]]%in%cacoi[,1] & subj.ind[[2]]%in%cacoi[,1])
indi[[2]] <- which(subj.ind[[1]]%in%cacoi[,2] & subj.ind[[2]]%in%cacoi[,2])
mresi[[1]] <- colSums(mres[indi[[1]],snp_indx])
mresi[[2]] <- colSums(mres[indi[[2]],snp_indx])
perm_matrix <- perm_matrix + 1*( myres_gap - matrix(mresi[[1]]-mresi[[2]],length(mresi[[1]]-mresi[[2]]),ncol(myres_gap))<= 0 )
}
perm_matrix <- cbind(snp_pos_plot,perm_matrix/nperm)
colnames(perm_matrix) <- c("snp_pos",colnames(perm_matrix)[-1])
perm_table <- data.table::as.data.table(perm_matrix)
pval <- data.table::melt(perm_table, id.vars = "snp_pos", variable.name = "case_control", value.name = "value")
pval[value==0,value:=1/(2*nperm)]
pval <- pval[case_control!=reference]
class(pval) <- append("test_snps", class(pval))
return(pval)
}
#' Permutation test for gwas object
#'
#' @param obj object of class gwas
#'
#' @param snp_start elect starting position of snps.
#' @param snp_end select ending position of snp.
#' @param nperm Number of permutations.
#' @param reference reference group of subjects in which we want to perform fisher test
#' @param ... other variables
#' @return the output will be a test_snps (data.table) object including 3 columns:
#' \dQuote{snp_pos}, \dQuote{case_control}, and \dQuote{value} which is
#' a p-values.
#'@examples
#'\donttest{
#'piggyback::pb_download(repo = "soroushmdg/gwid",tag = "v0.0.1",dest = tempdir())
#'ibd_data_file <- paste0(tempdir(),"//chr3.ibd")
#'genome_data_file <- paste0(tempdir(),"//chr3.gds")
#'phase_data_file <- paste0(tempdir(),"//chr3.vcf")
#'case_control_data_file <- paste0(tempdir(),"//case-cont-RA.withmap.Rda")
#'# case-control data
#'case_control <- gwid::case_control(case_control_rda = case_control_data_file)
#'names(case_control) #cases and controls group
#'summary(case_control) # in here, we only consider cases,cont1,cont2,cont3 groups in the study
#'case_control$cases[1:3] # first three subject names of cases group
#'# read SNP data (use SNPRelate to convert it to gds) and count number of minor alleles
#'snp_data_gds <- gwid::build_gwas(gds_data = genome_data_file,
#'caco = case_control,gwas_generator = TRUE)
#'class(snp_data_gds)
#'names(snp_data_gds)
#'head(snp_data_gds$snps) # it has information about counts of minor alleles in each location.
#'# read haplotype data (output of beagle)
#'haplotype_data <- gwid::build_phase(phased_vcf = phase_data_file,caco = case_control)
#'class(haplotype_data)
#'names(haplotype_data)
#'dim(haplotype_data$Hap.1) #22302 SNP and 1911 subjects
#'# read IBD data (output of Refined-IBD)
#'ibd_data <- gwid::build_gwid(ibd_data = ibd_data_file,gwas = snp_data_gds)
#'class(ibd_data)
#'ibd_data$ibd # refined IBD output
#'ibd_data$res # count number of IBD for each SNP location
#'# plot count of IBD in chromosome 3
#'plot(ibd_data,y = c("cases","cont1"),ly = FALSE)
#'# Further investigate location between 117M and 122M
#'# significant number of IBD's in group cases, compare to cont1, cont2 and cont3.
#'plot(ibd_data,y = c("cases","cont1"),snp_start = 119026294,snp_end = 120613594,ly = FALSE)
#'model_fisher <- gwid::fisher_test(ibd_data,case_control,reference = "cases",
#'snp_start = 119026294,snp_end = 120613594)
#'model_permutation <- permutation_test(ibd_data,snp_data_gds,
#'snp_start = 119026294,snp_end = 120613594,nperm=20,reference = "cases")
#'class(model_fisher)
#'plot(model_fisher, y = c("cases","cont1"),ly = FALSE)
#'hap_str <- gwid::haplotype_structure(ibd_data,phase = haplotype_data,w = 10,
#'snp_start = 119026294,snp_end = 120613594)
#'haplo_freq <- gwid::haplotype_frequency(hap_str)
#'plot(haplo_freq,y = c("cases", "cont1"),plot_type = "haplotype_structure_frequency",
#'nwin = 1, type = "version1",ly = FALSE)
#'}
#' @export
permutation_test.gwas <- function(obj, snp_start, snp_end,
nperm = 1000, reference = "cases", ...) {
if (missing(obj)) {
stop("please provide gwas object")
}
if (missing(snp_start)) {
snp_start <- snp_pos_total[1]
}
if (missing(snp_end)) {
snp_end <- snp_pos_total[length(snp_pos_total)]
}
value <- snp_pos <- case_control <- NULL
list_ind_result_snps <- which(unlist(lapply(obj, inherits, "result_snps")))
snp_pos_total <- unique(obj[[list_ind_result_snps]]$snp_pos)
snp_indx <- which(snp_pos_total >= snp_start & snp_pos_total <= snp_end)
leni <- diff(range(snp_indx)) + 1
snp_pos_plot <- snp_pos_total[snp_indx][1:leni]
caco <- obj$caco
total_subjects <- sapply(obj$caco, length)
myres <- as.matrix(data.table::dcast(obj$snps[snp_pos %in% snp_pos_plot], snp_pos ~ case_control, value.var = "value")[, -1, with = FALSE])
myres_gap <- myres[, reference] - myres
perm_matrix <- matrix(0, nrow = nrow(myres_gap), ncol = ncol(myres_gap))
gwas_temp <- obj
gwas_temp$caco <- list(all_subj = (unlist(gwas_temp$caco, use.names = FALSE)))
class(gwas_temp$caco) <- class(obj$caco)
mres <- do.call(what = "rbind", args = obj$smp.snp)
gwas_temp$smp.snp <- list(all_subj = mres)
for (i in 1:nperm) {
caco_temp <- obj$caco
cacoi <- split(sample(gwas_temp$caco$all_subj), rep(1:length(obj$caco), total_subjects))
names(cacoi) <- names(obj$caco)
new_snps <- sapply(cacoi, function(x) {
colSums(gwas_temp$smp.snp$all_subj[rownames(gwas_temp$smp.snp$all_subj) %in% x, snp_indx])
})
myres_gap_perm <- new_snps[, reference] - new_snps
perm_matrix <- perm_matrix + 1 * (myres_gap - myres_gap_perm <= 0)
}
perm_matrix <- cbind(snp_pos_plot, perm_matrix / nperm)
colnames(perm_matrix) <- c("snp_pos", colnames(perm_matrix)[-1])
perm_table <- data.table::as.data.table(perm_matrix)
pval <- data.table::melt(perm_table, id.vars = "snp_pos", variable.name = "case_control", value.name = "value")
pval[value==0,value:=1/(2*nperm)]
pval <- pval[case_control!=reference]
class(pval) <- append("test_snps",class(pval))
return(pval)
}
G_Test <- function (x, y = NULL, correct = c("none", "williams", "yates"),
p = rep(1/length(x), length(x)))
{
DNAME <- deparse(substitute(x))
if (is.data.frame(x))
x <- as.matrix(x)
if (is.matrix(x)) {
if (min(dim(x)) == 1)
x <- as.vector(x)
}
if (!is.matrix(x) && !is.null(y)) {
if (length(x) != length(y))
stop("x and y must have the same length")
DNAME <- paste(DNAME, "and", deparse(substitute(y)))
OK <- stats::complete.cases(x, y)
x <- as.factor(x[OK])
y <- as.factor(y[OK])
if ((nlevels(x) < 2) || (nlevels(y) < 2))
stop("x and y must have at least 2 levels")
x <- table(x, y)
}
if (any(x < 0) || any(is.na(x)))
stop("all entries of x must be nonnegative and finite")
if ((n <- sum(x)) == 0)
stop("at least one entry of x must be positive")
correct <- match.arg(correct)
if (is.matrix(x)) {
nrows <- nrow(x)
ncols <- ncol(x)
if (correct == "yates") {
if (dim(x)[1] != 2 || dim(x)[2] != 2)
stop("Yates' correction requires a 2 x 2 matrix")
if ((x[1, 1] * x[2, 2]) - (x[1, 2] * x[2, 1]) > 0) {
x <- x + 0.5
diag(x) <- diag(x) - 1
}
else {
x <- x - 0.5
diag(x) <- diag(x) + 1
}
}
sr <- apply(x, 1, sum)
sc <- apply(x, 2, sum)
E <- outer(sr, sc, "*")/n
g <- 0
for (i in 1:nrows) {
for (j in 1:ncols) {
if (x[i, j] != 0)
g <- g + x[i, j] * log(x[i, j]/E[i, j])
}
q <- 1
if (correct == "williams") {
row.tot <- col.tot <- 0
for (i in 1:nrows) {
row.tot <- row.tot + 1/(sum(x[i, ]))
}
for (j in 1:ncols) {
col.tot <- col.tot + 1/(sum(x[, j]))
}
q <- 1 + ((n * row.tot - 1) * (n * col.tot -
1))/(6 * n * (ncols - 1) * (nrows - 1))
}
STATISTIC <- G <- 2 * g/q
PARAMETER <- (nrow(x) - 1) * (ncol(x) - 1)
PVAL <- 1 - stats::pchisq(STATISTIC, df = PARAMETER)
if (correct == "none")
METHOD <- "Log likelihood ratio (G-test) test of independence without correction"
if (correct == "williams")
METHOD <- "Log likelihood ratio (G-test) test of independence with Williams' correction"
if (correct == "yates")
METHOD <- "Log likelihood ratio (G-test) test of independence with Yates' correction"
}
}
else {
METHOD <- "Log likelihood ratio (G-test) goodness of fit test"
if (length(x) == 1)
stop("x must at least have 2 elements")
if (length(x) != length(p))
stop("x and p must have the same number of elements")
E <- n * p
if (correct == "yates") {
if (length(x) != 2)
stop("Yates' correction requires 2 data values")
if ((x[1] - E[1]) > 0.25) {
x[1] <- x[1] - 0.5
x[2] <- x[2] + 0.5
}
else if ((E[1] - x[1]) > 0.25) {
x[1] <- x[1] + 0.5
x[2] <- x[2] - 0.5
}
}
names(E) <- names(x)
g <- 0
for (i in 1:length(x)) {
if (x[i] != 0)
g <- g + x[i] * log(x[i]/E[i])
}
q <- 1
if (correct == "williams") {
q <- 1 + (length(x) + 1)/(6 * n)
}
STATISTIC <- G <- 2 * g/q
PARAMETER <- length(x) - 1
PVAL <- stats::pchisq(STATISTIC, PARAMETER, lower.tail = FALSE)
}
names(STATISTIC) <- "G"
names(PARAMETER) <- "X-squared df"
names(PVAL) <- "p.value"
structure(list(statistic = STATISTIC, parameter = PARAMETER,
p.value = PVAL, method = METHOD, data.name = DNAME, observed = x,
expected = E), class = "htest")
}
#' perform gtest
#'
#' @param haplotype_structure object of a class
#'
#' @param ... other variables
#' @return the output will be a test_snps (data.table) object including 3 columns:
#' \dQuote{snp_pos}, \dQuote{case_control}, and \dQuote{value} which is
#' a p-values.
#' @export
gtest <- function(haplotype_structure, ...) {
UseMethod("gtest")
}
#' Perform G-test on haplotype structures extracted from \code{haplotype_structure} function
#'
#' @param haplotype_structure An object of class haplotype_structure. Output of
#' \code{haplotype_structure} function.
#'
#' @param reference reference group of subjects in which we want to perform G-test
#'
#' @param ... other variables
#'
#'@return the output will be a test_snps (data.table) object including 3 columns:
#' \dQuote{snp_pos}, \dQuote{case_control}, and \dQuote{value} which is
#' a p-values.
#'
#' @export
gtest.haplotype_structure <- function(haplotype_structure, reference, ...) {
if (missing(haplotype_structure) | missing(reference)){
stop("please provide function arguments")
}
pval <- matrix(0, nrow = length(unique(haplotype_structure$snp_pos)), ncol = length(unique(haplotype_structure$case_control)))
colnames(pval) <- levels(haplotype_structure$case_control)
idx <- which(colnames(pval) == reference)
for (j in (1:ncol(pval))[-idx]) {
caco <- c(reference, colnames(pval)[j])
case_control <- snp_pos <- structures <- window_number <- NULL
tmp <- haplotype_structure[case_control %in% caco, list(case_control, structures, window_number)]
tmp <- split(tmp, tmp[, window_number])
tmp <- lapply(tmp, function(x) {
y <- x[case_control %in% caco, !"window_number"]
y <- table(y)[caco, ]
})
pval[, j] <- unlist(lapply(tmp, function(x) {
y <- G_Test(x, correct = "williams")$p.value
y
}))
}
pval[, idx] <- 1
pval_table <- data.table::as.data.table(pval)[, snp_pos := sort(unique(haplotype_structure$snp_pos))]
pval_table <- data.table::melt(pval_table,
id.vars = "snp_pos", value.name = "value",
measure.vars = colnames(pval), variable.name = "case_control"
)
output <- data.table::copy(pval_table)
output <- output[case_control!=reference]
class(output) <- append("test_snps",class(output))
return(output)
}
#' fisher exact test for result_snps count data
#'
#' @param obj An object of class result_snps
#'
#' @param caco An object of class caco. Output of \code{case_control} function.
#' @param reference reference group of subjects in which we want to perform fisher test.
#' @param alternative indicates the alternative hypothesis and must be one of "two.sided", "greater" or "less".
#' You can specify just the initial letter. Only used in the 2 by 2 case
#' @param ... optional arguments to fisher.test
#' @return the output will be a test_snps (data.table) object including 3 columns:
#' \dQuote{snp_pos}, \dQuote{case_control}, and \dQuote{value} which is
#' a p-values.
#'
#'@examples
#'\donttest{
#'piggyback::pb_download(repo = "soroushmdg/gwid",tag = "v0.0.1",dest = tempdir())
#'ibd_data_file <- paste0(tempdir(),"//chr3.ibd")
#'genome_data_file <- paste0(tempdir(),"//chr3.gds")
#'phase_data_file <- paste0(tempdir(),"//chr3.vcf")
#'case_control_data_file <- paste0(tempdir(),"//case-cont-RA.withmap.Rda")
#'# case-control data
#'case_control <- gwid::case_control(case_control_rda = case_control_data_file)
#'names(case_control) #cases and controls group
#'summary(case_control) # in here, we only consider cases,cont1,cont2,cont3 groups in the study
#'case_control$cases[1:3] # first three subject names of cases group
#'# read SNP data (use SNPRelate to convert it to gds) and count number of minor alleles
#'snp_data_gds <- gwid::build_gwas(gds_data = genome_data_file,
#'caco = case_control,gwas_generator = TRUE)
#'class(snp_data_gds)
#'names(snp_data_gds)
#'head(snp_data_gds$snps) # it has information about counts of minor alleles in each location.
#'# read haplotype data (output of beagle)
#'haplotype_data <- gwid::build_phase(phased_vcf = phase_data_file,caco = case_control)
#'class(haplotype_data)
#'names(haplotype_data)
#'dim(haplotype_data$Hap.1) #22302 SNP and 1911 subjects
#'# read IBD data (output of Refined-IBD)
#'ibd_data <- gwid::build_gwid(ibd_data = ibd_data_file,gwas = snp_data_gds)
#'class(ibd_data)
#'ibd_data$ibd # refined IBD output
#'ibd_data$res # count number of IBD for each SNP location
#'# plot count of IBD in chromosome 3
#'plot(ibd_data,y = c("cases","cont1"),ly = FALSE)
#'# Further investigate location between 117M and 122M
#'# significant number of IBD's in group cases, compare to cont1, cont2 and cont3.
#'plot(ibd_data,y = c("cases","cont1"),snp_start = 119026294,snp_end = 120613594,ly = FALSE)
#'model_fisher <- gwid::fisher_test(ibd_data,case_control,reference = "cases",
#'snp_start = 119026294,snp_end = 120613594)
#'model_permutation <- permutation_test(ibd_data,snp_data_gds,
#'snp_start = 119026294,snp_end = 120613594,nperm=20,reference = "cases")
#'class(model_fisher)
#'plot(model_fisher, y = c("cases","cont1"),ly = FALSE)
#'hap_str <- gwid::haplotype_structure(ibd_data,phase = haplotype_data,w = 10,
#'snp_start = 119026294,snp_end = 120613594)
#'haplo_freq <- gwid::haplotype_frequency(hap_str)
#'plot(haplo_freq,y = c("cases", "cont1"),plot_type = "haplotype_structure_frequency",
#'nwin = 1, type = "version1",ly = FALSE)
#'}
#' @export
fisher_test.result_snps <- function(obj, caco, reference, alternative = c("two.sided","greater","less"), ...) {
if (missing(obj)){
stop("Please provide an object of class result_snps")
}
if (missing(caco)){
stop("Please provide an object of class caco")
}
if (missing(reference)){
stop("Please provide reference e.x. 'case1' ")
}
alternative <- match.arg(alternative)
res_mat <- data.table::dcast(obj, snp_pos ~ case_control , value.var = "value")
snp_pos <- res_mat[["snp_pos"]]
res_mat <- as.matrix(res_mat[,-1,with=FALSE])
#idx is column index of reference
idx <- which(colnames(res_mat) == reference)
cc <- array(0, dim = c(2, nrow(res_mat), ncol(res_mat)))
pval <- matrix(0, nrow = nrow(res_mat), ncol = ncol(res_mat))
cc[1, , ] <- res_mat
nums <- as.integer(summary(caco)[1:length(caco)]) # number of samples in each case_control
cc[2, , ] <- rep(nums * (nums - 1) / 2, each = nrow(res_mat)) - res_mat
for (i in 1:nrow(res_mat)) {
for (j in (1:ncol(res_mat))[-idx]) pval[i, j] <- fisher.test(cc[, i, c(idx, j)], alternative = alternative)$p.value
}
colnames(pval) <- colnames(res_mat)
pval[,idx] <- 1
pval_table <- data.table::as.data.table(pval)[,snp_pos:=snp_pos]
# transform pval to long format 3 columns (snp_pos, case_control, value)
pval_table <- data.table::melt(pval_table,id.vars="snp_pos",value.name = "value",
measure.vars=colnames(pval),variable.name="case_control")
pval_table <- pval_table[case_control!=reference]
class(pval_table) <- append("test_snps",class(pval_table))
return(pval_table)
}
#' Permutation test for `haplotype_structure` object
#'
#' @param obj object of class `haplotype_structure`
#'
#' @param nperm Number of permutations.
#' @param reference reference group of subjects in which we want to perform `gtest`
#' @param ... other variables
#' @return the output will be a test_snps (data.table) object including 3 columns:
#' \dQuote{snp_pos}, \dQuote{case_control}, and \dQuote{value} which is
#' a p-values.
#'@examples
#'\donttest{
#'piggyback::pb_download(repo = "soroushmdg/gwid",tag = "v0.0.1",dest = tempdir())
#'ibd_data_file <- paste0(tempdir(),"//chr3.ibd")
#'genome_data_file <- paste0(tempdir(),"//chr3.gds")
#'phase_data_file <- paste0(tempdir(),"//chr3.vcf")
#'case_control_data_file <- paste0(tempdir(),"//case-cont-RA.withmap.Rda")
#'# case-control data
#'case_control <- gwid::case_control(case_control_rda = case_control_data_file)
#'names(case_control) #cases and controls group
#'summary(case_control) # in here, we only consider cases,cont1,cont2,cont3 groups in the study
#'case_control$cases[1:3] # first three subject names of cases group
#'# read SNP data (use SNPRelate to convert it to gds) and count number of minor alleles
#'snp_data_gds <- gwid::build_gwas(gds_data = genome_data_file,
#'caco = case_control,gwas_generator = TRUE)
#'class(snp_data_gds)
#'names(snp_data_gds)
#'head(snp_data_gds$snps) # it has information about counts of minor alleles in each location.
#'# read haplotype data (output of beagle)
#'haplotype_data <- gwid::build_phase(phased_vcf = phase_data_file,caco = case_control)
#'class(haplotype_data)
#'names(haplotype_data)
#'dim(haplotype_data$Hap.1) #22302 SNP and 1911 subjects
#'# read IBD data (output of Refined-IBD)
#'ibd_data <- gwid::build_gwid(ibd_data = ibd_data_file,gwas = snp_data_gds)
#'class(ibd_data)
#'ibd_data$ibd # refined IBD output
#'ibd_data$res # count number of IBD for each SNP location
#'# plot count of IBD in chromosome 3
#'plot(ibd_data,y = c("cases","cont1"),ly = FALSE)
#'# Further investigate location between 117M and 122M
#'# significant number of IBD's in group cases, compare to cont1, cont2 and cont3.
#'plot(ibd_data,y = c("cases","cont1"),snp_start = 119026294,snp_end = 120613594,ly = FALSE)
#'model_fisher <- gwid::fisher_test(ibd_data,case_control,reference = "cases",
#'snp_start = 119026294,snp_end = 120613594)
#'model_permutation <- permutation_test(ibd_data,snp_data_gds,
#'snp_start = 119026294,snp_end = 120613594,nperm=20,reference = "cases")
#'class(model_fisher)
#'plot(model_fisher, y = c("cases","cont1"),ly = FALSE)
#'hap_str <- gwid::haplotype_structure(ibd_data,phase = haplotype_data,w = 10,
#'snp_start = 119026294,snp_end = 120613594)
#'haplo_freq <- gwid::haplotype_frequency(hap_str)
#'plot(haplo_freq,y = c("cases", "cont1"),plot_type = "haplotype_structure_frequency",
#'nwin = 1, type = "version1",ly = FALSE)
#'}
#' @export
permutation_test.haplotype_structure <- function(obj, nperm, reference, ...) {
if (missing(obj) | missing(reference)) {
stop("please provide function arguments")
}
value <- NULL
statistics <- matrix(0, nrow = length(unique(obj$snp_pos)), ncol = length(unique(obj$case_control)))
colnames(statistics) <- levels(obj$case_control)
idx <- which(colnames(statistics) == reference)
statistics_perm <- vector(mode = "list", length = nperm)
for (i in 1:nperm) {
statistics_perm[[i]] <- matrix(0, nrow = nrow(statistics), ncol = ncol(statistics))
}
Freq <- NULL
for (j in (1:ncol(statistics))[-idx]) {
caco <- c(reference, colnames(statistics)[j])
case_control <- snp_pos <- structures <- window_number <- NULL
tmp <- obj[case_control %in% caco, list(case_control, structures, window_number)]
tmp <- split(tmp, tmp[, window_number])
tmp <- lapply(tmp, function(x) {
y <- x[case_control %in% caco, !"window_number"]
y <- table(y)[caco, ]
})
statistics[, j] <- unlist(lapply(tmp, function(x) {
y <- G_Test(x, correct = "williams")$statistic
y
}))
for (i in 1:nperm) {
statistics_perm[[i]][, j] <- unlist(lapply(tmp, function(x) {
data_frame_tmp <- as.data.frame(x)
data_frame_tmp[, 3] <- sample(data_frame_tmp[, 3])
cn <- colnames(data_frame_tmp)
shuffled_table_obj <- xtabs(Freq ~ ., data_frame_tmp)
y <- G_Test(shuffled_table_obj, correct = "williams")$statistic
y
}))
}
}
pval <- matrix(0, nrow = nrow(statistics), ncol = ncol(statistics))
pval <- Reduce("+", lapply(statistics_perm, function(x) {
x > statistics
})) / nperm
colnames(pval) <- colnames(statistics)
pval[, idx] <- 1
pval_table <- data.table::as.data.table(pval)[, snp_pos := sort(unique(obj$snp_pos))]
pval_table <- data.table::melt(pval_table,
id.vars = "snp_pos", value.name = "value",
measure.vars = colnames(statistics), variable.name = "case_control"
)
output <- data.table::copy(pval_table)
output[value==0,value:=1/(2*nperm)]
output <- output[case_control!=reference]
class(output) <- append("test_snps", class(output))
return(output)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.