R/Dynamic_Window_Results_Summary.R
In xihaoli/STAARpipelineSummary: Summarization and Visualization of Analysis Results Generated by STAARpipeline
Defines functions
Dynamic_Window_Results_Summary
Documented in
Dynamic_Window_Results_Summary
#' Summarize the results of dynamic window analysis generated by \code{STAARpipeline} package and
#' perform conditional analysis for (unconditionally) significant genetic regions by adjusting for a given list of known variants
#'
#' The \code{Dynamic_Window_Results_Summary} function takes in the results of dynamic window analysis generated by \code{STAARpipeline} package,
#' the object from fitting the null model, and the set of known variants to be adjusted for in conditional analysis
#' to summarize the dynamic window analysis results and analyze the conditional association between a quantitative/dichotomous phenotype and
#' the rare variants in the unconditional significant genetic regions.
#' @param agds_dir a vector containing file directory of annotated GDS (aGDS) files for all chromosomes (1-22).
#' @param jobs_num a data frame containing the number of jobs for association analysis.
#' The data frame must include a column with the name "scang_num"
#' @param input_path file directory of the input dynamic window analysis results.
#' @param output_path file directory of the output summary results.
#' @param dynamic_window_results_name file names of the input dynamic window analysis results.
#' @param obj_nullmodel an object from fitting the null model, which is either the output from \code{fit_nullmodel} function in the \code{STAARpipeline} package,
#' or the output from \code{fitNullModel} function in the \code{GENESIS} package and transformed using the \code{genesis2staar_nullmodel} function in the \code{STAARpipeline} package.
#' @param known_loci a data frame of variants to be adjusted for in conditional analysis and should
#' contain 4 columns in the following order: chromosome (CHR), position (POS), reference allele (REF),
#' and alternative allele (ALT) (default = NULL).
#' @param method_cond a character value indicating the method for conditional analysis.
#' \code{optimal} refers to regressing residuals from the null model on \code{known_loci}
#' as well as all covariates used in fitting the null model (fully adjusted) and taking the residuals;
#' \code{naive} refers to regressing residuals from the null model on \code{known_loci}
#' and taking the residuals (default = \code{optimal}).
#' @param geno_missing_imputation method of handling missing genotypes. Either "mean" or "minor" (default = "mean").
#' @param variant_type type of variant included in the conditional analysis. Choice includes "SNV", "Indel", or "variant" (default = "SNV").
#' @param QC_label channel name of the QC label in the GDS/aGDS file (default = "annotation/filter").
#' @param Annotation_dir channel name of the annotations in the aGDS file \cr (default = "annotation/info/FunctionalAnnotation").
#' @param Annotation_name_catalog a data frame containing the name and the corresponding channel name in the aGDS file.
#' @param Use_annotation_weights use annotations as weights or not (default = FALSE).
#' @param Annotation_name a vector of annotation names used in SCANG-STAAR (default = NULL).
#' @param alpha threshod to control the genome-wise (family-wise) error rate (default = 0.05).
#' @return The function returns the following analysis results:
#' @return \code{SCANG_S_res_uncond_cond.Rdata} and \code{SCANG_S_res_uncond_cond.csv}: A matrix that summarized the unconditional and
#' conditional results of the significant regions (GWER<alpha) detected by the SCANG-STAAR-S procedure (conditional results available if known_loci is not a NULL),
#' including chromosome ("chr"), start position ("start_pos"), end position ("end_pos"), number of variants ("SNV_nos"),
#' family-wise/genome-wide error rate (GWER), unconditional STAAR-S p-value ("STAAR_S"), conditional STAAR-S p-value ("STAAR_S_cond"),
#' conditional ACAT-V p-value ("ACAT_V_cond"), conditional Burden p-value ("Burden_cond"), conditional SKAT p-value ("SKAT_cond"),
#' and conditional STAAR-O p-value ("STAAR_O_cond").
#' @return \code{SCANG_B_res_uncond_cond.Rdata} and \code{SCANG_B_res_uncond_cond.csv}: A matrix that summarized the unconditional and
#' conditional results of the significant regions detected by the SCANG-STAAR-B procedure (conditional results available if known_loci is not a NULL).
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_O_res_uncond_cond.Rdata} and \code{SCANG_O_res_uncond_cond.csv}: A matrix that summarized the unconditional and
#' conditional results of the significant regions detected by the SCANG-STAAR-O procedure (conditional results available if known_loci is not a NULL).
#' Details see SCANG-STAAR-S.
#' @return \code{results_dynamic_window.Rdata}: A Rdata file that summarized the significant regions detected by SCANG-STAAR procedure.
#' @return \code{SCANG_S_top1.Rdata} and \code{SCANG_S_top1.csv}: A matrix that summarized the top 1 unconditional region detected by SCANG-STAAR-S,
#' including the STAAR-S p-value ("STAAR_S"), chromosome ("chr"), start position ("start_pos"), end position ("end_pos"),
#' family-wise/genome-wide error rate (GWER) and the number of variants ("SNV_nos").
#' @return \code{SCANG_B_top1.Rdata} and \code{SCANG_B_top1.csv}: A matrix that summarized the top 1 unconditional region detected by SCANG-STAAR-B.
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_O_top1.Rdata} and \code{SCANG_O_top1.csv}: A matrix that summarized the top 1 unconditional region detected by SCANG-STAAR-O.
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_S_res.Rdata} and \code{SCANG_S_res.csv}: A matrix that summarized the significant regions (GWER<alpha) detected by SCANG-STAAR-S,
#' including the negative log transformation of STAAR-S p-value ("-logp"), chromosome ("chr"), start position ("start_pos"), end position ("end_pos"),
#' family-wise/genome-wide error rate (GWER) and the number of variants ("SNV_num").
#' @return \code{SCANG_B_res.Rdata} and \code{SCANG_B_res.csv}: A matrix that summarized the significant regions detected by SCANG-STAAR-B.
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_O_res.Rdata} and \code{SCANG_O_res.csv}: A matrix that summarized the significant regions detected by SCANG-STAAR-O.
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_S_res_cond.Rdata} and \code{SCANG_S_res_cond.csv}: A matrix that summarized the conditional p-values of the significant regions (GWER<alpha) detected by SCANG-STAAR-S,
#' including chromosome ("chr"), start position ("Start Loc"), end position ("End Loc"),
#' the number of variants ("#SNV"), annotation-weighted ACAT-V, Burden and SKAT conditional p-values,
#' and STAAR conditional p-values of the regions with GWER smaller than the threshold alpha (available if known_loci is not a NULL).
#' @return \code{SCANG_B_res_cond.Rdata} and \code{SCANG_B_res_cond.csv}: A matrix that summarized the conditional p-values of the significant regions (GWER<alpha) detected by SCANG-STAAR-B (available if known_loci is not a NULL),
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_O_res_cond.Rdata} and \code{SCANG_O_res_cond.csv}: A matrix that summarized the conditional p-values of the significant regions (GWER<alpha) detected by SCANG-STAAR-O (available if known_loci is not a NULL),
#' Details see SCANG-STAAR-S.
#' @references Li, Z., Li, X., et al. (2022). A framework for detecting
#' noncoding rare-variant associations of large-scale whole-genome sequencing
#' studies. \emph{Nature Methods}, \emph{19}(12), 1599-1611.
#' (\href{https://doi.org/10.1038/s41592-022-01640-x}{pub})
#' @export
Dynamic_Window_Results_Summary <- function(agds_dir,jobs_num,input_path,output_path,dynamic_window_results_name,
obj_nullmodel,known_loci=NULL,
method_cond=c("optimal","naive"),
QC_label="annotation/filter",geno_missing_imputation=c("mean","minor"),variant_type=c("SNV","Indel","variant"),
Annotation_dir="annotation/info/FunctionalAnnotation",Annotation_name_catalog,
Use_annotation_weights=FALSE,Annotation_name=NULL,alpha=0.05){
## evaluate choices
method_cond <- match.arg(method_cond)
variant_type <- match.arg(variant_type)
geno_missing_imputation <- match.arg(geno_missing_imputation)
##### SCANG_O
SCANG_O_res <- c()
SCANG_O_top1 <- c()
SCANG_O_emthr <- c()
##### SCANG_S
SCANG_S_res <- c()
SCANG_S_top1 <- c()
SCANG_S_emthr <- c()
##### SCANG_B
SCANG_B_res <- c()
SCANG_B_top1 <- c()
SCANG_B_emthr <- c()
for(i in 1:sum(jobs_num$scang_num))
{
print(i)
results_scang <- get(load(paste0(input_path,dynamic_window_results_name,"_",i,".Rdata")))
# SCANG-O
SCANG_O_emthr <- rbind(SCANG_O_emthr,results_scang$SCANG_O_emthr)
SCANG_O_res <- rbind(SCANG_O_res,results_scang$SCANG_O_sig)
SCANG_O_top1 <- rbind(SCANG_O_top1,results_scang$SCANG_O_top1)
# SCANG-S
SCANG_S_emthr <- rbind(SCANG_S_emthr,results_scang$SCANG_S_emthr)
SCANG_S_res <- rbind(SCANG_S_res,results_scang$SCANG_S_sig)
SCANG_S_top1 <- rbind(SCANG_S_top1,results_scang$SCANG_S_top1)
# SCANG-B
SCANG_B_emthr <- rbind(SCANG_B_emthr,results_scang$SCANG_B_emthr)
SCANG_B_res <- rbind(SCANG_B_res,results_scang$SCANG_B_sig)
SCANG_B_top1 <- rbind(SCANG_B_top1,results_scang$SCANG_B_top1)
}
colnames(SCANG_O_res) <- colnames(SCANG_O_top1)
colnames(SCANG_S_res) <- colnames(SCANG_S_top1)
colnames(SCANG_B_res) <- colnames(SCANG_B_top1)
############ SCANG-O
# empirical threshold
emthr_genome <- apply(SCANG_O_emthr,2,max)
thr <- quantile(emthr_genome,1-alpha)
SCANG_O_res <- SCANG_O_res[SCANG_O_res[,1]>thr,,drop=FALSE]
if(length(SCANG_O_res)>=6)
{
for(i in 1:dim(SCANG_O_res)[1])
{
SCANG_O_res[i,5] <- mean(emthr_genome>SCANG_O_res[i,1])
}
SCANG_O_res <- SCANG_O_res[order(SCANG_O_res[,2]*1e10+SCANG_O_res[,3]),]
}
SCANG_O_top1 <- SCANG_O_top1[which.max(SCANG_O_top1[,1]),]
SCANG_O_top1[5] <- mean(emthr_genome>SCANG_O_top1[1])
thr_SCANG_O <- thr
emthr_genome_SCANG_O <- emthr_genome
############ SCANG-S
# empirical threshold
emthr_genome <- apply(SCANG_S_emthr,2,max)
thr <- quantile(emthr_genome,1-alpha)
SCANG_S_res <- SCANG_S_res[SCANG_S_res[,1]>thr,,drop=FALSE]
if(length(SCANG_S_res)>=6)
{
for(i in 1:dim(SCANG_S_res)[1])
{
SCANG_S_res[i,5] <- mean(emthr_genome>SCANG_S_res[i,1])
}
SCANG_S_res <- SCANG_S_res[order(SCANG_S_res[,2]*1e10+SCANG_S_res[,3]),]
}
SCANG_S_top1 <- SCANG_S_top1[which.max(SCANG_S_top1[,1]),]
SCANG_S_top1[5] <- mean(emthr_genome>SCANG_S_top1[1])
thr_SCANG_S <- thr
emthr_genome_SCANG_S <- emthr_genome
############ SCANG-B
# empirical threshold
emthr_genome <- apply(SCANG_B_emthr,2,max)
thr <- quantile(emthr_genome,1-alpha)
SCANG_B_res <- SCANG_B_res[SCANG_B_res[,1]>thr,,drop=FALSE]
if(length(SCANG_B_res)>=6)
{
for(i in 1:dim(SCANG_B_res)[1])
{
SCANG_B_res[i,5] <- mean(emthr_genome>SCANG_B_res[i,1])
}
SCANG_B_res <- SCANG_B_res[order(SCANG_B_res[,2]*1e10+SCANG_B_res[,3]),]
}
SCANG_B_top1 <- SCANG_B_top1[which.max(SCANG_B_top1[,1]),]
SCANG_B_top1[5] <- mean(emthr_genome>SCANG_B_top1[1])
thr_SCANG_B <- thr
emthr_genome_SCANG_B <- emthr_genome
SCANG_res <- list(SCANG_B_top1 = SCANG_B_top1, SCANG_B_emthr = emthr_genome_SCANG_B, SCANG_B_sig = SCANG_B_res,
SCANG_S_top1 = SCANG_S_top1, SCANG_S_emthr = emthr_genome_SCANG_S, SCANG_S_sig = SCANG_S_res,
SCANG_O_top1 = SCANG_O_top1, SCANG_O_emthr = emthr_genome_SCANG_O, SCANG_O_sig = SCANG_O_res,alpha=alpha)
# Output path
save(SCANG_res,file=paste0(output_path,"results_dynamic_window.Rdata"))
# SCANG_S
save(SCANG_S_res,file=paste0(output_path,"SCANG_S_res.Rdata"))
write.csv(SCANG_S_res,paste0(output_path,"SCANG_S_res.csv"))
# SCANG_O
save(SCANG_O_res,file=paste0(output_path,"SCANG_O_res.Rdata"))
write.csv(SCANG_O_res,paste0(output_path,"SCANG_O_res.csv"))
# SCANG_B
save(SCANG_B_res,file=paste0(output_path,"SCANG_B_res.Rdata"))
write.csv(SCANG_B_res,paste0(output_path,"SCANG_B_res.csv"))
### TOP 1 Results
# SCANG_S
SCANG_S_top1 <- SCANG_res$SCANG_S_top1
SCANG_S_top1[1] <- exp(-SCANG_S_top1[1])
SCANG_S_top1 <- matrix(SCANG_S_top1,nrow=1)
colnames(SCANG_S_top1) <- c("STAAR_S","chr","start_pos","end_pos","GWER","SNV_nos")
save(SCANG_S_top1,file=paste0(output_path,"SCANG_S_top1.Rdata"))
write.csv(SCANG_S_top1,paste0(output_path,"SCANG_S_top1.csv"))
# SCANG_B
SCANG_B_top1 <- SCANG_res$SCANG_B_top1
SCANG_B_top1[1] <- exp(-SCANG_B_top1[1])
SCANG_B_top1 <- matrix(SCANG_B_top1,nrow=1)
colnames(SCANG_B_top1) <- c("STAAR_B","chr","start_pos","end_pos","GWER","SNV_nos")
save(SCANG_B_top1,file=paste0(output_path,"SCANG_B_top1.Rdata"))
write.csv(SCANG_B_top1,paste0(output_path,"SCANG_B_top1.csv"))
# SCANG_O
SCANG_O_top1 <- SCANG_res$SCANG_O_top1
SCANG_O_top1[1] <- exp(-SCANG_O_top1[1])
SCANG_O_top1 <- matrix(SCANG_O_top1,nrow=1)
colnames(SCANG_O_top1) <- c("SCANG_STAAR_O","chr","start_pos","end_pos","GWER","SNV_nos")
save(SCANG_O_top1,file=paste0(output_path,"SCANG_O_top1.Rdata"))
write.csv(SCANG_O_top1,paste0(output_path,"SCANG_O_top1.csv"))
########################################################################
# Conditional Analysis
########################################################################
if(length(known_loci)!=0)
{
################### SCANG-S
SCANG_S_res <- SCANG_res$SCANG_S_sig
if(class(SCANG_S_res)[1]!="matrix")
{
SCANG_S_res <- matrix(SCANG_S_res,nrow=1)
}
SCANG_S_res_cond <- c()
if(length(SCANG_S_res)!=0)
{
for(kk in 1:dim(SCANG_S_res)[1])
{
chr <- as.numeric(SCANG_S_res[kk,2])
start_loc <- as.numeric(SCANG_S_res[kk,3])
end_loc <- as.numeric(SCANG_S_res[kk,4])
gds.path <- agds_dir[chr]
genofile <- seqOpen(gds.path)
res_cond <- Sliding_Window_cond(chr=chr,genofile=genofile,obj_nullmodel=obj_nullmodel,
start_loc=start_loc,end_loc=end_loc,known_loci=known_loci,method_cond=method_cond,
QC_label=QC_label,variant_type=variant_type,geno_missing_imputation=geno_missing_imputation,
Annotation_name_catalog=Annotation_name_catalog,Annotation_dir=Annotation_dir,
Use_annotation_weights=Use_annotation_weights,Annotation_name=Annotation_name)
SCANG_S_res_cond <- rbind(SCANG_S_res_cond,res_cond)
seqClose(genofile)
}
}
save(SCANG_S_res_cond,file=paste0(output_path,"SCANG_S_res_cond.Rdata"))
write.csv(SCANG_S_res_cond,paste0(output_path,"SCANG_S_res_cond.csv"))
### Uncond_Cond
SCANG_S_res_uncond_cond <- c()
if(length(SCANG_S_res)!=0)
{
## Uncond
SCANG_S_res[,1] <- exp(-SCANG_S_res[,1])
colnames(SCANG_S_res) <- c("STAAR_S","chr","start_pos","end_pos","GWER","SNV_nos")
SCANG_S_res <- SCANG_S_res[,c(2:4,6,5,1),drop=FALSE]
## SCANG_S_cond
STAAR_S <- SCANG_S_res_cond[,c("STAAR-S(1,25)","STAAR-S(1,1)"),drop=FALSE]
STAAR_S <- as.matrix(STAAR_S)
STAAR_S_cond <- apply(STAAR_S,1,function(z) CCT(as.numeric(z)))
STAAR_cond <- SCANG_S_res_cond[,c("ACAT-V(1,25)","Burden(1,1)","SKAT(1,25)","STAAR-O"),drop=FALSE]
## Uncond_Cond
SCANG_S_res_uncond_cond <- cbind(SCANG_S_res,STAAR_S_cond)
SCANG_S_res_uncond_cond <- cbind(SCANG_S_res_uncond_cond,STAAR_cond)
colnames(SCANG_S_res_uncond_cond)[(dim(SCANG_S_res_uncond_cond)[2]-3):dim(SCANG_S_res_uncond_cond)[2]] <- c("ACAT_V_cond","Burden_cond","SKAT_cond","STAAR_O_cond")
}
save(SCANG_S_res_uncond_cond,file=paste0(output_path,"SCANG_S_res_uncond_cond.Rdata"))
write.csv(SCANG_S_res_uncond_cond,paste0(output_path,"SCANG_S_res_uncond_cond.csv"))
######################## SCANG-O
SCANG_O_res <- SCANG_res$SCANG_O_sig
if(class(SCANG_O_res)[1]!="matrix")
{
SCANG_O_res <- matrix(SCANG_O_res,nrow=1)
}
SCANG_O_res_cond <- c()
if(length(SCANG_O_res)!=0)
{
for(kk in 1:dim(SCANG_O_res)[1])
{
chr <- as.numeric(SCANG_O_res[kk,2])
start_loc <- as.numeric(SCANG_O_res[kk,3])
end_loc <- as.numeric(SCANG_O_res[kk,4])
gds.path <- agds_dir[chr]
genofile <- seqOpen(gds.path)
res_cond <- Sliding_Window_cond(chr=chr,genofile=genofile,obj_nullmodel=obj_nullmodel,
start_loc=start_loc,end_loc=end_loc,known_loci=known_loci,method_cond=method_cond,
QC_label=QC_label,variant_type=variant_type,geno_missing_imputation=geno_missing_imputation,
Annotation_name_catalog=Annotation_name_catalog,Annotation_dir=Annotation_dir,
Use_annotation_weights=Use_annotation_weights,Annotation_name=Annotation_name)
SCANG_O_res_cond <- rbind(SCANG_O_res_cond,res_cond)
seqClose(genofile)
}
}
save(SCANG_O_res_cond,file=paste0(output_path,"SCANG_O_res_cond.Rdata"))
write.csv(SCANG_O_res_cond,paste0(output_path,"SCANG_O_res_cond.csv"))
### Uncond_Cond
SCANG_O_res_uncond_cond <- c()
if(length(SCANG_O_res)!=0)
{
## Uncond
SCANG_O_res[,1] <- exp(-SCANG_O_res[,1])
colnames(SCANG_O_res) <- c("SCANG_STAAR_O","chr","start_pos","end_pos","GWER","SNV_nos")
SCANG_O_res <- SCANG_O_res[,c(2:4,6,5,1),drop=FALSE]
## SCANG_O_cond
STAAR_O <- SCANG_O_res_cond[,c("STAAR-S(1,25)","STAAR-S(1,1)","STAAR-B(1,25)","STAAR-B(1,1)"),drop=FALSE]
STAAR_O <- as.matrix(STAAR_O)
SCANG_STAAR_O_cond <- apply(STAAR_O,1,function(z) CCT(as.numeric(z)))
STAAR_cond <- SCANG_O_res_cond[,c("ACAT-V(1,25)","Burden(1,1)","SKAT(1,25)","STAAR-O"),drop=FALSE]
SCANG_O_res_uncond_cond <- cbind(SCANG_O_res,SCANG_STAAR_O_cond)
SCANG_O_res_uncond_cond <- cbind(SCANG_O_res_uncond_cond,STAAR_cond)
colnames(SCANG_O_res_uncond_cond)[(dim(SCANG_O_res_uncond_cond)[2]-3):dim(SCANG_O_res_uncond_cond)[2]] <- c("ACAT_V_cond","Burden_cond","SKAT_cond","STAAR_O_cond")
}
save(SCANG_O_res_uncond_cond,file=paste0(output_path,"SCANG_O_res_uncond_cond.Rdata"))
write.csv(SCANG_O_res_uncond_cond,paste0(output_path,"SCANG_O_res_uncond_cond.csv"))
######################## SCANG-B
SCANG_B_res <- SCANG_res$SCANG_B_sig
if(class(SCANG_B_res)[1]!="matrix")
{
SCANG_B_res <- matrix(SCANG_B_res,nrow=1)
}
SCANG_B_res_cond <- c()
if(length(SCANG_B_res)!=0)
{
for(kk in 1:dim(SCANG_B_res)[1])
{
chr <- as.numeric(SCANG_B_res[kk,2])
start_loc <- as.numeric(SCANG_B_res[kk,3])
end_loc <- as.numeric(SCANG_B_res[kk,4])
gds.path <- agds_dir[chr]
genofile <- seqOpen(gds.path)
res_cond <- Sliding_Window_cond(chr=chr,genofile=genofile,obj_nullmodel=obj_nullmodel,
start_loc=start_loc,end_loc=end_loc,known_loci=known_loci,method_cond=method_cond,
QC_label=QC_label,variant_type=variant_type,geno_missing_imputation=geno_missing_imputation,
Annotation_name_catalog=Annotation_name_catalog,Annotation_dir=Annotation_dir,
Use_annotation_weights=Use_annotation_weights,Annotation_name=Annotation_name)
SCANG_B_res_cond <- rbind(SCANG_B_res_cond,res_cond)
seqClose(genofile)
}
}
save(SCANG_B_res_cond,file=paste0(output_path,"SCANG_B_res_cond.Rdata"))
write.csv(SCANG_B_res_cond,paste0(output_path,"SCANG_B_res_cond.csv"))
### Uncond_Cond
SCANG_B_res_uncond_cond <- c()
if(length(SCANG_B_res)!=0)
{
## Uncond
SCANG_B_res[,1] <- exp(-SCANG_B_res[,1])
colnames(SCANG_B_res) <- c("STAAR_B","chr","start_pos","end_pos","GWER","SNV_nos")
SCANG_B_res <- SCANG_B_res[,c(2:4,6,5,1),drop=FALSE]
## SCANG_B_cond
STAAR_B <- SCANG_B_res_cond[,c("STAAR-B(1,25)","STAAR-B(1,1)"),drop=FALSE]
STAAR_B <- as.matrix(STAAR_B)
STAAR_B_cond <- apply(STAAR_B,1,function(z) CCT(as.numeric(z)))
STAAR_cond <- SCANG_B_res_cond[,c("ACAT-V(1,25)","Burden(1,1)","SKAT(1,25)","STAAR-O"),drop=FALSE]
SCANG_B_res_uncond_cond <- cbind(SCANG_B_res,STAAR_B_cond)
SCANG_B_res_uncond_cond <- cbind(SCANG_B_res_uncond_cond,STAAR_cond)
colnames(SCANG_B_res_uncond_cond)[(dim(SCANG_B_res_uncond_cond)[2]-3):dim(SCANG_B_res_uncond_cond)[2]] <- c("ACAT_V_cond","Burden_cond","SKAT_cond","STAAR_O_cond")
}
save(SCANG_B_res_uncond_cond,file=paste0(output_path,"SCANG_B_res_uncond_cond.Rdata"))
write.csv(SCANG_B_res_uncond_cond,paste0(output_path,"SCANG_B_res_uncond_cond.csv"))
}
}
xihaoli/STAARpipelineSummary documentation built on Oct. 20, 2024, 9:35 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions Dynamic_Window_Results_Summary
Documented in Dynamic_Window_Results_Summary
#' Summarize the results of dynamic window analysis generated by \code{STAARpipeline} package and
#' perform conditional analysis for (unconditionally) significant genetic regions by adjusting for a given list of known variants
#'
#' The \code{Dynamic_Window_Results_Summary} function takes in the results of dynamic window analysis generated by \code{STAARpipeline} package,
#' the object from fitting the null model, and the set of known variants to be adjusted for in conditional analysis
#' to summarize the dynamic window analysis results and analyze the conditional association between a quantitative/dichotomous phenotype and
#' the rare variants in the unconditional significant genetic regions.
#' @param agds_dir a vector containing file directory of annotated GDS (aGDS) files for all chromosomes (1-22).
#' @param jobs_num a data frame containing the number of jobs for association analysis.
#' The data frame must include a column with the name "scang_num"
#' @param input_path file directory of the input dynamic window analysis results.
#' @param output_path file directory of the output summary results.
#' @param dynamic_window_results_name file names of the input dynamic window analysis results.
#' @param obj_nullmodel an object from fitting the null model, which is either the output from \code{fit_nullmodel} function in the \code{STAARpipeline} package,
#' or the output from \code{fitNullModel} function in the \code{GENESIS} package and transformed using the \code{genesis2staar_nullmodel} function in the \code{STAARpipeline} package.
#' @param known_loci a data frame of variants to be adjusted for in conditional analysis and should
#' contain 4 columns in the following order: chromosome (CHR), position (POS), reference allele (REF),
#' and alternative allele (ALT) (default = NULL).
#' @param method_cond a character value indicating the method for conditional analysis.
#' \code{optimal} refers to regressing residuals from the null model on \code{known_loci}
#' as well as all covariates used in fitting the null model (fully adjusted) and taking the residuals;
#' \code{naive} refers to regressing residuals from the null model on \code{known_loci}
#' and taking the residuals (default = \code{optimal}).
#' @param geno_missing_imputation method of handling missing genotypes. Either "mean" or "minor" (default = "mean").
#' @param variant_type type of variant included in the conditional analysis. Choice includes "SNV", "Indel", or "variant" (default = "SNV").
#' @param QC_label channel name of the QC label in the GDS/aGDS file (default = "annotation/filter").
#' @param Annotation_dir channel name of the annotations in the aGDS file \cr (default = "annotation/info/FunctionalAnnotation").
#' @param Annotation_name_catalog a data frame containing the name and the corresponding channel name in the aGDS file.
#' @param Use_annotation_weights use annotations as weights or not (default = FALSE).
#' @param Annotation_name a vector of annotation names used in SCANG-STAAR (default = NULL).
#' @param alpha threshod to control the genome-wise (family-wise) error rate (default = 0.05).
#' @return The function returns the following analysis results:
#' @return \code{SCANG_S_res_uncond_cond.Rdata} and \code{SCANG_S_res_uncond_cond.csv}: A matrix that summarized the unconditional and
#' conditional results of the significant regions (GWER<alpha) detected by the SCANG-STAAR-S procedure (conditional results available if known_loci is not a NULL),
#' including chromosome ("chr"), start position ("start_pos"), end position ("end_pos"), number of variants ("SNV_nos"),
#' family-wise/genome-wide error rate (GWER), unconditional STAAR-S p-value ("STAAR_S"), conditional STAAR-S p-value ("STAAR_S_cond"),
#' conditional ACAT-V p-value ("ACAT_V_cond"), conditional Burden p-value ("Burden_cond"), conditional SKAT p-value ("SKAT_cond"),
#' and conditional STAAR-O p-value ("STAAR_O_cond").
#' @return \code{SCANG_B_res_uncond_cond.Rdata} and \code{SCANG_B_res_uncond_cond.csv}: A matrix that summarized the unconditional and
#' conditional results of the significant regions detected by the SCANG-STAAR-B procedure (conditional results available if known_loci is not a NULL).
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_O_res_uncond_cond.Rdata} and \code{SCANG_O_res_uncond_cond.csv}: A matrix that summarized the unconditional and
#' conditional results of the significant regions detected by the SCANG-STAAR-O procedure (conditional results available if known_loci is not a NULL).
#' Details see SCANG-STAAR-S.
#' @return \code{results_dynamic_window.Rdata}: A Rdata file that summarized the significant regions detected by SCANG-STAAR procedure.
#' @return \code{SCANG_S_top1.Rdata} and \code{SCANG_S_top1.csv}: A matrix that summarized the top 1 unconditional region detected by SCANG-STAAR-S,
#' including the STAAR-S p-value ("STAAR_S"), chromosome ("chr"), start position ("start_pos"), end position ("end_pos"),
#' family-wise/genome-wide error rate (GWER) and the number of variants ("SNV_nos").
#' @return \code{SCANG_B_top1.Rdata} and \code{SCANG_B_top1.csv}: A matrix that summarized the top 1 unconditional region detected by SCANG-STAAR-B.
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_O_top1.Rdata} and \code{SCANG_O_top1.csv}: A matrix that summarized the top 1 unconditional region detected by SCANG-STAAR-O.
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_S_res.Rdata} and \code{SCANG_S_res.csv}: A matrix that summarized the significant regions (GWER<alpha) detected by SCANG-STAAR-S,
#' including the negative log transformation of STAAR-S p-value ("-logp"), chromosome ("chr"), start position ("start_pos"), end position ("end_pos"),
#' family-wise/genome-wide error rate (GWER) and the number of variants ("SNV_num").
#' @return \code{SCANG_B_res.Rdata} and \code{SCANG_B_res.csv}: A matrix that summarized the significant regions detected by SCANG-STAAR-B.
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_O_res.Rdata} and \code{SCANG_O_res.csv}: A matrix that summarized the significant regions detected by SCANG-STAAR-O.
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_S_res_cond.Rdata} and \code{SCANG_S_res_cond.csv}: A matrix that summarized the conditional p-values of the significant regions (GWER<alpha) detected by SCANG-STAAR-S,
#' including chromosome ("chr"), start position ("Start Loc"), end position ("End Loc"),
#' the number of variants ("#SNV"), annotation-weighted ACAT-V, Burden and SKAT conditional p-values,
#' and STAAR conditional p-values of the regions with GWER smaller than the threshold alpha (available if known_loci is not a NULL).
#' @return \code{SCANG_B_res_cond.Rdata} and \code{SCANG_B_res_cond.csv}: A matrix that summarized the conditional p-values of the significant regions (GWER<alpha) detected by SCANG-STAAR-B (available if known_loci is not a NULL),
#' Details see SCANG-STAAR-S.
#' @return \code{SCANG_O_res_cond.Rdata} and \code{SCANG_O_res_cond.csv}: A matrix that summarized the conditional p-values of the significant regions (GWER<alpha) detected by SCANG-STAAR-O (available if known_loci is not a NULL),
#' Details see SCANG-STAAR-S.
#' @references Li, Z., Li, X., et al. (2022). A framework for detecting
#' noncoding rare-variant associations of large-scale whole-genome sequencing
#' studies. \emph{Nature Methods}, \emph{19}(12), 1599-1611.
#' (\href{https://doi.org/10.1038/s41592-022-01640-x}{pub})
#' @export
Dynamic_Window_Results_Summary <- function(agds_dir,jobs_num,input_path,output_path,dynamic_window_results_name,
obj_nullmodel,known_loci=NULL,
method_cond=c("optimal","naive"),
QC_label="annotation/filter",geno_missing_imputation=c("mean","minor"),variant_type=c("SNV","Indel","variant"),
Annotation_dir="annotation/info/FunctionalAnnotation",Annotation_name_catalog,
Use_annotation_weights=FALSE,Annotation_name=NULL,alpha=0.05){
## evaluate choices
method_cond <- match.arg(method_cond)
variant_type <- match.arg(variant_type)
geno_missing_imputation <- match.arg(geno_missing_imputation)
##### SCANG_O
SCANG_O_res <- c()
SCANG_O_top1 <- c()
SCANG_O_emthr <- c()
##### SCANG_S
SCANG_S_res <- c()
SCANG_S_top1 <- c()
SCANG_S_emthr <- c()
##### SCANG_B
SCANG_B_res <- c()
SCANG_B_top1 <- c()
SCANG_B_emthr <- c()
for(i in 1:sum(jobs_num$scang_num))
{
print(i)
results_scang <- get(load(paste0(input_path,dynamic_window_results_name,"_",i,".Rdata")))
# SCANG-O
SCANG_O_emthr <- rbind(SCANG_O_emthr,results_scang$SCANG_O_emthr)
SCANG_O_res <- rbind(SCANG_O_res,results_scang$SCANG_O_sig)
SCANG_O_top1 <- rbind(SCANG_O_top1,results_scang$SCANG_O_top1)
# SCANG-S
SCANG_S_emthr <- rbind(SCANG_S_emthr,results_scang$SCANG_S_emthr)
SCANG_S_res <- rbind(SCANG_S_res,results_scang$SCANG_S_sig)
SCANG_S_top1 <- rbind(SCANG_S_top1,results_scang$SCANG_S_top1)
# SCANG-B
SCANG_B_emthr <- rbind(SCANG_B_emthr,results_scang$SCANG_B_emthr)
SCANG_B_res <- rbind(SCANG_B_res,results_scang$SCANG_B_sig)
SCANG_B_top1 <- rbind(SCANG_B_top1,results_scang$SCANG_B_top1)
}
colnames(SCANG_O_res) <- colnames(SCANG_O_top1)
colnames(SCANG_S_res) <- colnames(SCANG_S_top1)
colnames(SCANG_B_res) <- colnames(SCANG_B_top1)
############ SCANG-O
# empirical threshold
emthr_genome <- apply(SCANG_O_emthr,2,max)
thr <- quantile(emthr_genome,1-alpha)
SCANG_O_res <- SCANG_O_res[SCANG_O_res[,1]>thr,,drop=FALSE]
if(length(SCANG_O_res)>=6)
{
for(i in 1:dim(SCANG_O_res)[1])
{
SCANG_O_res[i,5] <- mean(emthr_genome>SCANG_O_res[i,1])
}
SCANG_O_res <- SCANG_O_res[order(SCANG_O_res[,2]*1e10+SCANG_O_res[,3]),]
}
SCANG_O_top1 <- SCANG_O_top1[which.max(SCANG_O_top1[,1]),]
SCANG_O_top1[5] <- mean(emthr_genome>SCANG_O_top1[1])
thr_SCANG_O <- thr
emthr_genome_SCANG_O <- emthr_genome
############ SCANG-S
# empirical threshold
emthr_genome <- apply(SCANG_S_emthr,2,max)
thr <- quantile(emthr_genome,1-alpha)
SCANG_S_res <- SCANG_S_res[SCANG_S_res[,1]>thr,,drop=FALSE]
if(length(SCANG_S_res)>=6)
{
for(i in 1:dim(SCANG_S_res)[1])
{
SCANG_S_res[i,5] <- mean(emthr_genome>SCANG_S_res[i,1])
}
SCANG_S_res <- SCANG_S_res[order(SCANG_S_res[,2]*1e10+SCANG_S_res[,3]),]
}
SCANG_S_top1 <- SCANG_S_top1[which.max(SCANG_S_top1[,1]),]
SCANG_S_top1[5] <- mean(emthr_genome>SCANG_S_top1[1])
thr_SCANG_S <- thr
emthr_genome_SCANG_S <- emthr_genome
############ SCANG-B
# empirical threshold
emthr_genome <- apply(SCANG_B_emthr,2,max)
thr <- quantile(emthr_genome,1-alpha)
SCANG_B_res <- SCANG_B_res[SCANG_B_res[,1]>thr,,drop=FALSE]
if(length(SCANG_B_res)>=6)
{
for(i in 1:dim(SCANG_B_res)[1])
{
SCANG_B_res[i,5] <- mean(emthr_genome>SCANG_B_res[i,1])
}
SCANG_B_res <- SCANG_B_res[order(SCANG_B_res[,2]*1e10+SCANG_B_res[,3]),]
}
SCANG_B_top1 <- SCANG_B_top1[which.max(SCANG_B_top1[,1]),]
SCANG_B_top1[5] <- mean(emthr_genome>SCANG_B_top1[1])
thr_SCANG_B <- thr
emthr_genome_SCANG_B <- emthr_genome
SCANG_res <- list(SCANG_B_top1 = SCANG_B_top1, SCANG_B_emthr = emthr_genome_SCANG_B, SCANG_B_sig = SCANG_B_res,
SCANG_S_top1 = SCANG_S_top1, SCANG_S_emthr = emthr_genome_SCANG_S, SCANG_S_sig = SCANG_S_res,
SCANG_O_top1 = SCANG_O_top1, SCANG_O_emthr = emthr_genome_SCANG_O, SCANG_O_sig = SCANG_O_res,alpha=alpha)
# Output path
save(SCANG_res,file=paste0(output_path,"results_dynamic_window.Rdata"))
# SCANG_S
save(SCANG_S_res,file=paste0(output_path,"SCANG_S_res.Rdata"))
write.csv(SCANG_S_res,paste0(output_path,"SCANG_S_res.csv"))
# SCANG_O
save(SCANG_O_res,file=paste0(output_path,"SCANG_O_res.Rdata"))
write.csv(SCANG_O_res,paste0(output_path,"SCANG_O_res.csv"))
# SCANG_B
save(SCANG_B_res,file=paste0(output_path,"SCANG_B_res.Rdata"))
write.csv(SCANG_B_res,paste0(output_path,"SCANG_B_res.csv"))
### TOP 1 Results
# SCANG_S
SCANG_S_top1 <- SCANG_res$SCANG_S_top1
SCANG_S_top1[1] <- exp(-SCANG_S_top1[1])
SCANG_S_top1 <- matrix(SCANG_S_top1,nrow=1)
colnames(SCANG_S_top1) <- c("STAAR_S","chr","start_pos","end_pos","GWER","SNV_nos")
save(SCANG_S_top1,file=paste0(output_path,"SCANG_S_top1.Rdata"))
write.csv(SCANG_S_top1,paste0(output_path,"SCANG_S_top1.csv"))
# SCANG_B
SCANG_B_top1 <- SCANG_res$SCANG_B_top1
SCANG_B_top1[1] <- exp(-SCANG_B_top1[1])
SCANG_B_top1 <- matrix(SCANG_B_top1,nrow=1)
colnames(SCANG_B_top1) <- c("STAAR_B","chr","start_pos","end_pos","GWER","SNV_nos")
save(SCANG_B_top1,file=paste0(output_path,"SCANG_B_top1.Rdata"))
write.csv(SCANG_B_top1,paste0(output_path,"SCANG_B_top1.csv"))
# SCANG_O
SCANG_O_top1 <- SCANG_res$SCANG_O_top1
SCANG_O_top1[1] <- exp(-SCANG_O_top1[1])
SCANG_O_top1 <- matrix(SCANG_O_top1,nrow=1)
colnames(SCANG_O_top1) <- c("SCANG_STAAR_O","chr","start_pos","end_pos","GWER","SNV_nos")
save(SCANG_O_top1,file=paste0(output_path,"SCANG_O_top1.Rdata"))
write.csv(SCANG_O_top1,paste0(output_path,"SCANG_O_top1.csv"))
########################################################################
# Conditional Analysis
########################################################################
if(length(known_loci)!=0)
{
################### SCANG-S
SCANG_S_res <- SCANG_res$SCANG_S_sig
if(class(SCANG_S_res)[1]!="matrix")
{
SCANG_S_res <- matrix(SCANG_S_res,nrow=1)
}
SCANG_S_res_cond <- c()
if(length(SCANG_S_res)!=0)
{
for(kk in 1:dim(SCANG_S_res)[1])
{
chr <- as.numeric(SCANG_S_res[kk,2])
start_loc <- as.numeric(SCANG_S_res[kk,3])
end_loc <- as.numeric(SCANG_S_res[kk,4])
gds.path <- agds_dir[chr]
genofile <- seqOpen(gds.path)
res_cond <- Sliding_Window_cond(chr=chr,genofile=genofile,obj_nullmodel=obj_nullmodel,
start_loc=start_loc,end_loc=end_loc,known_loci=known_loci,method_cond=method_cond,
QC_label=QC_label,variant_type=variant_type,geno_missing_imputation=geno_missing_imputation,
Annotation_name_catalog=Annotation_name_catalog,Annotation_dir=Annotation_dir,
Use_annotation_weights=Use_annotation_weights,Annotation_name=Annotation_name)
SCANG_S_res_cond <- rbind(SCANG_S_res_cond,res_cond)
seqClose(genofile)
}
}
save(SCANG_S_res_cond,file=paste0(output_path,"SCANG_S_res_cond.Rdata"))
write.csv(SCANG_S_res_cond,paste0(output_path,"SCANG_S_res_cond.csv"))
### Uncond_Cond
SCANG_S_res_uncond_cond <- c()
if(length(SCANG_S_res)!=0)
{
## Uncond
SCANG_S_res[,1] <- exp(-SCANG_S_res[,1])
colnames(SCANG_S_res) <- c("STAAR_S","chr","start_pos","end_pos","GWER","SNV_nos")
SCANG_S_res <- SCANG_S_res[,c(2:4,6,5,1),drop=FALSE]
## SCANG_S_cond
STAAR_S <- SCANG_S_res_cond[,c("STAAR-S(1,25)","STAAR-S(1,1)"),drop=FALSE]
STAAR_S <- as.matrix(STAAR_S)
STAAR_S_cond <- apply(STAAR_S,1,function(z) CCT(as.numeric(z)))
STAAR_cond <- SCANG_S_res_cond[,c("ACAT-V(1,25)","Burden(1,1)","SKAT(1,25)","STAAR-O"),drop=FALSE]
## Uncond_Cond
SCANG_S_res_uncond_cond <- cbind(SCANG_S_res,STAAR_S_cond)
SCANG_S_res_uncond_cond <- cbind(SCANG_S_res_uncond_cond,STAAR_cond)
colnames(SCANG_S_res_uncond_cond)[(dim(SCANG_S_res_uncond_cond)[2]-3):dim(SCANG_S_res_uncond_cond)[2]] <- c("ACAT_V_cond","Burden_cond","SKAT_cond","STAAR_O_cond")
}
save(SCANG_S_res_uncond_cond,file=paste0(output_path,"SCANG_S_res_uncond_cond.Rdata"))
write.csv(SCANG_S_res_uncond_cond,paste0(output_path,"SCANG_S_res_uncond_cond.csv"))
######################## SCANG-O
SCANG_O_res <- SCANG_res$SCANG_O_sig
if(class(SCANG_O_res)[1]!="matrix")
{
SCANG_O_res <- matrix(SCANG_O_res,nrow=1)
}
SCANG_O_res_cond <- c()
if(length(SCANG_O_res)!=0)
{
for(kk in 1:dim(SCANG_O_res)[1])
{
chr <- as.numeric(SCANG_O_res[kk,2])
start_loc <- as.numeric(SCANG_O_res[kk,3])
end_loc <- as.numeric(SCANG_O_res[kk,4])
gds.path <- agds_dir[chr]
genofile <- seqOpen(gds.path)
res_cond <- Sliding_Window_cond(chr=chr,genofile=genofile,obj_nullmodel=obj_nullmodel,
start_loc=start_loc,end_loc=end_loc,known_loci=known_loci,method_cond=method_cond,
QC_label=QC_label,variant_type=variant_type,geno_missing_imputation=geno_missing_imputation,
Annotation_name_catalog=Annotation_name_catalog,Annotation_dir=Annotation_dir,
Use_annotation_weights=Use_annotation_weights,Annotation_name=Annotation_name)
SCANG_O_res_cond <- rbind(SCANG_O_res_cond,res_cond)
seqClose(genofile)
}
}
save(SCANG_O_res_cond,file=paste0(output_path,"SCANG_O_res_cond.Rdata"))
write.csv(SCANG_O_res_cond,paste0(output_path,"SCANG_O_res_cond.csv"))
### Uncond_Cond
SCANG_O_res_uncond_cond <- c()
if(length(SCANG_O_res)!=0)
{
## Uncond
SCANG_O_res[,1] <- exp(-SCANG_O_res[,1])
colnames(SCANG_O_res) <- c("SCANG_STAAR_O","chr","start_pos","end_pos","GWER","SNV_nos")
SCANG_O_res <- SCANG_O_res[,c(2:4,6,5,1),drop=FALSE]
## SCANG_O_cond
STAAR_O <- SCANG_O_res_cond[,c("STAAR-S(1,25)","STAAR-S(1,1)","STAAR-B(1,25)","STAAR-B(1,1)"),drop=FALSE]
STAAR_O <- as.matrix(STAAR_O)
SCANG_STAAR_O_cond <- apply(STAAR_O,1,function(z) CCT(as.numeric(z)))
STAAR_cond <- SCANG_O_res_cond[,c("ACAT-V(1,25)","Burden(1,1)","SKAT(1,25)","STAAR-O"),drop=FALSE]
SCANG_O_res_uncond_cond <- cbind(SCANG_O_res,SCANG_STAAR_O_cond)
SCANG_O_res_uncond_cond <- cbind(SCANG_O_res_uncond_cond,STAAR_cond)
colnames(SCANG_O_res_uncond_cond)[(dim(SCANG_O_res_uncond_cond)[2]-3):dim(SCANG_O_res_uncond_cond)[2]] <- c("ACAT_V_cond","Burden_cond","SKAT_cond","STAAR_O_cond")
}
save(SCANG_O_res_uncond_cond,file=paste0(output_path,"SCANG_O_res_uncond_cond.Rdata"))
write.csv(SCANG_O_res_uncond_cond,paste0(output_path,"SCANG_O_res_uncond_cond.csv"))
######################## SCANG-B
SCANG_B_res <- SCANG_res$SCANG_B_sig
if(class(SCANG_B_res)[1]!="matrix")
{
SCANG_B_res <- matrix(SCANG_B_res,nrow=1)
}
SCANG_B_res_cond <- c()
if(length(SCANG_B_res)!=0)
{
for(kk in 1:dim(SCANG_B_res)[1])
{
chr <- as.numeric(SCANG_B_res[kk,2])
start_loc <- as.numeric(SCANG_B_res[kk,3])
end_loc <- as.numeric(SCANG_B_res[kk,4])
gds.path <- agds_dir[chr]
genofile <- seqOpen(gds.path)
res_cond <- Sliding_Window_cond(chr=chr,genofile=genofile,obj_nullmodel=obj_nullmodel,
start_loc=start_loc,end_loc=end_loc,known_loci=known_loci,method_cond=method_cond,
QC_label=QC_label,variant_type=variant_type,geno_missing_imputation=geno_missing_imputation,
Annotation_name_catalog=Annotation_name_catalog,Annotation_dir=Annotation_dir,
Use_annotation_weights=Use_annotation_weights,Annotation_name=Annotation_name)
SCANG_B_res_cond <- rbind(SCANG_B_res_cond,res_cond)
seqClose(genofile)
}
}
save(SCANG_B_res_cond,file=paste0(output_path,"SCANG_B_res_cond.Rdata"))
write.csv(SCANG_B_res_cond,paste0(output_path,"SCANG_B_res_cond.csv"))
### Uncond_Cond
SCANG_B_res_uncond_cond <- c()
if(length(SCANG_B_res)!=0)
{
## Uncond
SCANG_B_res[,1] <- exp(-SCANG_B_res[,1])
colnames(SCANG_B_res) <- c("STAAR_B","chr","start_pos","end_pos","GWER","SNV_nos")
SCANG_B_res <- SCANG_B_res[,c(2:4,6,5,1),drop=FALSE]
## SCANG_B_cond
STAAR_B <- SCANG_B_res_cond[,c("STAAR-B(1,25)","STAAR-B(1,1)"),drop=FALSE]
STAAR_B <- as.matrix(STAAR_B)
STAAR_B_cond <- apply(STAAR_B,1,function(z) CCT(as.numeric(z)))
STAAR_cond <- SCANG_B_res_cond[,c("ACAT-V(1,25)","Burden(1,1)","SKAT(1,25)","STAAR-O"),drop=FALSE]
SCANG_B_res_uncond_cond <- cbind(SCANG_B_res,STAAR_B_cond)
SCANG_B_res_uncond_cond <- cbind(SCANG_B_res_uncond_cond,STAAR_cond)
colnames(SCANG_B_res_uncond_cond)[(dim(SCANG_B_res_uncond_cond)[2]-3):dim(SCANG_B_res_uncond_cond)[2]] <- c("ACAT_V_cond","Burden_cond","SKAT_cond","STAAR_O_cond")
}
save(SCANG_B_res_uncond_cond,file=paste0(output_path,"SCANG_B_res_uncond_cond.Rdata"))
write.csv(SCANG_B_res_uncond_cond,paste0(output_path,"SCANG_B_res_uncond_cond.csv"))
}
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.