Last updated: 2020-03-18
Checks: 7 0
Knit directory: m6AQTL_reproducibleDocument/
This reproducible R Markdown analysis was created with workflowr (version 1.6.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20200317) was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish or wflow_git_commit). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Ignored files:
Ignored: .Rhistory
Ignored: .Rproj.user/
Note that any generated files, e.g. HTML, png, CSS, etc., are not included in this status report because it is ok for generated content to have uncommitted changes.
These are the previous versions of the R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.
| File | Version | Author | Date | Message |
|---|---|---|---|---|
| Rmd | ed546ad | kevinlkx | 2020-03-18 | wflow_publish(“analysis/fine_mapping_m6AQTLs.Rmd”) |
| html | 2e78d20 | kevinlkx | 2020-03-18 | Build site. |
| Rmd | 5d10f06 | kevinlkx | 2020-03-18 | wflow_publish(“analysis/fine_mapping_m6AQTLs.Rmd”) |
m6AQTL_workflowr/analysis/prepare_susie_input_allPeaks.RmdsuppressMessages(library(dplyr))
suppressMessages(library(doParallel))
suppressMessages(library(susieR))
suppressMessages(library(vcfR))
get_genotype_dosage <- function(vcf_file, SNP_pos, iPeak){
library(vcfR)
## Helper function to convert genotype into dosage.
.genoDosage <- function(x){
return( stringr::str_count(x,"1") )
}
SNP_chr <- sapply(strsplit(as.character(SNP_pos), ":"), "[[", 1)[1]
SNP_pos_range <- range(as.integer(sapply(strsplit(as.character(SNP_pos), ":"), "[[", 2)))
SNP_range <- data.frame(chr = SNP_chr, start = SNP_pos_range[1], end = SNP_pos_range[2])
# tmpVcf <- tempfile(fileext = ".vcf.gz")
# SNP_chr_pos_range <- paste0(SNP_range$chr, ":", SNP_range$start, "-", SNP_range$end)
# system(paste("bcftools view -r", SNP_chr_pos_range, vcf_file, "-O z -o", tmpVcf))
# geno.vcf <- try( read.vcfR(file = tmpVcf, verbose = F), silent = T)
# unlink(tmpVcf) # remove the temp file to free space
## Use unix command line to accesss the genotype from vcf file. This is for fast data access.
tmpVcf <- tempfile(fileext = ".vcf.gz")
system(paste0("zcat ",vcf_file," | awk 'NR==1 {print $0} (/^#CHROM/){print $0} (!/^#/ && $2 >= ", SNP_range$start," && $2 <= ",SNP_range$end," ) {print $0}'| gzip > ",tmpVcf))
geno.vcf <- try( read.vcfR( file =tmpVcf, verbose = F ) , silent = T)
### This is to handle a wired error in the read.vcfR function.
if(class(geno.vcf) == "try-error"){
system(paste0("zcat ",vcf_file," | awk '/^#/ {print $0} (!/^#/ && $2 >= ",SNP_range$start," && $2 <= ",SNP_range$end," ) {print $0}'| gzip >",tmpVcf))
geno.vcf <- read.vcfR( file =tmpVcf, verbose = F )
}
unlink(tmpVcf) # remove the temp file to free space.
## check if genotype available for this peak
if( nrow(geno.vcf@fix) == 0 ){
cat(paste("Warning:", iPeak, "no SNPs in vcf! \n"))
return(NULL)
}
if( nrow(geno.vcf@fix) > length(SNP_pos)){
cat("Warning:", iPeak, "#SNPs in vcf > #SNP from summary stats! \n")
}else if( nrow(geno.vcf@fix) < length(SNP_pos)){
cat("Warning:", iPeak, "#SNPs in vcf < #SNP from summary stats!! \n")
}
## Get genotype as dosage format and filter for MAF
if( unique(geno.vcf@gt[,"FORMAT"]) == "DS" ){
## Directly extract dosage
geno <- if(nrow(geno.vcf@fix) == 1){
t(apply(extract.gt(geno.vcf, element = 'DS' ),2,as.numeric ) )
}else{
apply(extract.gt(geno.vcf, element = 'DS' ),2,as.numeric )
}
# rownames(geno) <- geno.vcf@fix[,"ID"]
rownames(geno) <- paste0(geno.vcf@fix[,"CHROM"], ":", geno.vcf@fix[,"POS"])
## filter out any genotype that has MAF<0.05
MAF <- apply(geno,1,function(x) !any(table(round(x) )>0.95*ncol(geno)) )
geno <- geno[MAF,]
geno.vcf <- geno.vcf[MAF,]
}else if(unique(geno.vcf@gt[,"FORMAT"]) == "GT"){
geno.vcf <- geno.vcf[is.biallelic(geno.vcf),]
## get genotype as Dosage
tmp_geno <- extract.gt(geno.vcf, element = 'GT' )
geno <- t( apply( tmp_geno ,1, .genoDosage ) )
colnames(geno) <- colnames(tmp_geno)
rownames(geno) <- paste0(geno.vcf@fix[,"CHROM"], ":", geno.vcf@fix[,"POS"])
## filter out any genotype that has MAF<0.05
MAF <- apply(geno,1,function(x) !any(table(x)>0.95*ncol(geno)) )
geno <- geno[MAF,]
geno.vcf <- geno.vcf[MAF,]
}
if( nrow(geno) == 0 ){ return(NULL) } ## skip this iteration if no genotype left after filtering
return(list(geno = geno, geno.vcf = geno.vcf))
}
### Load m6A-QTL full individual data
m6A_version <- "jointPeak_threshold5_MeRIP_HISAT2Map"
m6A_phenotype_name <- "m6APeak_logOR_GC.IP.adjusted_qqnorm"
num_PCs <- 15
dir_fineMapping <- paste0("/project2/xinhe/m6A/data_shared/m6A_QTLs/", m6A_version, "/fineMapping/", m6A_phenotype_name, "/m6AQTL_allPeaks")
dir.create(dir_fineMapping, showWarnings = F, recursive = T)
dir_genotype <- "/project2/xinhe/m6A/m6A_seq/m6A_QTL/genotypes/hg19/YRI/sixtySample_hg19_201903"
dir_m6A_phenotype <- paste0("/project2/xinhe/m6A/m6A_seq/m6A_QTL/peakcalling/", m6A_version, "/", m6A_phenotype_name)
cat("m6A_version:", m6A_version, "\n")
cat("m6A_phenotype_name:", m6A_phenotype_name, "\n")
cat("Directory of fine-mapping result:", dir_fineMapping, "\n")
dir_m6AQTL <- paste0("/project2/xinhe/m6A/data_shared/m6A_QTLs/", m6A_version, "/fastQTL_cis_100kb/", m6A_phenotype_name)
m6AQTL <- readRDS(paste0(dir_m6AQTL, "/m6AQTL.", m6A_phenotype_name, ".", num_PCs, "PCs.fastQTL.nominals.rds"))
m6A_phenotype <- read.table(paste0(dir_m6A_phenotype, "/", m6A_phenotype_name, ".fastQTL.txt.gz"), header = T, comment.char = "", stringsAsFactors = F)
colnames(m6A_phenotype)[1] <- "chr"
m6A_phenotype_matrix <- t(as.matrix(m6A_phenotype[, -c(1:4)]))
colnames(m6A_phenotype_matrix) <- m6A_phenotype$PEAK
PCs <- read.table(paste0(dir_m6A_phenotype, "/PCs/", num_PCs, "PCs.txt.gz"), header = T, comment.char = "", stringsAsFactors = F)
PCs_matrix <- t(as.matrix(PCs[,-1]))
colnames(PCs_matrix) <- PCs[,1]
### Preparing Susie input data
m6AQTL_PEAKs <- unique(m6AQTL$PEAK)
cat("Prepare Susie input for", length(m6AQTL_PEAKs), "peaks. \n")
registerDoParallel(cores=min(15, detectCores()-1))
foreach( iPeak=m6AQTL_PEAKs, .combine = function(a,b){ NULL } ) %dopar% {
cat("Peak: ",iPeak, "\n")
peak_bed <- m6A_phenotype[m6A_phenotype$PEAK == iPeak, 1:4]
peakStat <- dplyr::filter(m6AQTL, PEAK == iPeak)
# phenotype vector: n x 1
y <- m6A_phenotype_matrix[, iPeak]
# genotype matrix: n x p, p: number of SNPs
vcf_file <- paste0(dir_genotype, "/Imputed.sixty.hg19.", peak_bed$chr, ".vcf.gz")
geno.l <- get_genotype_dosage(vcf_file, peakStat$SNP, iPeak)
X <- t(as.matrix(geno.l$geno))
## sort individuals
X <- X[names(y), ]
PCs_matrix <- PCs_matrix[names(y), ]
## select common SNPs and sort SNPs
SNP_order <- intersect(colnames(X), peakStat$SNP)
X <- X[, match(SNP_order, colnames(X))]
peakStat <- peakStat[match(SNP_order, peakStat$SNP), ]
est <- peakStat[,c("SNP", "beta", "se", "pvalue")]
dir.create(paste0(dir_fineMapping, "/INPUT/", iPeak), showWarnings = F, recursive = T)
write.table(est, gzfile(paste0(dir_fineMapping, "/INPUT/", iPeak, "/m6AQTLs.est.txt.gz")), sep = "\t", col.names = F, row.names = F, quote = F)
saveRDS(list(X = X, y = y, PCs = PCs_matrix), file = paste0(dir_fineMapping, "/INPUT/", iPeak, "/m6AQTL_full_data.RDS"))
}
cat("Finished preparing susie input. \n")
## Check for any peaks without data
toFinish <- foreach( iPeak = m6AQTL_PEAKs, .combine = c )%dopar%{
if( file.exists( paste0(dir_fineMapping, "/INPUT/", iPeak, "/m6AQTL_full_data.RDS") ) ){
return(NULL)
}else{
return(iPeak)
}
}
cat("Peaks without data: \n")
print(toFinish)m6AQTL_workflowr/analysis/susie_allPeaks_L3_estPriorVar.Rmd# Fine-mapping m6AQTLs using SuSiE on m6A all peaks (L = 3, estimate prior variance)
suppressMessages(library(dplyr))
suppressMessages(library(doParallel))
suppressMessages(library(susieR))
m6A_version <- "jointPeak_threshold5_MeRIP_HISAT2Map"
m6A_phenotype_name <- "m6APeak_logOR_GC.IP.adjusted_qqnorm"
num_PCs <- 15
L_Susie <- 3
susie_setting <- "Susie_L3_estPriorVar"
dir_fineMapping <- paste0("/project2/xinhe/m6A/data_shared/m6A_QTLs/", m6A_version, "/fineMapping/", m6A_phenotype_name, "/m6AQTL_allPeaks")
dir_susie_output <- paste0(dir_fineMapping, "/", susie_setting)
dir.create(dir_susie_output, showWarnings = F, recursive = T)
num_cores <- min(15, detectCores() - 1)
cat("m6A_version:", m6A_version, "\n")
cat("m6A_phenotype_name:", m6A_phenotype_name, "\n")
cat("Susie setting:", susie_setting, "\n")
cat(num_PCs, "PCs included \n")
cat("Directory of SuSiE input data:", paste0(dir_fineMapping, "/INPUT/"), "\n")
cat("Directory of SuSiE output results:", dir_susie_output, "\n")
## Load m6AQTL summary statistics
dir_m6AQTL <- paste0("/project2/xinhe/m6A/data_shared/m6A_QTLs/", m6A_version, "/fastQTL_cis_100kb/", m6A_phenotype_name)
m6AQTL <- readRDS(paste0(dir_m6AQTL, "/m6AQTL.", m6A_phenotype_name, ".", num_PCs, "PCs.fastQTL.nominals.rds"))
m6AQTL_PEAKs <- unique(m6AQTL$PEAK)
peak_map <- data.frame(peakName = paste0("Peak",1:length(m6AQTL_PEAKs)), peakID = m6AQTL_PEAKs, row.names = m6AQTL_PEAKs, stringsAsFactors = FALSE)
rm(m6AQTL)
### check for peaks without individual data
cat(length(m6AQTL_PEAKs), "peaks \n")
toFinish <- foreach( iPeak = m6AQTL_PEAKs, .combine = c )%dopar%{
if( file.exists( paste0(dir_fineMapping, "/INPUT/", iPeak, "/m6AQTL_full_data.RDS") ) ){
return(NULL)
}else{
return(iPeak)
}
}
cat("Peaks without individual data: \n")
print(toFinish)
## Run Susie with and without prior, and save PVE and zscores
cat("Finemapping: ",length(peak_map$peakID), "peaks \n")
dir.create(paste0(dir_susie_output, "/Peaks/"), showWarnings = F, recursive = T)
registerDoParallel(num_cores)
res <- foreach( iPeak=peak_map$peakID )%dopar%{
# cat(iPeak, "\r")
data <- readRDS(paste0(dir_fineMapping, "/INPUT/", iPeak, "/m6AQTL_full_data.RDS"))
# est <- read.table(gzfile(paste0(dir_fineMapping, "/INPUT/", iPeak, "/m6AQTLs.est.txt.gz")), sep = "\t")
# colnames(est) <- c("SNP", "beta", "se", "pvalue")
X <- data$X
y <- data$y
PCs <- data$PCs
# covariate matrix: n x (k+1), including intercept
Z <- cbind(intercept = 1, PCs)
######## removing covariates from X ########
A <- crossprod(Z)
# chol decomposition for (Z'Z)^(-1)
R <- chol(solve(A))
W <- R %*% t(Z) %*% X
# Remove Covariates from X
Xhat <- as.matrix(X - Z %*% crossprod(R,W))
# LD
R_LD <- cor(Xhat)
######## removing covariates from Y ########
A <- crossprod(Z)
SZy <- solve(A, drop(y %*% Z))
yhat <- y - drop(Z %*% SZy)
## compute PVE
reg <- susieR:::univariate_regression(Xhat, yhat)
z_score <- reg$betahat/reg$sebetahat
names(z_score) <- colnames(Xhat)
top_idx <- which.max(abs(z_score))
PVE <- var(Xhat[,top_idx] * reg$betahat[top_idx]) / var(yhat)
cat(PVE, "\n")
## without prior (uniform prior)
fitted_susie_noPrior <- susie(Xhat, yhat,
L = L_Susie,
prior_weights = NULL,
estimate_residual_variance = TRUE,
estimate_prior_variance = TRUE)
# susie_plot(fitted_susie_noPrior, y="PIP")
## RNA feature prior
prior <- read.table( paste0("/project2/xinhe/m6A/data_shared/m6A_QTLs/linear_model/annotation/RNAannot.prior/",iPeak,".prior"))
colnames(prior) <- c("SNP", "prior")
rownames(prior) <- prior$SNP
prior <- prior[colnames(Xhat),] ## match SNP order
if(nrow(prior) != ncol(Xhat) | any(colnames(Xhat) != prior$SNP)){cat("check prior for:", iPeak, "! \n")}
fitted_susie_RNAprior <- susie(Xhat, yhat,
L = L_Susie,
prior_weights = prior$prior,
estimate_residual_variance = TRUE,
estimate_prior_variance = TRUE)
# susie_plot(fitted_susie_RNAprior, y="PIP")
## Distance only prior
prior <- read.table( paste0("/project2/xinhe/m6A/data_shared/m6A_QTLs/linear_model/annotation/DistanceOnly.prior/",iPeak,".prior"))
colnames(prior) <- c("SNP", "prior")
rownames(prior) <- prior$SNP
prior <- prior[colnames(Xhat),] ## match SNP order
if(nrow(prior) != ncol(Xhat) | any(colnames(Xhat) != prior$SNP)){cat("check prior for:", iPeak, "! \n")}
fitted_susie_distPrior <- susie(Xhat, yhat,
L = L_Susie,
prior_weights = prior$prior,
estimate_residual_variance = TRUE,
estimate_prior_variance = TRUE)
# susie_plot(fitted_susie_distPrior, y="PIP")
SusieResult_peak = list(PEAK = as.character(iPeak),
z_score = z_score,
PVE = PVE,
fitted_susie_noPrior = fitted_susie_noPrior,
fitted_susie_RNAprior = fitted_susie_RNAprior,
fitted_susie_distPrior = fitted_susie_distPrior)
# saveRDS(SusieResult_peak, file = paste0(dir_susie_output, "/Peaks/", peak_map[as.character(iPeak),"peakName"],".RDS") )
SusieResult_peak
}
names(res) <- as.character(peak_map$peakID)
saveRDS(res, file = paste0(dir_susie_output, "/fitted_susie_allPeaks.RDS") )
saveRDS(peak_map, file = paste0(dir_susie_output, "/peak_map.RDS") )
cat("Done. \n")
### PVE for all peaks
PVE_all <- sapply(1:length(res), function(i) res[[i]]$PVE)
summary(PVE_all)m6AQTL_workflowr/analysis/combine_susie_results_allPeaks_L3.RmdsuppressMessages(library(dplyr))
suppressMessages(library(doParallel))
suppressMessages(library(susieR))
## Combine peak results into full statistics for all peak-SNP pairs
combine_susie_results <- function(fitted_susie_allPeaks, peakIDs, m6AQTL, m6A_phenotype, num_cores = 5){
registerDoParallel(num_cores)
SusieResult <- foreach( iPeak = peakIDs, .combine = rbind)%dopar%{
SusieResult_peak <- fitted_susie_allPeaks[[iPeak]]
Zscore <- SusieResult_peak$z_score
SNP <- names(Zscore)
if(susie_prior == "Susie_RNAfeature_prior"){
fitted_susie <- SusieResult_peak$fitted_susie_RNAprior
prior <- read.table( paste0("/project2/xinhe/m6A/data_shared/m6A_QTLs/linear_model/annotation/RNAannot.prior/",iPeak,".prior"))
colnames(prior) <- c("SNP", "prior")
rownames(prior) <- prior$SNP
prior <- prior[SNP, "prior"] ## match SNP order
}else if (susie_prior == "Susie_DistanceOnly_prior"){
fitted_susie <- SusieResult_peak$fitted_susie_distPrior
prior <- read.table( paste0("/project2/xinhe/m6A/data_shared/m6A_QTLs/linear_model/annotation/DistanceOnly.prior/",iPeak,".prior"))
colnames(prior) <- c("SNP", "prior")
rownames(prior) <- prior$SNP
prior <- prior[SNP, "prior"] ## match SNP order
}else if (susie_prior == "Susie_noPrior"){
fitted_susie <- SusieResult_peak$fitted_susie_noPrior
prior <- NA
}else{
stop("please choose the correct susie_prior!")
}
tmpOut <- data.frame(SNP = SNP, PEAK = iPeak, PIP = fitted_susie$pip, CS = NA, Zscore = Zscore, prior = prior)
idx_inCS <- unlist(fitted_susie$sets$cs)
if( !is.null(idx_inCS) ){ tmpOut$CS[idx_inCS] <- rep(names(fitted_susie$sets$cs),lapply(fitted_susie$sets$cs, length)) }
tmpOut[order(tmpOut$PIP,decreasing = T),]
# ## plot Susie result
# b <- rep(0,length(Zscore))
# b[which.max(abs(Zscore))] <- 1
# susie_plot(fitted_susie, y="PIP", b = b, max_cs=0.4, main = paste('SuSiE CS purity:', paste(round(fitted_susie$sets$purity[,1],2), collapse=';')))
}
matchId <- match(paste0(SusieResult$SNP,SusieResult$PEAK), paste0(m6AQTL$SNP,m6AQTL$PEAK))
SusieResult$beta <- m6AQTL$beta[matchId]
SusieResult$std.err <- m6AQTL$se[matchId]
SusieResult$sdY <- NA
for(iPeak in peakIDs){
tmpY <- m6A_phenotype[match(iPeak, m6A_phenotype$PEAK), -c(1:4)]
SusieResult$sdY[which(SusieResult$PEAK == iPeak )] <- sd(tmpY)
}
return(SusieResult)
}
## settings
m6A_version <- "jointPeak_threshold5_MeRIP_HISAT2Map"
m6A_phenotype_name <- "m6APeak_logOR_GC.IP.adjusted_qqnorm"
num_PCs <- 15
L_Susie <- 3
dir_fineMapping <- paste0("/project2/xinhe/m6A/data_shared/m6A_QTLs/", m6A_version, "/fineMapping/", m6A_phenotype_name, "/m6AQTL_allPeaks")
num_cores <- min(6, detectCores() - 1)
registerDoParallel(num_cores)
cat("m6A_version:", m6A_version, "\n")
cat("m6A_phenotype_name:", m6A_phenotype_name, "\n")
cat(num_PCs, "PCs included \n")
cat("Directory of SuSiE results:", dir_fineMapping, "\n")
cat("L_Susie =", L_Susie, "\n")
## Load m6A-QTL data
dir_m6AQTL <- paste0("/project2/xinhe/m6A/data_shared/m6A_QTLs/", m6A_version, "/fastQTL_cis_100kb/", m6A_phenotype_name)
m6AQTL <- readRDS(paste0(dir_m6AQTL, "/m6AQTL.", m6A_phenotype_name, ".", num_PCs, "PCs.fastQTL.nominals.rds"))
m6AQTLs_PermPass <- readRDS(paste0(dir_m6AQTL, "/m6AQTL.", m6A_phenotype_name, ".", num_PCs, "PCs.fastQTL.permutations.rds"))
m6AQTLs_PermPass <- m6AQTLs_PermPass[!is.na(m6AQTLs_PermPass$qvalue), ]
dir_m6A_phenotype <- paste0("/project2/xinhe/m6A/m6A_seq/m6A_QTL/peakcalling/", m6A_version, "/", m6A_phenotype_name)
m6A_phenotype <- read.table(paste0(dir_m6A_phenotype, "/", m6A_phenotype_name, ".fastQTL.txt.gz"), header = T, comment.char = "", stringsAsFactors = F)
colnames(m6A_phenotype)[1] <- "chr"
## Susie estimated prior variance
### load Susie results for all peaks
susie_setting <- paste0("Susie_L", L_Susie, "_estPriorVar")
print(susie_setting)
dir_susie_output <- paste0(dir_fineMapping, "/", susie_setting)
peak_map <- readRDS(paste0(dir_susie_output, "/peak_map.RDS") )
fitted_susie_allPeaks <- readRDS(paste0(dir_susie_output, "/fitted_susie_allPeaks.RDS") )
if(nrow(peak_map) != length(fitted_susie_allPeaks)){
stop("unequal number of peaks \n")
}
### Combine peak results into full statistics for all peak-SNP pairs
#### without prior (uniform prior)
susie_prior <- "Susie_noPrior"
print(susie_prior)
dir_output <- paste0(dir_fineMapping, "/", susie_setting, "/", susie_prior)
dir.create(dir_output, showWarnings = F, recursive = T)
for(thresh_qvalue_set in c(0.1, 0.2, 0.3)){
ePeaks_set <- as.character(m6AQTLs_PermPass[m6AQTLs_PermPass$qvalue < thresh_qvalue_set, "PEAK"])
cat(length(ePeaks_set), "peaks with qvalue < ", thresh_qvalue_set, "\n")
SusieResult <- combine_susie_results(fitted_susie_allPeaks, ePeaks_set, m6AQTL, m6A_phenotype, num_cores = num_cores)
cat(length(unique(SusieResult[!is.na(SusieResult$CS), "PEAK"])), "peaks have CS \n")
saveRDS(SusieResult, file = paste0(dir_output, "/SusieResult_qvalue_", thresh_qvalue_set,".RDS"))
}
#### RNA feature prior
susie_prior <- "Susie_RNAfeature_prior"
print(susie_prior)
dir_output <- paste0(dir_fineMapping, "/", susie_setting, "/", susie_prior)
dir.create(dir_output, showWarnings = F, recursive = T)
for(thresh_qvalue_set in c(0.1, 0.2, 0.3)){
ePeaks_set <- as.character(m6AQTLs_PermPass[m6AQTLs_PermPass$qvalue < thresh_qvalue_set, "PEAK"])
cat(length(ePeaks_set), "peaks with qvalue < ", thresh_qvalue_set, "\n")
SusieResult <- combine_susie_results(fitted_susie_allPeaks, ePeaks_set, m6AQTL, m6A_phenotype, num_cores = num_cores)
cat(length(unique(SusieResult[!is.na(SusieResult$CS), "PEAK"])), "peaks have CS \n")
saveRDS(SusieResult, file = paste0(dir_output, "/SusieResult_qvalue_", thresh_qvalue_set,".RDS"))
}
#### Distance only prior
susie_prior <- "Susie_DistanceOnly_prior"
print(susie_prior)
dir_output <- paste0(dir_fineMapping, "/", susie_setting, "/", susie_prior)
dir.create(dir_output, showWarnings = F, recursive = T)
for(thresh_qvalue_set in c(0.1, 0.2, 0.3)){
ePeaks_set <- as.character(m6AQTLs_PermPass[m6AQTLs_PermPass$qvalue < thresh_qvalue_set, "PEAK"])
cat(length(ePeaks_set), "peaks with qvalue < ", thresh_qvalue_set, "\n")
SusieResult <- combine_susie_results(fitted_susie_allPeaks, ePeaks_set, m6AQTL, m6A_phenotype, num_cores = num_cores)
cat(length(unique(SusieResult[!is.na(SusieResult$CS), "PEAK"])), "peaks have CS \n")
saveRDS(SusieResult, file = paste0(dir_output, "/SusieResult_qvalue_", thresh_qvalue_set,".RDS"))
}
sessionInfo()R version 3.5.1 (2018-07-02)
Platform: x86_64-pc-linux-gnu (64-bit)
Running under: Scientific Linux 7.4 (Nitrogen)
Matrix products: default
BLAS/LAPACK: /software/openblas-0.2.19-el7-x86_64/lib/libopenblas_haswellp-r0.2.19.so
locale:
[1] LC_CTYPE=en_US.UTF-8 LC_NUMERIC=C
[3] LC_TIME=en_US.UTF-8 LC_COLLATE=en_US.UTF-8
[5] LC_MONETARY=en_US.UTF-8 LC_MESSAGES=en_US.UTF-8
[7] LC_PAPER=en_US.UTF-8 LC_NAME=C
[9] LC_ADDRESS=C LC_TELEPHONE=C
[11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] workflowr_1.6.0
loaded via a namespace (and not attached):
[1] Rcpp_1.0.3 rprojroot_1.3-2 digest_0.6.23 later_1.0.0
[5] R6_2.4.1 backports_1.1.5 git2r_0.26.1.9000 magrittr_1.5
[9] evaluate_0.14 stringi_1.4.5 rlang_0.4.4 fs_1.3.1
[13] promises_1.1.0 whisker_0.4 rmarkdown_2.1 tools_3.5.1
[17] stringr_1.4.0 glue_1.3.1 httpuv_1.5.2 xfun_0.12
[21] yaml_2.2.0 compiler_3.5.1 htmltools_0.4.0 knitr_1.28