Last updated: 2021-09-09
Checks: 6 1
Knit directory: ctwas_applied/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish
to commit the R Markdown file and build the HTML.
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(20210726)
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 results in this page were generated with repository version 59e5f4d. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
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:
Unstaged changes:
Modified: analysis/ukb-d-30500_irnt_Liver.Rmd
Modified: analysis/ukb-d-30500_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30600_irnt_Liver.Rmd
Modified: analysis/ukb-d-30600_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30610_irnt_Liver.Rmd
Modified: analysis/ukb-d-30610_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30620_irnt_Liver.Rmd
Modified: analysis/ukb-d-30620_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30630_irnt_Liver.Rmd
Modified: analysis/ukb-d-30630_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30640_irnt_Liver.Rmd
Modified: analysis/ukb-d-30640_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30650_irnt_Liver.Rmd
Modified: analysis/ukb-d-30650_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30660_irnt_Liver.Rmd
Modified: analysis/ukb-d-30660_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30670_irnt_Liver.Rmd
Modified: analysis/ukb-d-30670_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30680_irnt_Liver.Rmd
Modified: analysis/ukb-d-30690_irnt_Liver.Rmd
Modified: analysis/ukb-d-30690_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30700_irnt_Liver.Rmd
Modified: analysis/ukb-d-30700_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30710_irnt_Liver.Rmd
Modified: analysis/ukb-d-30710_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30720_irnt_Liver.Rmd
Modified: analysis/ukb-d-30720_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30730_irnt_Liver.Rmd
Modified: analysis/ukb-d-30740_irnt_Liver.Rmd
Modified: analysis/ukb-d-30740_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30750_irnt_Liver.Rmd
Modified: analysis/ukb-d-30750_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30760_irnt_Liver.Rmd
Modified: analysis/ukb-d-30760_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30770_irnt_Liver.Rmd
Modified: analysis/ukb-d-30770_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30780_irnt_Liver.Rmd
Modified: analysis/ukb-d-30780_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30790_irnt_Liver.Rmd
Modified: analysis/ukb-d-30800_irnt_Liver.Rmd
Modified: analysis/ukb-d-30800_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30810_irnt_Liver.Rmd
Modified: analysis/ukb-d-30820_irnt_Liver.Rmd
Modified: analysis/ukb-d-30820_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30830_irnt_Liver.Rmd
Modified: analysis/ukb-d-30830_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30840_irnt_Liver.Rmd
Modified: analysis/ukb-d-30840_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30850_irnt_Liver.Rmd
Modified: analysis/ukb-d-30850_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30860_irnt_Liver.Rmd
Modified: analysis/ukb-d-30860_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30870_irnt_Liver.Rmd
Modified: analysis/ukb-d-30870_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30880_irnt_Liver.Rmd
Modified: analysis/ukb-d-30880_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30890_irnt_Liver.Rmd
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 repository in which changes were made to the R Markdown (analysis/ukb-d-30880_irnt_Whole_Blood.Rmd
) and HTML (docs/ukb-d-30880_irnt_Whole_Blood.html
) files. If you’ve configured a remote Git repository (see ?wflow_git_remote
), click on the hyperlinks in the table below to view the files as they were in that past version.
File | Version | Author | Date | Message |
---|---|---|---|---|
Rmd | cbf7408 | wesleycrouse | 2021-09-08 | adding enrichment to reports |
html | cbf7408 | wesleycrouse | 2021-09-08 | adding enrichment to reports |
Rmd | 4970e3e | wesleycrouse | 2021-09-08 | updating reports |
html | 4970e3e | wesleycrouse | 2021-09-08 | updating reports |
Rmd | dfd2b5f | wesleycrouse | 2021-09-07 | regenerating reports |
html | dfd2b5f | wesleycrouse | 2021-09-07 | regenerating reports |
Rmd | 61b53b3 | wesleycrouse | 2021-09-06 | updated PVE calculation |
html | 61b53b3 | wesleycrouse | 2021-09-06 | updated PVE calculation |
Rmd | 837dd01 | wesleycrouse | 2021-09-01 | adding additional fixedsigma report |
Rmd | d0a5417 | wesleycrouse | 2021-08-30 | adding new reports to the index |
Rmd | 0922de7 | wesleycrouse | 2021-08-18 | updating all reports with locus plots |
html | 0922de7 | wesleycrouse | 2021-08-18 | updating all reports with locus plots |
html | 1c62980 | wesleycrouse | 2021-08-11 | Updating reports |
Rmd | 5981e80 | wesleycrouse | 2021-08-11 | Adding more reports |
html | 5981e80 | wesleycrouse | 2021-08-11 | Adding more reports |
Rmd | da9f015 | wesleycrouse | 2021-08-07 | adding more reports |
html | da9f015 | wesleycrouse | 2021-08-07 | adding more reports |
These are the results of a ctwas
analysis of the UK Biobank trait Urate (quantile)
using Whole_Blood
gene weights.
The GWAS was conducted by the Neale Lab, and the biomarker traits we analyzed are discussed here. Summary statistics were obtained from IEU OpenGWAS using GWAS ID: ukb-d-30880_irnt
. Results were obtained from from IEU rather than Neale Lab because they are in a standardard format (GWAS VCF). Note that 3 of the 34 biomarker traits were not available from IEU and were excluded from analysis.
The weights are mashr GTEx v8 models on Whole_Blood
eQTL obtained from PredictDB. We performed a full harmonization of the variants, including recovering strand ambiguous variants. This procedure is discussed in a separate document. (TO-DO: add report that describes harmonization)
LD matrices were computed from a 10% subset of Neale lab subjects. Subjects were matched using the plate and well information from genotyping. We included only biallelic variants with MAF>0.01 in the original Neale Lab GWAS. (TO-DO: add more details [number of subjects, variants, etc])
TO-DO: add enhanced QC reporting (total number of weights, why each variant was missing for all genes)
qclist_all <- list()
qc_files <- paste0(results_dir, "/", list.files(results_dir, pattern="exprqc.Rd"))
for (i in 1:length(qc_files)){
load(qc_files[i])
chr <- unlist(strsplit(rev(unlist(strsplit(qc_files[i], "_")))[1], "[.]"))[1]
qclist_all[[chr]] <- cbind(do.call(rbind, lapply(qclist,unlist)), as.numeric(substring(chr,4)))
}
qclist_all <- data.frame(do.call(rbind, qclist_all))
colnames(qclist_all)[ncol(qclist_all)] <- "chr"
rm(qclist, wgtlist, z_gene_chr)
#number of imputed weights
nrow(qclist_all)
[1] 11095
#number of imputed weights by chromosome
table(qclist_all$chr)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
1129 747 624 400 479 621 560 383 404 430 682 652 192 362 331
16 17 18 19 20 21 22
551 725 159 911 313 130 310
#proportion of imputed weights without missing variants
mean(qclist_all$nmiss==0)
[1] 0.762776
#load ctwas results
ctwas_res <- data.table::fread(paste0(results_dir, "/", analysis_id, "_ctwas.susieIrss.txt"))
#make unique identifier for regions
ctwas_res$region_tag <- paste(ctwas_res$region_tag1, ctwas_res$region_tag2, sep="_")
#compute PVE for each gene/SNP
ctwas_res$PVE = ctwas_res$susie_pip*ctwas_res$mu2/sample_size #check PVE calculation
#separate gene and SNP results
ctwas_gene_res <- ctwas_res[ctwas_res$type == "gene", ]
ctwas_gene_res <- data.frame(ctwas_gene_res)
ctwas_snp_res <- ctwas_res[ctwas_res$type == "SNP", ]
ctwas_snp_res <- data.frame(ctwas_snp_res)
#add gene information to results
sqlite <- RSQLite::dbDriver("SQLite")
db = RSQLite::dbConnect(sqlite, paste0("/project2/compbio/predictdb/mashr_models/mashr_", weight, ".db"))
query <- function(...) RSQLite::dbGetQuery(db, ...)
gene_info <- query("select gene, genename from extra")
gene_info <- query("select gene, genename, gene_type from extra")
RSQLite::dbDisconnect(db)
ctwas_gene_res <- cbind(ctwas_gene_res, gene_info[sapply(ctwas_gene_res$id, match, gene_info$gene), c("genename", "gene_type")])
#add z score to results
load(paste0(results_dir, "/", analysis_id, "_expr_z_gene.Rd"))
ctwas_gene_res$z <- z_gene[ctwas_gene_res$id,]$z
#load(paste0(results_dir, "/", analysis_id, "_expr_z_snp.Rd")) #for new version, stored after harmonization
z_snp <- readRDS(paste0(results_dir, "/", trait_id, ".RDS")) #for old version, unharmonized
z_snp <- z_snp[z_snp$id %in% ctwas_snp_res$id,] #subset snps to those included in analysis, note some are duplicated, need to match which allele was used
ctwas_snp_res$z <- z_snp$z[match(ctwas_snp_res$id, z_snp$id)] #for duplicated snps, this arbitrarily uses the first allele
ctwas_snp_res$z_flag <- as.numeric(ctwas_snp_res$id %in% z_snp$id[duplicated(z_snp$id)]) #mark the unclear z scores, flag=1
#formatting and rounding for tables
ctwas_gene_res$z <- round(ctwas_gene_res$z,2)
ctwas_snp_res$z <- round(ctwas_snp_res$z,2)
ctwas_gene_res$susie_pip <- round(ctwas_gene_res$susie_pip,3)
ctwas_snp_res$susie_pip <- round(ctwas_snp_res$susie_pip,3)
ctwas_gene_res$mu2 <- round(ctwas_gene_res$mu2,2)
ctwas_snp_res$mu2 <- round(ctwas_snp_res$mu2,2)
ctwas_gene_res$PVE <- signif(ctwas_gene_res$PVE, 2)
ctwas_snp_res$PVE <- signif(ctwas_snp_res$PVE, 2)
#merge gene and snp results with added information
ctwas_gene_res$z_flag=NA
ctwas_snp_res$genename=NA
ctwas_snp_res$gene_type=NA
ctwas_res <- rbind(ctwas_gene_res,
ctwas_snp_res[,colnames(ctwas_gene_res)])
#store columns to report
report_cols <- colnames(ctwas_gene_res)[!(colnames(ctwas_gene_res) %in% c("type", "region_tag1", "region_tag2", "cs_index", "gene_type", "z_flag", "id", "chrom", "pos"))]
first_cols <- c("genename", "region_tag")
report_cols <- c(first_cols, report_cols[!(report_cols %in% first_cols)])
report_cols_snps <- c("id", report_cols[-1])
#get number of SNPs from s1 results; adjust for thin
ctwas_res_s1 <- data.table::fread(paste0(results_dir, "/", analysis_id, "_ctwas.s1.susieIrss.txt"))
n_snps <- sum(ctwas_res_s1$type=="SNP")/thin
rm(ctwas_res_s1)
library(ggplot2)
library(cowplot)
********************************************************
Note: As of version 1.0.0, cowplot does not change the
default ggplot2 theme anymore. To recover the previous
behavior, execute:
theme_set(theme_cowplot())
********************************************************
load(paste0(results_dir, "/", analysis_id, "_ctwas.s2.susieIrssres.Rd"))
df <- data.frame(niter = rep(1:ncol(group_prior_rec), 2),
value = c(group_prior_rec[1,], group_prior_rec[2,]),
group = rep(c("Gene", "SNP"), each = ncol(group_prior_rec)))
df$group <- as.factor(df$group)
df$value[df$group=="SNP"] <- df$value[df$group=="SNP"]*thin #adjust parameter to account for thin argument
p_pi <- ggplot(df, aes(x=niter, y=value, group=group)) +
geom_line(aes(color=group)) +
geom_point(aes(color=group)) +
xlab("Iteration") + ylab(bquote(pi)) +
ggtitle("Prior mean") +
theme_cowplot()
df <- data.frame(niter = rep(1:ncol(group_prior_var_rec), 2),
value = c(group_prior_var_rec[1,], group_prior_var_rec[2,]),
group = rep(c("Gene", "SNP"), each = ncol(group_prior_var_rec)))
df$group <- as.factor(df$group)
p_sigma2 <- ggplot(df, aes(x=niter, y=value, group=group)) +
geom_line(aes(color=group)) +
geom_point(aes(color=group)) +
xlab("Iteration") + ylab(bquote(sigma^2)) +
ggtitle("Prior variance") +
theme_cowplot()
plot_grid(p_pi, p_sigma2)
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
#estimated group prior
estimated_group_prior <- group_prior_rec[,ncol(group_prior_rec)]
names(estimated_group_prior) <- c("gene", "snp")
estimated_group_prior["snp"] <- estimated_group_prior["snp"]*thin #adjust parameter to account for thin argument
print(estimated_group_prior)
gene snp
0.008022976 0.000204890
#estimated group prior variance
estimated_group_prior_var <- group_prior_var_rec[,ncol(group_prior_var_rec)]
names(estimated_group_prior_var) <- c("gene", "snp")
print(estimated_group_prior_var)
gene snp
20.81243 18.20088
#report sample size
print(sample_size)
[1] 343836
#report group size
group_size <- c(nrow(ctwas_gene_res), n_snps)
print(group_size)
[1] 11095 8697330
#estimated group PVE
estimated_group_pve <- estimated_group_prior_var*estimated_group_prior*group_size/sample_size #check PVE calculation
names(estimated_group_pve) <- c("gene", "snp")
print(estimated_group_pve)
gene snp
0.005388083 0.094329551
#compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))
[1] 0.01825411 0.71407143
#distribution of PIPs
hist(ctwas_gene_res$susie_pip, xlim=c(0,1), main="Distribution of Gene PIPs")
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
#genes with PIP>0.8 or 20 highest PIPs
head(ctwas_gene_res[order(-ctwas_gene_res$susie_pip),report_cols], max(sum(ctwas_gene_res$susie_pip>0.8), 20))
genename region_tag susie_pip mu2 PVE z
10765 ZDHHC18 1_18 1.000 142.25 4.1e-04 12.20
9813 MUC1 1_77 1.000 246.05 7.2e-04 -18.55
3752 KCNK17 6_30 0.971 39.18 1.1e-04 -6.59
10816 MBD5 2_88 0.969 81.94 2.3e-04 9.34
6598 NRG1 8_31 0.966 71.88 2.0e-04 -8.59
5012 TRIM29 11_72 0.941 47.36 1.3e-04 7.00
8523 ZNF217 20_31 0.932 25.40 6.9e-05 -4.69
5874 ANO7 2_143 0.913 37.89 1.0e-04 -6.01
6726 ZNF276 16_54 0.894 31.37 8.2e-05 -5.44
9820 TLCD2 17_2 0.891 25.23 6.5e-05 -4.79
8177 THBS3 1_77 0.890 216.19 5.6e-04 16.74
4484 CSNK1G2 19_2 0.859 27.97 7.0e-05 5.00
10355 SLC39A10 2_116 0.846 27.35 6.7e-05 -4.94
4505 AGAP3 7_94 0.823 25.21 6.0e-05 -4.56
9325 ZNF816 19_37 0.803 33.38 7.8e-05 -5.73
920 STXBP2 19_7 0.798 23.95 5.6e-05 4.38
1916 TRIM35 8_27 0.772 24.88 5.6e-05 4.74
2540 NECTIN1 11_72 0.753 22.58 4.9e-05 -4.47
5564 ATP1B1 1_82 0.751 24.21 5.3e-05 4.27
12181 EGLN2 19_30 0.746 24.64 5.3e-05 -4.12
#plot PIP vs effect size
plot(ctwas_gene_res$susie_pip, ctwas_gene_res$mu2, xlab="PIP", ylab="mu^2", main="Gene PIPs vs Effect Size")
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
#genes with 20 largest effect sizes
head(ctwas_gene_res[order(-ctwas_gene_res$mu2),report_cols],20)
genename region_tag susie_pip mu2 PVE z
6501 PDIA4 7_92 0.000 7907.58 0.0e+00 2.15
168 SPRTN 1_118 0.000 5408.85 1.6e-10 3.40
3138 EXOC8 1_118 0.000 3884.27 0.0e+00 2.52
4687 TMEM60 7_49 0.000 2850.57 0.0e+00 2.04
9157 ZNF518B 4_11 0.000 1433.04 0.0e+00 58.12
2475 SLC2A9 4_11 0.000 1391.98 0.0e+00 53.26
10421 ZNF786 7_92 0.000 1280.98 0.0e+00 -1.75
8305 ZNF282 7_92 0.000 1229.89 0.0e+00 -1.63
10367 ZNF398 7_92 0.000 948.52 0.0e+00 -1.27
3140 TSNAX 1_118 0.000 758.59 0.0e+00 1.97
11094 APTR 7_49 0.000 494.47 0.0e+00 -2.79
7145 DISC1 1_118 0.000 459.53 0.0e+00 -1.78
98 PHTF2 7_49 0.000 415.54 0.0e+00 0.42
10882 ZNF425 7_92 0.000 381.88 0.0e+00 -0.60
2737 TRIM38 6_20 0.005 315.61 4.5e-06 -23.44
9813 MUC1 1_77 1.000 246.05 7.2e-04 -18.55
10864 SPDYC 11_36 0.001 230.89 9.2e-07 -9.34
8177 THBS3 1_77 0.890 216.19 5.6e-04 16.74
2953 NRBP1 2_16 0.021 215.65 1.3e-05 15.35
11023 SIPA1 11_36 0.000 214.21 8.6e-11 12.60
#genes with 20 highest pve
head(ctwas_gene_res[order(-ctwas_gene_res$PVE),report_cols],20)
genename region_tag susie_pip mu2 PVE z
9813 MUC1 1_77 1.000 246.05 7.2e-04 -18.55
8177 THBS3 1_77 0.890 216.19 5.6e-04 16.74
10765 ZDHHC18 1_18 1.000 142.25 4.1e-04 12.20
10816 MBD5 2_88 0.969 81.94 2.3e-04 9.34
6598 NRG1 8_31 0.966 71.88 2.0e-04 -8.59
5012 TRIM29 11_72 0.941 47.36 1.3e-04 7.00
3752 KCNK17 6_30 0.971 39.18 1.1e-04 -6.59
6089 FADS1 11_34 0.679 52.66 1.0e-04 7.27
5874 ANO7 2_143 0.913 37.89 1.0e-04 -6.01
171 UQCRC1 3_34 0.687 43.25 8.6e-05 -6.52
3544 FAM35A 10_57 0.740 39.04 8.4e-05 -6.62
6726 ZNF276 16_54 0.894 31.37 8.2e-05 -5.44
9325 ZNF816 19_37 0.803 33.38 7.8e-05 -5.73
12389 RP11-686O6.2 2_120 0.359 71.39 7.5e-05 8.04
10021 ZKSCAN4 6_22 0.155 154.44 7.0e-05 12.61
4484 CSNK1G2 19_2 0.859 27.97 7.0e-05 5.00
8523 ZNF217 20_31 0.932 25.40 6.9e-05 -4.69
10355 SLC39A10 2_116 0.846 27.35 6.7e-05 -4.94
4126 QRICH2 17_42 0.671 34.09 6.7e-05 6.55
9820 TLCD2 17_2 0.891 25.23 6.5e-05 -4.79
#genes with 20 largest z scores
head(ctwas_gene_res[order(-abs(ctwas_gene_res$z)),report_cols],20)
genename region_tag susie_pip mu2 PVE z
9157 ZNF518B 4_11 0.000 1433.04 0.0e+00 58.12
2475 SLC2A9 4_11 0.000 1391.98 0.0e+00 53.26
2737 TRIM38 6_20 0.005 315.61 4.5e-06 -23.44
9813 MUC1 1_77 1.000 246.05 7.2e-04 -18.55
8177 THBS3 1_77 0.890 216.19 5.6e-04 16.74
187 HFE 6_20 0.005 112.31 1.7e-06 -16.15
12299 U91328.19 6_20 0.007 112.39 2.4e-06 15.64
8691 MAP3K11 11_36 0.000 166.15 3.7e-15 15.38
2953 NRBP1 2_16 0.021 215.65 1.3e-05 15.35
2956 SNX17 2_16 0.030 196.19 1.7e-05 14.39
9933 BTN3A2 6_20 0.010 88.78 2.7e-06 -13.24
12027 CTA-14H9.5 6_20 0.016 91.31 4.3e-06 -13.00
9350 HIST1H2BC 6_20 0.005 70.83 1.1e-06 -12.78
10021 ZKSCAN4 6_22 0.155 154.44 7.0e-05 12.61
11023 SIPA1 11_36 0.000 214.21 8.6e-11 12.60
6712 ZSCAN12 6_22 0.104 154.07 4.7e-05 -12.57
10765 ZDHHC18 1_18 1.000 142.25 4.1e-04 12.20
12301 RP1-86C11.7 6_21 0.000 200.97 1.7e-07 -12.06
10762 ASAH2B 10_33 0.015 128.27 5.4e-06 11.56
9293 R3HDM2 12_36 0.040 114.76 1.3e-05 -11.33
#set nominal signifiance threshold for z scores
alpha <- 0.05
#bonferroni adjusted threshold for z scores
sig_thresh <- qnorm(1-(alpha/nrow(ctwas_gene_res)/2), lower=T)
#Q-Q plot for z scores
obs_z <- ctwas_gene_res$z[order(ctwas_gene_res$z)]
exp_z <- qnorm((1:nrow(ctwas_gene_res))/nrow(ctwas_gene_res))
plot(exp_z, obs_z, xlab="Expected z", ylab="Observed z", main="Gene z score Q-Q plot")
abline(a=0,b=1)
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
#plot z score vs PIP
plot(abs(ctwas_gene_res$z), ctwas_gene_res$susie_pip, xlab="abs(z)", ylab="PIP")
abline(v=sig_thresh, col="red", lty=2)
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
#proportion of significant z scores
mean(abs(ctwas_gene_res$z) > sig_thresh)
[1] 0.02631816
#genes with most significant z scores
head(ctwas_gene_res[order(-abs(ctwas_gene_res$z)),report_cols],20)
genename region_tag susie_pip mu2 PVE z
9157 ZNF518B 4_11 0.000 1433.04 0.0e+00 58.12
2475 SLC2A9 4_11 0.000 1391.98 0.0e+00 53.26
2737 TRIM38 6_20 0.005 315.61 4.5e-06 -23.44
9813 MUC1 1_77 1.000 246.05 7.2e-04 -18.55
8177 THBS3 1_77 0.890 216.19 5.6e-04 16.74
187 HFE 6_20 0.005 112.31 1.7e-06 -16.15
12299 U91328.19 6_20 0.007 112.39 2.4e-06 15.64
8691 MAP3K11 11_36 0.000 166.15 3.7e-15 15.38
2953 NRBP1 2_16 0.021 215.65 1.3e-05 15.35
2956 SNX17 2_16 0.030 196.19 1.7e-05 14.39
9933 BTN3A2 6_20 0.010 88.78 2.7e-06 -13.24
12027 CTA-14H9.5 6_20 0.016 91.31 4.3e-06 -13.00
9350 HIST1H2BC 6_20 0.005 70.83 1.1e-06 -12.78
10021 ZKSCAN4 6_22 0.155 154.44 7.0e-05 12.61
11023 SIPA1 11_36 0.000 214.21 8.6e-11 12.60
6712 ZSCAN12 6_22 0.104 154.07 4.7e-05 -12.57
10765 ZDHHC18 1_18 1.000 142.25 4.1e-04 12.20
12301 RP1-86C11.7 6_21 0.000 200.97 1.7e-07 -12.06
10762 ASAH2B 10_33 0.015 128.27 5.4e-06 11.56
9293 R3HDM2 12_36 0.040 114.76 1.3e-05 -11.33
ctwas_gene_res_sortz <- ctwas_gene_res[order(-abs(ctwas_gene_res$z)),]
n_plots <- 5
for (region_tag_plot in head(unique(ctwas_gene_res_sortz$region_tag), n_plots)){
ctwas_res_region <- ctwas_res[ctwas_res$region_tag==region_tag_plot,]
start <- min(ctwas_res_region$pos)
end <- max(ctwas_res_region$pos)
ctwas_res_region <- ctwas_res_region[order(ctwas_res_region$pos),]
ctwas_res_region_gene <- ctwas_res_region[ctwas_res_region$type=="gene",]
ctwas_res_region_snp <- ctwas_res_region[ctwas_res_region$type=="SNP",]
#region name
print(paste0("Region: ", region_tag_plot))
#table of genes in region
print(ctwas_res_region_gene[,report_cols])
par(mfrow=c(4,1))
#gene z scores
plot(ctwas_res_region_gene$pos, abs(ctwas_res_region_gene$z), xlab="Position", ylab="abs(gene_z)", xlim=c(start,end),
ylim=c(0,max(sig_thresh, abs(ctwas_res_region_gene$z))),
main=paste0("Region: ", region_tag_plot))
abline(h=sig_thresh,col="red",lty=2)
#significance threshold for SNPs
alpha_snp <- 5*10^(-8)
sig_thresh_snp <- qnorm(1-alpha_snp/2, lower=T)
#snp z scores
plot(ctwas_res_region_snp$pos, abs(ctwas_res_region_snp$z), xlab="Position", ylab="abs(snp_z)",xlim=c(start,end),
ylim=c(0,max(sig_thresh_snp, max(abs(ctwas_res_region_snp$z)))))
abline(h=sig_thresh_snp,col="purple",lty=2)
#gene pips
plot(ctwas_res_region_gene$pos, ctwas_res_region_gene$susie_pip, xlab="Position", ylab="Gene PIP", xlim=c(start,end), ylim=c(0,1))
abline(h=0.8,col="blue",lty=2)
#snp pips
plot(ctwas_res_region_snp$pos, ctwas_res_region_snp$susie_pip, xlab="Position", ylab="SNP PIP", xlim=c(start,end), ylim=c(0,1))
abline(h=0.8,col="blue",lty=2)
}
[1] "Region: 4_11"
genename region_tag susie_pip mu2 PVE z
2475 SLC2A9 4_11 0 1391.98 0 53.26
9157 ZNF518B 4_11 0 1433.04 0 58.12
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 6_20"
genename region_tag susie_pip mu2 PVE z
2737 TRIM38 6_20 0.005 315.61 4.5e-06 -23.44
12299 U91328.19 6_20 0.007 112.39 2.4e-06 15.64
12332 U91328.22 6_20 0.007 13.56 2.7e-07 1.70
12483 HIST1H3A 6_20 0.005 57.42 8.6e-07 -8.42
10043 HIST1H1C 6_20 0.008 39.66 9.5e-07 4.31
187 HFE 6_20 0.005 112.31 1.7e-06 -16.15
10373 HIST1H4C 6_20 0.012 15.23 5.3e-07 -1.83
10005 HIST1H1T 6_20 0.008 84.91 1.9e-06 11.23
9350 HIST1H2BC 6_20 0.005 70.83 1.1e-06 -12.78
9348 HIST1H2AC 6_20 0.005 52.10 8.0e-07 -10.66
6686 HIST1H2BD 6_20 0.007 33.25 7.2e-07 -5.23
12425 HIST1H2BE 6_20 0.016 30.38 1.4e-06 3.69
12528 HIST1H2BF 6_20 0.005 6.80 9.8e-08 -1.88
10342 HIST1H2AD 6_20 0.011 10.83 3.4e-07 -0.12
12518 HIST1H4E 6_20 0.010 10.30 3.0e-07 0.02
12520 HIST1H2AE 6_20 0.005 6.42 1.0e-07 0.04
12444 HIST1H3E 6_20 0.021 21.39 1.3e-06 2.68
12482 HIST1H2BH 6_20 0.005 13.99 2.0e-07 -1.79
12409 HIST1H3G 6_20 0.007 19.06 3.9e-07 -5.30
6688 HIST1H4H 6_20 0.005 11.12 1.8e-07 2.02
2669 BTN3A3 6_20 0.005 6.39 9.0e-08 0.57
9933 BTN3A2 6_20 0.010 88.78 2.7e-06 -13.24
3723 BTN2A2 6_20 0.084 51.74 1.3e-05 7.22
308 BTN3A1 6_20 0.005 8.04 1.2e-07 -3.09
2771 BTN2A1 6_20 0.009 12.36 3.1e-07 -1.72
12027 CTA-14H9.5 6_20 0.016 91.31 4.3e-06 -13.00
9548 HMGN4 6_20 0.009 9.61 2.4e-07 0.93
11308 HCG11 6_20 0.008 8.72 2.0e-07 -1.60
5867 ABT1 6_20 0.056 30.31 4.9e-06 4.09
9407 ZNF322 6_20 0.007 17.98 3.5e-07 5.16
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 1_77"
genename region_tag susie_pip mu2 PVE z
7201 PMVK 1_77 0.006 33.56 5.9e-07 4.16
7202 PBXIP1 1_77 0.006 33.05 5.6e-07 4.11
6909 SHC1 1_77 0.001 5.76 1.1e-08 0.78
8673 CKS1B 1_77 0.001 5.69 1.1e-08 -0.73
6907 ZBTB7B 1_77 0.002 22.16 1.3e-07 -3.16
7204 DCST2 1_77 0.058 64.79 1.1e-05 7.14
7205 DCST1 1_77 0.058 64.79 1.1e-05 -7.14
5634 ADAM15 1_77 0.001 19.76 7.6e-08 -0.59
5644 EFNA3 1_77 0.001 6.76 1.3e-08 1.16
8179 EFNA1 1_77 0.002 37.49 2.5e-07 4.74
8178 SLC50A1 1_77 0.001 7.25 1.5e-08 -1.62
8670 MTX1 1_77 0.001 81.54 1.7e-07 -9.72
9813 MUC1 1_77 1.000 246.05 7.2e-04 -18.55
8177 THBS3 1_77 0.890 216.19 5.6e-04 16.74
9103 GBA 1_77 0.001 14.72 3.6e-08 3.90
6918 FAM189B 1_77 0.002 11.67 5.6e-08 -0.20
5649 HCN3 1_77 0.001 68.05 2.9e-07 1.50
6917 FDPS 1_77 0.002 72.97 4.5e-07 0.34
4425 DAP3 1_77 0.001 8.78 1.6e-08 0.93
7207 YY1AP1 1_77 0.001 38.40 1.4e-07 5.21
4434 SYT11 1_77 0.087 90.55 2.3e-05 -1.33
5648 RIT1 1_77 0.001 25.86 5.1e-08 -1.76
4426 KIAA0907 1_77 0.002 22.14 1.1e-07 -2.57
3099 ARHGEF2 1_77 0.001 62.77 1.7e-07 -1.56
7227 SSR2 1_77 0.001 11.61 2.7e-08 -2.36
3100 LAMTOR2 1_77 0.001 7.44 1.8e-08 -0.72
4430 RAB25 1_77 0.001 23.70 4.4e-08 0.87
10224 SEMA4A 1_77 0.001 12.77 2.5e-08 -3.38
6920 PMF1 1_77 0.008 34.91 7.7e-07 3.77
11586 BGLAP 1_77 0.009 28.86 7.4e-07 -2.38
6919 PAQR6 1_77 0.003 19.18 1.7e-07 -1.66
7226 TMEM79 1_77 0.003 20.70 1.9e-07 -2.09
10714 SMG5 1_77 0.003 20.03 1.7e-07 -2.06
10641 GLMP 1_77 0.003 19.54 1.6e-07 -2.03
7224 TSACC 1_77 0.001 5.49 1.1e-08 0.50
7225 CCT3 1_77 0.001 4.92 9.4e-09 -0.21
3101 MEF2D 1_77 0.001 8.93 2.4e-08 -0.77
9652 IQGAP3 1_77 0.001 5.08 9.7e-09 0.20
10047 TTC24 1_77 0.006 22.57 3.9e-07 -2.27
7210 NAXE 1_77 0.001 11.00 4.0e-08 1.23
4428 BCAN 1_77 0.001 8.55 2.5e-08 -1.03
5586 RRNAD1 1_77 0.001 6.48 1.5e-08 -0.68
5588 ISG20L2 1_77 0.001 6.48 1.5e-08 0.68
5590 HDGF 1_77 0.001 6.75 1.5e-08 0.61
10592 NTRK1 1_77 0.001 7.80 2.0e-08 -0.81
314 SH2D2A 1_77 0.001 6.20 1.3e-08 0.52
10039 PEAR1 1_77 0.001 11.90 4.8e-08 -1.26
4429 ARHGEF11 1_77 0.001 4.94 9.3e-09 0.21
12235 RP11-85G21.3 1_77 0.001 6.58 1.5e-08 0.64
5585 FCRL5 1_77 0.001 5.17 9.9e-09 -0.21
6928 FCRL3 1_77 0.001 6.32 1.4e-08 0.79
4432 FCRL2 1_77 0.002 13.06 6.1e-08 -1.52
7243 FCRL1 1_77 0.001 7.02 1.8e-08 -1.02
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 11_36"
genename region_tag susie_pip mu2 PVE z
11809 EEF1G 11_36 0.000 23.73 3.7e-15 -2.04
6091 EML3 11_36 0.000 30.93 1.0e-14 -2.40
9762 INTS5 11_36 0.000 23.89 3.8e-15 -2.05
6097 B3GAT3 11_36 0.000 21.41 2.6e-15 -1.91
1234 GANAB 11_36 0.000 73.57 2.1e-12 -3.89
11114 METTL12 11_36 0.000 33.31 1.4e-14 -2.50
12568 C11orf98 11_36 0.000 33.31 1.4e-14 2.50
10878 UQCC3 11_36 0.000 88.37 1.2e-11 4.29
7034 LBHD1 11_36 0.000 42.50 4.7e-14 -2.88
7033 GNG3 11_36 0.000 6.87 1.8e-16 0.68
8012 BSCL2 11_36 0.000 84.73 7.9e-12 4.19
7035 TTC9C 11_36 0.000 86.28 9.4e-12 -4.23
9809 TMEM179B 11_36 0.000 19.93 2.1e-15 1.83
8091 TMEM223 11_36 0.000 4.87 1.0e-16 -0.14
7037 STX5 11_36 0.000 15.30 9.8e-16 1.52
12607 RP11-727F15.14 11_36 0.000 15.25 9.7e-16 1.52
8013 HRASLS5 11_36 0.000 82.85 6.3e-12 4.14
4490 RARRES3 11_36 0.000 4.77 1.0e-16 0.03
9004 PLA2G16 11_36 0.000 33.71 1.5e-14 -2.52
9726 ATL3 11_36 0.000 11.97 5.4e-16 -1.26
4489 RTN3 11_36 0.000 14.15 8.0e-16 1.44
8014 C11orf84 11_36 0.000 7.59 2.2e-16 -0.79
803 MARK2 11_36 0.000 57.37 3.0e-13 3.40
2552 NAA40 11_36 0.000 114.95 2.5e-10 4.92
3901 PRDX5 11_36 0.000 6.77 1.8e-16 0.41
8986 COX8A 11_36 0.000 88.61 1.2e-11 -4.29
7973 OTUB1 11_36 0.000 88.68 1.2e-11 -4.29
3904 FLRT1 11_36 0.000 6.36 1.4e-16 2.37
4486 MACROD1 11_36 0.000 66.93 7.7e-13 5.16
6113 TRPT1 11_36 0.000 6.82 1.6e-16 1.40
8709 VEGFB 11_36 0.000 37.90 6.2e-15 2.93
2502 DNAJC4 11_36 0.000 14.28 5.2e-16 3.66
8708 FKBP2 11_36 0.000 81.64 3.3e-12 5.84
6114 PLCB3 11_36 0.000 38.19 3.1e-14 -4.91
15 BAD 11_36 0.000 6.89 1.7e-16 1.19
8666 ESRRA 11_36 0.000 24.31 4.4e-15 2.72
8660 TRMT112 11_36 0.000 32.42 5.8e-15 -3.51
8027 CCDC88B 11_36 0.000 22.67 7.5e-15 3.68
7041 RPS6KA4 11_36 0.000 26.66 2.2e-15 -5.96
11371 AP003774.6 11_36 0.000 5.92 1.3e-16 -2.02
9452 AP003774.4 11_36 0.000 11.56 3.8e-16 -1.23
8025 SF1 11_36 0.000 8.62 2.5e-16 1.32
4531 MEN1 11_36 0.000 65.18 9.2e-12 0.63
8432 CDC42BPG 11_36 0.000 125.35 1.1e-10 -8.15
2508 EHD1 11_36 0.000 116.13 2.5e-11 -6.82
2507 ATG2A 11_36 0.000 105.14 9.6e-12 -6.50
718 PPP2R5B 11_36 0.000 162.62 2.2e-09 8.17
8026 MAJIN 11_36 0.000 184.60 1.5e-08 -8.68
11026 ARL2 11_36 0.000 127.27 1.1e-09 4.34
8024 BATF2 11_36 0.000 9.13 2.0e-16 3.95
2504 SNX15 11_36 0.007 207.40 3.9e-06 -6.51
6116 TM7SF2 11_36 0.000 38.85 3.0e-15 -6.59
8786 ZNHIT2 11_36 0.000 38.54 2.9e-15 -6.57
6115 FAU 11_36 0.000 169.49 7.0e-08 6.98
10864 SPDYC 11_36 0.001 230.89 9.2e-07 -9.34
240 CAPN1 11_36 0.000 48.33 5.1e-14 4.29
238 POLA2 11_36 0.000 9.30 2.2e-16 3.75
4530 DPF2 11_36 0.000 10.75 5.4e-16 2.49
8745 TIGD3 11_36 0.000 17.13 1.5e-15 -1.57
7038 SLC25A45 11_36 0.000 10.39 5.2e-16 -1.80
11663 NEAT1 11_36 0.000 29.53 2.5e-14 -1.22
8022 LTBP3 11_36 0.000 18.39 3.1e-15 -1.32
9040 FAM89B 11_36 0.000 11.21 9.9e-16 -0.29
8693 KCNK7 11_36 0.000 114.27 7.4e-14 9.88
8691 MAP3K11 11_36 0.000 166.15 3.7e-15 15.38
10388 PCNX3 11_36 0.000 45.87 7.7e-14 5.06
11023 SIPA1 11_36 0.000 214.21 8.6e-11 12.60
8652 RELA 11_36 0.000 47.96 3.0e-15 5.28
8647 KAT5 11_36 0.000 65.35 2.1e-12 -4.87
8636 RNASEH2C 11_36 0.000 77.38 1.7e-12 6.70
11839 RP11-770G2.2 11_36 0.000 37.20 8.3e-16 6.82
11803 AP5B1 11_36 0.000 24.81 1.1e-14 0.42
8610 MUS81 11_36 0.000 63.38 3.0e-12 -0.55
8613 CFL1 11_36 0.000 57.18 9.0e-13 -0.96
8599 EFEMP2 11_36 0.000 7.74 1.6e-16 -1.54
8590 CTSW 11_36 0.000 22.38 1.4e-15 -2.31
8586 FIBP 11_36 0.000 85.14 8.9e-11 1.45
8924 CCDC85B 11_36 0.000 22.56 6.5e-16 4.98
8918 C11orf68 11_36 0.000 44.46 8.1e-15 5.40
8915 DRAP1 11_36 0.000 55.80 2.4e-15 8.69
8909 TSGA10IP 11_36 0.000 33.16 1.0e-15 7.26
8890 BANF1 11_36 0.000 60.39 1.3e-15 9.54
8888 CST6 11_36 0.000 24.19 4.8e-15 1.02
1166 SF3B2 11_36 0.000 16.65 4.7e-16 3.24
8863 PACS1 11_36 0.000 41.86 6.2e-15 -5.26
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 2_16"
genename region_tag susie_pip mu2 PVE z
11045 SLC35F6 2_16 0.015 18.26 7.8e-07 -3.70
3366 TMEM214 2_16 0.015 7.99 3.5e-07 1.10
5074 EMILIN1 2_16 0.021 29.91 1.8e-06 4.84
5061 KHK 2_16 0.013 6.76 2.6e-07 -0.99
5059 CGREF1 2_16 0.011 4.82 1.5e-07 -0.06
5070 PREB 2_16 0.011 8.91 3.0e-07 -2.00
5076 ATRAID 2_16 0.018 77.61 4.0e-06 -9.11
1090 CAD 2_16 0.046 21.48 2.9e-06 -3.07
5071 SLC5A6 2_16 0.018 16.85 8.9e-07 2.99
7303 UCN 2_16 0.016 20.68 9.6e-07 -4.41
2952 GTF3C2 2_16 0.017 21.62 1.1e-06 4.48
2956 SNX17 2_16 0.030 196.19 1.7e-05 14.39
7304 ZNF513 2_16 0.011 96.76 3.2e-06 10.03
2953 NRBP1 2_16 0.021 215.65 1.3e-05 15.35
5057 IFT172 2_16 0.027 32.83 2.5e-06 -5.41
1087 GCKR 2_16 0.029 34.30 2.9e-06 5.50
10613 GPN1 2_16 0.029 43.15 3.6e-06 6.44
9018 CCDC121 2_16 0.187 27.27 1.5e-05 1.38
6660 BRE 2_16 0.038 30.25 3.3e-06 4.67
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
#snps with PIP>0.8 or 20 highest PIPs
head(ctwas_snp_res[order(-ctwas_snp_res$susie_pip),report_cols_snps],
max(sum(ctwas_snp_res$susie_pip>0.8), 20))
id region_tag susie_pip mu2 PVE z
56720 rs766167074 1_118 1.000 5751.16 1.7e-02 -2.94
72399 rs780093 2_16 1.000 441.15 1.3e-03 -22.26
100027 rs882759 2_70 1.000 56.04 1.6e-04 -7.38
113237 rs7565788 2_103 1.000 140.49 4.1e-04 -11.91
113299 rs6433115 2_103 1.000 50.38 1.5e-04 7.40
115592 rs863678 2_106 1.000 102.08 3.0e-04 6.01
117533 rs11690832 2_110 1.000 59.14 1.7e-04 8.04
126605 rs2068218 2_128 1.000 44.83 1.3e-04 -6.60
198856 rs141435299 4_10 1.000 514.28 1.5e-03 -1.03
199060 rs57136958 4_11 1.000 377.71 1.1e-03 -10.08
199116 rs13115469 4_11 1.000 9484.50 2.8e-02 133.38
199151 rs3775948 4_11 1.000 9730.85 2.8e-02 131.05
199177 rs75968456 4_11 1.000 660.26 1.9e-03 -2.69
199480 rs6831973 4_12 1.000 85.49 2.5e-04 -9.83
199506 rs76285604 4_12 1.000 134.17 3.9e-04 -11.85
199593 rs142309009 4_12 1.000 73.37 2.1e-04 -9.06
225352 rs149027545 4_59 1.000 2233.20 6.5e-03 53.88
225418 rs10022462 4_60 1.000 99.49 2.9e-04 -12.99
229800 rs35518360 4_67 1.000 45.37 1.3e-04 -6.68
276354 rs255749 5_31 1.000 71.94 2.1e-04 8.64
277125 rs10077826 5_33 1.000 35.87 1.0e-04 -5.73
281562 rs10942549 5_43 1.000 170.65 5.0e-04 -15.21
318371 rs630258 6_7 1.000 163.02 4.7e-04 -17.15
324838 rs7754961 6_19 1.000 123.00 3.6e-04 15.01
324877 rs12213398 6_19 1.000 45.86 1.3e-04 -4.69
325521 rs6908155 6_21 1.000 63.06 1.8e-04 -2.88
329970 rs2856992 6_27 1.000 50.60 1.5e-04 -5.75
358655 rs10782229 6_84 1.000 42.78 1.2e-04 5.77
373878 rs78148157 7_3 1.000 114.56 3.3e-04 -9.55
373879 rs13241427 7_3 1.000 86.95 2.5e-04 8.10
400322 rs761767938 7_49 1.000 6271.03 1.8e-02 -3.88
400330 rs1544459 7_49 1.000 6306.74 1.8e-02 -4.00
449033 rs2941484 8_56 1.000 159.21 4.6e-04 15.49
483652 rs56030777 9_25 1.000 105.17 3.1e-04 -4.38
497215 rs1800977 9_53 1.000 48.11 1.4e-04 -6.96
524162 rs35182775 10_33 1.000 204.88 6.0e-04 15.09
527297 rs11510917 10_39 1.000 214.66 6.2e-04 -19.07
527310 rs1171619 10_39 1.000 332.94 9.7e-04 21.17
534222 rs149429992 10_52 1.000 7903.90 2.3e-02 2.53
545268 rs10886117 10_72 1.000 54.68 1.6e-04 7.40
562711 rs369062552 11_21 1.000 184.91 5.4e-04 11.72
562721 rs34830202 11_21 1.000 193.74 5.6e-04 -12.09
598886 rs11056397 12_13 1.000 35.29 1.0e-04 -5.84
611235 rs7397189 12_36 1.000 338.53 9.8e-04 -20.09
627821 rs653178 12_67 1.000 169.46 4.9e-04 -14.16
630757 rs2701194 12_74 1.000 71.82 2.1e-04 8.03
635090 rs76734539 12_82 1.000 60.40 1.8e-04 7.62
647865 rs566812111 13_25 1.000 2967.96 8.6e-03 2.56
677618 rs72681869 14_20 1.000 65.36 1.9e-04 -8.26
711310 rs2472297 15_35 1.000 66.83 1.9e-04 -9.72
711552 rs145727191 15_35 1.000 51.54 1.5e-04 8.91
711581 rs2955742 15_36 1.000 42.84 1.2e-04 7.22
719058 rs59646751 15_48 1.000 149.60 4.4e-04 14.47
732653 rs12927956 16_27 1.000 69.00 2.0e-04 7.55
760194 rs3794748 17_32 1.000 171.13 5.0e-04 -13.59
802776 rs56287732 19_23 1.000 40.84 1.2e-04 5.41
824742 rs209955 20_32 1.000 36.46 1.1e-04 5.89
836541 rs219783 21_16 1.000 88.59 2.6e-04 -9.60
868178 rs75246752 1_73 1.000 308.73 9.0e-04 -18.17
868325 rs185073199 1_73 1.000 281.42 8.2e-04 -18.84
930928 rs140927145 7_92 1.000 8892.68 2.6e-02 -2.70
940154 rs145700013 7_94 1.000 58.23 1.7e-04 -2.96
957040 rs542984928 11_36 1.000 241.27 7.0e-04 23.70
969052 rs6581124 12_35 1.000 109.18 3.2e-04 -10.40
12937 rs10788884 1_30 0.999 41.08 1.2e-04 6.38
507175 rs12380852 9_73 0.999 34.36 1.0e-04 5.71
552287 rs3842748 11_2 0.999 85.12 2.5e-04 -8.12
647869 rs12430288 13_25 0.999 2995.88 8.7e-03 2.63
755421 rs2688 17_22 0.999 32.30 9.4e-05 -5.49
957069 rs12363578 11_36 0.999 495.30 1.4e-03 -26.87
199276 rs4140694 4_11 0.998 864.91 2.5e-03 16.31
199599 rs61795273 4_12 0.998 54.59 1.6e-04 -7.25
534224 rs2152629 10_52 0.998 7922.61 2.3e-02 2.47
571707 rs1203614 11_37 0.998 50.15 1.5e-04 6.02
761515 rs11650989 17_36 0.998 39.50 1.1e-04 7.88
775289 rs527616 18_14 0.998 31.04 9.0e-05 -5.57
806242 rs814573 19_31 0.998 30.49 8.8e-05 -5.56
243285 rs4521364 4_95 0.997 56.37 1.6e-04 -5.84
314173 rs139078584 5_106 0.997 31.64 9.2e-05 6.53
617714 rs1848968 12_48 0.997 45.96 1.3e-04 -6.75
318194 rs200823080 6_6 0.996 33.74 9.8e-05 5.68
551072 rs2767419 10_85 0.996 34.60 1.0e-04 -5.65
179427 rs145422957 3_92 0.995 31.62 9.2e-05 -5.47
763623 rs113408695 17_39 0.995 31.74 9.2e-05 -5.55
204001 rs358256 4_20 0.994 33.89 9.8e-05 5.69
325354 rs13191326 6_21 0.994 231.67 6.7e-04 13.59
930924 rs6954405 7_92 0.994 8872.27 2.6e-02 -2.37
86985 rs7561263 2_46 0.993 35.22 1.0e-04 -6.03
173059 rs11712964 3_78 0.993 30.33 8.8e-05 5.45
713099 rs7174325 15_38 0.993 28.69 8.3e-05 4.89
826463 rs1407040 20_34 0.993 29.66 8.6e-05 -5.21
84547 rs778147 2_40 0.992 48.03 1.4e-04 6.79
512387 rs72777711 10_10 0.991 30.41 8.8e-05 -5.28
802630 rs71176165 19_23 0.991 69.48 2.0e-04 -9.09
719065 rs8037855 15_48 0.990 64.88 1.9e-04 10.82
92136 rs10196697 2_54 0.989 32.94 9.5e-05 -5.64
406723 rs45467892 7_61 0.989 41.19 1.2e-04 -6.54
715842 rs113404146 15_43 0.989 34.94 1.0e-04 5.81
400326 rs11972122 7_49 0.988 5827.76 1.7e-02 -3.92
692444 rs73349296 14_50 0.988 43.25 1.2e-04 -6.59
449071 rs2941432 8_56 0.987 69.94 2.0e-04 11.61
177014 rs115604285 3_87 0.986 32.60 9.3e-05 6.16
684650 rs10151620 14_34 0.985 36.60 1.0e-04 6.06
589995 rs73022311 11_77 0.984 26.57 7.6e-05 -4.89
738367 rs72799826 16_38 0.984 32.73 9.4e-05 -5.84
391262 rs700752 7_34 0.983 48.42 1.4e-04 6.67
844931 rs740219 22_10 0.982 35.74 1.0e-04 6.06
230575 rs2903386 4_69 0.979 34.27 9.8e-05 5.82
540386 rs7900763 10_64 0.978 26.87 7.6e-05 -4.72
394783 rs140420256 7_39 0.976 25.43 7.2e-05 4.74
609668 rs1878234 12_31 0.976 29.89 8.5e-05 -5.23
210003 rs12639940 4_32 0.974 25.93 7.3e-05 -4.41
330172 rs1126511 6_27 0.974 50.54 1.4e-04 5.39
630766 rs80019595 12_74 0.974 54.21 1.5e-04 -6.72
802814 rs889140 19_23 0.971 29.70 8.4e-05 -4.20
572155 rs6591334 11_37 0.970 38.19 1.1e-04 -5.74
712955 rs8024096 15_37 0.970 28.38 8.0e-05 5.11
992645 rs3760994 19_2 0.970 46.54 1.3e-04 -6.66
150789 rs62259692 3_36 0.969 35.30 9.9e-05 7.65
204798 rs34811474 4_21 0.969 35.24 9.9e-05 -5.82
429068 rs17151140 8_13 0.969 55.65 1.6e-04 5.11
738028 rs56259873 16_37 0.969 90.29 2.5e-04 10.62
775330 rs162000 18_14 0.968 25.57 7.2e-05 4.95
237565 rs62323480 4_83 0.966 25.77 7.2e-05 4.40
247282 rs139212650 4_102 0.965 24.93 7.0e-05 4.61
6284 rs10917555 1_13 0.964 25.92 7.3e-05 -4.75
34463 rs1979575 1_75 0.964 25.34 7.1e-05 -4.39
422135 rs11761498 7_98 0.961 32.78 9.2e-05 5.26
479025 rs10964603 9_16 0.961 25.17 7.0e-05 -4.52
277286 rs536916238 5_33 0.959 31.98 8.9e-05 -0.30
922914 rs56142449 6_61 0.958 47.74 1.3e-04 7.05
438509 rs111965375 8_34 0.955 25.87 7.2e-05 5.22
570268 rs199725747 11_32 0.955 31.05 8.6e-05 5.34
225372 rs28366540 4_59 0.952 558.20 1.5e-03 -33.87
333088 rs10223666 6_34 0.952 191.43 5.3e-04 14.34
70389 rs62112223 2_12 0.951 30.04 8.3e-05 -5.32
826332 rs73185042 20_34 0.950 25.34 7.0e-05 -4.69
316516 rs1272694 6_3 0.949 34.82 9.6e-05 -5.95
798130 rs11668601 19_14 0.948 34.35 9.5e-05 7.11
326301 rs13437375 6_23 0.944 29.96 8.2e-05 2.72
49268 rs114934802 1_104 0.941 23.80 6.5e-05 -3.09
277268 rs173964 5_33 0.941 104.26 2.9e-04 8.09
463617 rs921719 8_83 0.940 27.15 7.4e-05 -5.17
868345 rs139428292 1_73 0.940 75.54 2.1e-04 -9.43
140620 rs6778028 3_12 0.937 25.47 6.9e-05 4.59
424506 rs117950418 8_4 0.937 24.64 6.7e-05 -4.62
683084 rs34528648 14_33 0.935 32.20 8.8e-05 5.45
721644 rs2601781 16_4 0.933 25.40 6.9e-05 4.71
719041 rs75422555 15_47 0.932 31.13 8.4e-05 -6.41
439799 rs12543287 8_37 0.931 34.57 9.4e-05 -5.67
157749 rs56324130 3_49 0.927 23.89 6.4e-05 -4.43
218953 rs144193499 4_48 0.927 25.06 6.8e-05 4.55
505329 rs201421930 9_69 0.924 34.47 9.3e-05 -5.88
551018 rs75184896 10_84 0.922 26.72 7.2e-05 5.47
193668 rs7642977 3_118 0.921 29.16 7.8e-05 5.14
352883 rs12196331 6_71 0.921 28.67 7.7e-05 5.24
718742 rs116887089 15_47 0.920 24.30 6.5e-05 4.45
236735 rs72680231 4_81 0.919 24.59 6.6e-05 -4.66
185840 rs1290790 3_104 0.916 69.97 1.9e-04 6.15
610857 rs11608918 12_33 0.916 26.13 7.0e-05 4.71
235872 rs41278087 4_79 0.915 24.02 6.4e-05 -4.60
464936 rs36041912 8_85 0.911 24.43 6.5e-05 -4.56
522457 rs58434594 10_30 0.911 27.86 7.4e-05 4.83
154979 rs56145049 3_43 0.910 24.81 6.6e-05 -4.67
271612 rs13170671 5_23 0.909 70.56 1.9e-04 8.54
315427 rs6873880 5_109 0.904 27.94 7.3e-05 -5.03
346960 rs118165878 6_58 0.901 23.85 6.2e-05 4.44
211738 rs112161979 4_35 0.899 23.68 6.2e-05 4.56
229866 rs13140033 4_68 0.898 24.39 6.4e-05 -4.46
66812 rs139638572 2_6 0.896 28.26 7.4e-05 4.44
150972 rs2276816 3_36 0.892 32.58 8.5e-05 -0.96
72596 rs7606480 2_19 0.891 26.38 6.8e-05 -4.89
737597 rs12599580 16_36 0.891 27.48 7.1e-05 5.66
783917 rs784218 18_30 0.891 27.10 7.0e-05 4.29
711522 rs553274247 15_35 0.890 122.43 3.2e-04 -5.27
280971 rs6886422 5_42 0.886 25.31 6.5e-05 -4.52
447667 rs34650318 8_53 0.884 24.60 6.3e-05 -4.51
189108 rs17593458 3_110 0.880 25.58 6.5e-05 4.52
287010 rs3952745 5_53 0.874 32.71 8.3e-05 -6.44
8429 rs7516039 1_20 0.873 25.70 6.5e-05 4.71
324791 rs62392365 6_19 0.873 72.25 1.8e-04 11.55
910347 rs78595810 3_33 0.873 37.99 9.6e-05 6.24
354984 rs1997649 6_76 0.870 23.80 6.0e-05 4.28
572496 rs10752584 11_38 0.867 30.20 7.6e-05 4.73
613671 rs113897089 12_40 0.866 29.67 7.5e-05 5.22
544001 rs11196217 10_70 0.864 25.06 6.3e-05 -4.56
198439 rs28680668 4_9 0.863 33.18 8.3e-05 5.66
810362 rs2003700 20_1 0.863 24.39 6.1e-05 4.41
552274 rs11042594 11_2 0.862 67.85 1.7e-04 6.63
956689 rs78915580 11_36 0.862 65.08 1.6e-04 -7.09
792984 rs2074457 19_4 0.861 25.90 6.5e-05 4.44
5089 rs4336844 1_11 0.857 34.05 8.5e-05 4.45
571730 rs509533 11_37 0.855 45.23 1.1e-04 7.93
959380 rs5792371 11_36 0.853 267.54 6.6e-04 18.90
413023 rs75520353 7_74 0.852 24.12 6.0e-05 4.36
227402 rs11722010 4_63 0.851 25.08 6.2e-05 -4.49
698690 rs62004133 15_8 0.851 25.21 6.2e-05 -4.56
741707 rs76862947 16_44 0.849 139.45 3.4e-04 -11.76
26425 rs9432440 1_58 0.845 25.44 6.3e-05 4.55
483621 rs10971408 9_25 0.845 28.05 6.9e-05 -3.63
99801 rs2311597 2_70 0.838 59.90 1.5e-04 -9.93
40871 rs604388 1_87 0.836 24.72 6.0e-05 4.43
446348 rs71553284 8_50 0.836 25.43 6.2e-05 4.39
785487 rs17773471 18_33 0.835 35.77 8.7e-05 7.37
103829 rs71862935 2_79 0.832 27.40 6.6e-05 -4.89
737963 rs72799382 16_37 0.832 42.87 1.0e-04 -6.78
110874 rs12467636 2_96 0.830 35.47 8.6e-05 5.79
122243 rs56059523 2_120 0.829 29.74 7.2e-05 4.32
302294 rs356486 5_82 0.828 29.34 7.1e-05 5.01
150575 rs78342753 3_35 0.824 34.63 8.3e-05 5.19
407387 rs7803747 7_63 0.824 26.62 6.4e-05 5.00
601668 rs4077798 12_18 0.824 30.58 7.3e-05 -5.30
399594 rs1179610 7_48 0.820 26.15 6.2e-05 4.72
805138 rs11671669 19_30 0.816 26.29 6.2e-05 4.49
580024 rs57569860 11_52 0.815 23.94 5.7e-05 3.53
139938 rs9846767 3_11 0.811 25.13 5.9e-05 4.59
87000 rs72837990 2_46 0.810 24.87 5.9e-05 -4.42
435464 rs4871905 8_24 0.810 131.51 3.1e-04 12.10
824772 rs570322378 20_32 0.809 27.77 6.5e-05 4.45
57987 rs678615 1_122 0.807 40.29 9.5e-05 -6.17
640192 rs9315009 13_9 0.807 36.78 8.6e-05 5.93
547754 rs10794177 10_78 0.806 52.00 1.2e-04 7.11
547877 rs1693628 10_78 0.802 26.48 6.2e-05 -4.38
#plot PIP vs effect size
plot(ctwas_snp_res$susie_pip, ctwas_snp_res$mu2, xlab="PIP", ylab="mu^2", main="SNP PIPs vs Effect Size")
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
#SNPs with 50 largest effect sizes
head(ctwas_snp_res[order(-ctwas_snp_res$mu2),report_cols_snps],50)
id region_tag susie_pip mu2 PVE z
199151 rs3775948 4_11 1.000 9730.85 2.8e-02 131.05
199118 rs6838021 4_11 0.000 9648.29 0.0e+00 133.98
199119 rs6823324 4_11 0.000 9629.18 0.0e+00 -133.15
199122 rs11723439 4_11 0.000 9489.36 0.0e+00 -132.41
199116 rs13115469 4_11 1.000 9484.50 2.8e-02 133.38
930928 rs140927145 7_92 1.000 8892.68 2.6e-02 -2.70
930924 rs6954405 7_92 0.994 8872.27 2.6e-02 -2.37
930923 rs6953180 7_92 0.299 8870.22 7.7e-03 -2.31
930922 rs11971515 7_92 0.090 8868.53 2.3e-03 -2.27
930921 rs11974627 7_92 0.101 8867.27 2.6e-03 -2.27
930926 rs528981137 7_92 0.006 8864.62 1.6e-04 -2.16
930927 rs549027364 7_92 0.007 8864.41 1.7e-04 -2.16
199141 rs7439210 4_11 0.000 8719.19 0.0e+00 -128.12
930918 rs7807788 7_92 0.000 8649.45 8.7e-08 -2.43
930912 rs11546289 7_92 0.000 8625.48 1.9e-09 -2.37
930913 rs13243678 7_92 0.000 8616.95 7.4e-10 -2.36
930914 rs11546290 7_92 0.000 8611.36 2.4e-10 -2.34
930910 rs6952916 7_92 0.000 8570.29 4.4e-12 -2.29
930907 rs7793916 7_92 0.000 8558.99 3.0e-12 -2.30
930908 rs7781312 7_92 0.000 8539.69 4.1e-13 -2.28
930906 rs13247593 7_92 0.000 8530.79 8.3e-12 -2.38
930909 rs6965276 7_92 0.000 8484.08 6.6e-15 -2.28
930931 rs112579924 7_92 0.000 8444.92 0.0e+00 -1.71
930905 rs568222600 7_92 0.000 8406.53 3.5e-17 -2.28
930936 rs12667507 7_92 0.000 8396.38 0.0e+00 2.10
930903 rs112799047 7_92 0.000 8281.08 0.0e+00 -2.31
930904 rs10233906 7_92 0.000 8270.55 0.0e+00 -2.31
930901 rs112684174 7_92 0.000 8261.81 0.0e+00 -2.32
930902 rs112634954 7_92 0.000 8253.40 0.0e+00 -2.32
930896 rs6965440 7_92 0.000 8240.31 0.0e+00 -2.33
930895 rs6464930 7_92 0.000 8228.70 0.0e+00 -2.35
930890 rs6978068 7_92 0.000 8187.98 0.0e+00 -2.28
930881 rs10233535 7_92 0.000 8087.97 0.0e+00 -2.31
930882 rs56932055 7_92 0.000 8085.60 0.0e+00 -2.32
930886 rs6958855 7_92 0.000 8084.60 0.0e+00 -2.34
930889 rs6963547 7_92 0.000 8081.91 0.0e+00 -2.37
930871 rs10269104 7_92 0.000 8080.73 0.0e+00 -2.30
930887 rs1551926 7_92 0.000 8058.67 0.0e+00 2.45
199209 rs11723742 4_11 0.000 7995.98 0.0e+00 -116.64
930915 rs2306169 7_92 0.000 7930.34 0.0e+00 2.54
534224 rs2152629 10_52 0.998 7922.61 2.3e-02 2.47
930917 rs2306170 7_92 0.000 7913.74 0.0e+00 2.33
534220 rs7913261 10_52 0.562 7905.82 1.3e-02 2.51
534222 rs149429992 10_52 1.000 7903.90 2.3e-02 2.53
199293 rs17389602 4_11 0.000 7847.08 0.0e+00 -117.34
199295 rs78917351 4_11 0.000 7835.51 0.0e+00 -117.31
534225 rs1360953 10_52 0.000 7825.99 1.6e-08 2.50
534230 rs76574695 10_52 0.000 7824.36 1.4e-07 2.40
930925 rs6953222 7_92 0.000 7824.06 0.0e+00 -1.43
930920 rs34909003 7_92 0.000 7772.17 0.0e+00 -1.86
#SNPs with 50 highest pve
head(ctwas_snp_res[order(-ctwas_snp_res$PVE),report_cols_snps],50)
id region_tag susie_pip mu2 PVE z
199116 rs13115469 4_11 1.000 9484.50 0.02800 133.38
199151 rs3775948 4_11 1.000 9730.85 0.02800 131.05
930924 rs6954405 7_92 0.994 8872.27 0.02600 -2.37
930928 rs140927145 7_92 1.000 8892.68 0.02600 -2.70
534222 rs149429992 10_52 1.000 7903.90 0.02300 2.53
534224 rs2152629 10_52 0.998 7922.61 0.02300 2.47
400322 rs761767938 7_49 1.000 6271.03 0.01800 -3.88
400330 rs1544459 7_49 1.000 6306.74 0.01800 -4.00
56720 rs766167074 1_118 1.000 5751.16 0.01700 -2.94
400326 rs11972122 7_49 0.988 5827.76 0.01700 -3.92
534220 rs7913261 10_52 0.562 7905.82 0.01300 2.51
647869 rs12430288 13_25 0.999 2995.88 0.00870 2.63
647865 rs566812111 13_25 1.000 2967.96 0.00860 2.56
930923 rs6953180 7_92 0.299 8870.22 0.00770 -2.31
225352 rs149027545 4_59 1.000 2233.20 0.00650 53.88
56718 rs2486737 1_118 0.370 5713.07 0.00620 -3.17
56719 rs971534 1_118 0.358 5713.05 0.00590 -3.16
56717 rs10489611 1_118 0.304 5712.90 0.00500 -3.16
56714 rs2790891 1_118 0.290 5712.52 0.00480 -3.16
56715 rs2491405 1_118 0.290 5712.52 0.00480 -3.16
56711 rs2256908 1_118 0.277 5712.49 0.00460 -3.16
647868 rs1579715 13_25 0.388 2965.08 0.00330 -2.77
56726 rs2248646 1_118 0.181 5710.15 0.00300 -3.14
56727 rs2211176 1_118 0.173 5710.35 0.00290 -3.14
56728 rs2790882 1_118 0.173 5710.35 0.00290 -3.14
930921 rs11974627 7_92 0.101 8867.27 0.00260 -2.27
199276 rs4140694 4_11 0.998 864.91 0.00250 16.31
930922 rs11971515 7_92 0.090 8868.53 0.00230 -2.27
199177 rs75968456 4_11 1.000 660.26 0.00190 -2.69
56707 rs1076804 1_118 0.089 5703.74 0.00150 -3.14
198856 rs141435299 4_10 1.000 514.28 0.00150 -1.03
225372 rs28366540 4_59 0.952 558.20 0.00150 -33.87
957069 rs12363578 11_36 0.999 495.30 0.00140 -26.87
72399 rs780093 2_16 1.000 441.15 0.00130 -22.26
56729 rs1416913 1_118 0.071 5703.06 0.00120 -3.13
56723 rs2739509 1_118 0.070 5609.25 0.00110 -3.30
199060 rs57136958 4_11 1.000 377.71 0.00110 -10.08
198851 rs146530806 4_10 0.572 606.87 0.00100 6.59
198849 rs144362537 4_10 0.558 607.49 0.00099 6.57
611235 rs7397189 12_36 1.000 338.53 0.00098 -20.09
527310 rs1171619 10_39 1.000 332.94 0.00097 21.17
198850 rs34658640 4_10 0.536 607.42 0.00095 6.57
868178 rs75246752 1_73 1.000 308.73 0.00090 -18.17
56732 rs2790874 1_118 0.051 5702.25 0.00085 -3.12
868325 rs185073199 1_73 1.000 281.42 0.00082 -18.84
957040 rs542984928 11_36 1.000 241.27 0.00070 23.70
325354 rs13191326 6_21 0.994 231.67 0.00067 13.59
959380 rs5792371 11_36 0.853 267.54 0.00066 18.90
527297 rs11510917 10_39 1.000 214.66 0.00062 -19.07
524162 rs35182775 10_33 1.000 204.88 0.00060 15.09
#SNPs with 50 largest z scores
head(ctwas_snp_res[order(-abs(ctwas_snp_res$z)),report_cols_snps],50)
id region_tag susie_pip mu2 PVE z
199118 rs6838021 4_11 0 9648.29 0.000 133.98
199116 rs13115469 4_11 1 9484.50 0.028 133.38
199119 rs6823324 4_11 0 9629.18 0.000 -133.15
199122 rs11723439 4_11 0 9489.36 0.000 -132.41
199151 rs3775948 4_11 1 9730.85 0.028 131.05
199141 rs7439210 4_11 0 8719.19 0.000 -128.12
199293 rs17389602 4_11 0 7847.08 0.000 -117.34
199295 rs78917351 4_11 0 7835.51 0.000 -117.31
199209 rs11723742 4_11 0 7995.98 0.000 -116.64
199347 rs11722185 4_11 0 7644.92 0.000 -115.64
199333 rs11727390 4_11 0 7650.25 0.000 -115.57
199339 rs4697745 4_11 0 7032.95 0.000 -110.85
199370 rs546391476 4_11 0 6948.16 0.000 -110.51
199321 rs10489070 4_11 0 6827.95 0.000 -109.37
199159 rs9291642 4_11 0 5586.71 0.000 103.74
199256 rs4697717 4_11 0 6463.32 0.000 -103.72
199278 rs887732 4_11 0 6247.19 0.000 -103.53
199181 rs7349721 4_11 0 5360.29 0.000 100.68
199125 rs7375643 4_11 0 4136.86 0.000 -90.99
199145 rs11723970 4_11 0 4130.56 0.000 90.75
199124 rs7375599 4_11 0 4108.12 0.000 -90.71
199377 rs5856057 4_11 0 4544.79 0.000 -86.60
199139 rs6449177 4_11 0 3551.11 0.000 -85.49
199152 rs34501273 4_11 0 3886.64 0.000 84.99
199154 rs3733586 4_11 0 3887.21 0.000 84.99
199155 rs35438220 4_11 0 3886.06 0.000 84.99
199158 rs12507330 4_11 0 3885.14 0.000 84.97
199149 rs17187075 4_11 0 3737.74 0.000 84.09
199121 rs12498956 4_11 0 3434.53 0.000 -83.90
199163 rs3756236 4_11 0 3719.29 0.000 83.79
199178 rs11727199 4_11 0 3719.96 0.000 83.20
199179 rs10939665 4_11 0 3717.49 0.000 83.17
199180 rs12508991 4_11 0 3718.21 0.000 83.17
199192 rs35501905 4_11 0 3715.15 0.000 83.12
199113 rs35955619 4_11 0 3189.77 0.000 -82.89
199131 rs6844316 4_11 0 3332.72 0.000 -82.84
199136 rs4295261 4_11 0 3334.21 0.000 -82.84
199129 rs34297373 4_11 0 3331.51 0.000 -82.82
199132 rs28837683 4_11 0 3327.14 0.000 -82.81
199138 rs7434391 4_11 0 3334.52 0.000 -82.79
199175 rs13148371 4_11 0 3696.12 0.000 82.73
199123 rs7376155 4_11 0 3308.95 0.000 -82.67
199127 rs4314284 4_11 0 3308.99 0.000 -82.67
199196 rs7678211 4_11 0 3824.66 0.000 82.60
199217 rs4235356 4_11 0 3341.83 0.000 -80.91
199166 rs6827785 4_11 0 3247.91 0.000 79.26
199191 rs13120348 4_11 0 3236.71 0.000 78.51
199183 rs13122689 4_11 0 3227.14 0.000 78.50
199184 rs12504565 4_11 0 3227.57 0.000 78.50
199176 rs12506122 4_11 0 3240.61 0.000 78.48
#GO enrichment analysis
library(enrichR)
Welcome to enrichR
Checking connection ...
Enrichr ... Connection is Live!
FlyEnrichr ... Connection is available!
WormEnrichr ... Connection is available!
YeastEnrichr ... Connection is available!
FishEnrichr ... Connection is available!
dbs <- c("GO_Biological_Process_2021", "GO_Cellular_Component_2021", "GO_Molecular_Function_2021")
genes <- ctwas_gene_res$genename[ctwas_gene_res$susie_pip>0.8]
#number of genes for gene set enrichment
length(genes)
[1] 15
if (length(genes)>0){
GO_enrichment <- enrichr(genes, dbs)
for (db in dbs){
print(db)
df <- GO_enrichment[[db]]
df <- df[df$Adjusted.P.value<0.05,c("Term", "Overlap", "Adjusted.P.value", "Genes")]
print(df)
}
#DisGeNET enrichment
# devtools::install_bitbucket("ibi_group/disgenet2r")
library(disgenet2r)
disgenet_api_key <- get_disgenet_api_key(
email = "wesleycrouse@gmail.com",
password = "uchicago1" )
Sys.setenv(DISGENET_API_KEY= disgenet_api_key)
res_enrich <-disease_enrichment(entities=genes, vocabulary = "HGNC",
database = "CURATED" )
df <- res_enrich@qresult[1:10, c("Description", "FDR", "Ratio", "BgRatio")]
print(df)
#WebGestalt enrichment
library(WebGestaltR)
background <- ctwas_gene_res$genename
#listGeneSet()
databases <- c("pathway_KEGG", "disease_GLAD4U", "disease_OMIM")
enrichResult <- WebGestaltR(enrichMethod="ORA", organism="hsapiens",
interestGene=genes, referenceGene=background,
enrichDatabase=databases, interestGeneType="genesymbol",
referenceGeneType="genesymbol", isOutput=F)
print(enrichResult[,c("description", "size", "overlap", "FDR", "database", "userId")])
}
Uploading data to Enrichr... Done.
Querying GO_Biological_Process_2021... Done.
Querying GO_Cellular_Component_2021... Done.
Querying GO_Molecular_Function_2021... Done.
Parsing results... Done.
[1] "GO_Biological_Process_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)
[1] "GO_Cellular_Component_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)
[1] "GO_Molecular_Function_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)
ZNF276 gene(s) from the input list not found in DisGeNET CURATEDZDHHC18 gene(s) from the input list not found in DisGeNET CURATEDANO7 gene(s) from the input list not found in DisGeNET CURATEDTLCD2 gene(s) from the input list not found in DisGeNET CURATEDKCNK17 gene(s) from the input list not found in DisGeNET CURATEDAGAP3 gene(s) from the input list not found in DisGeNET CURATEDCSNK1G2 gene(s) from the input list not found in DisGeNET CURATED
Description FDR Ratio BgRatio
61 Medullary cystic kidney disease 1 0.02061431 1/8 1/9703
64 Mental Retardation, Autosomal Dominant 1 0.02061431 1/8 1/9703
72 2q23.1 microdeletion syndrome 0.02061431 1/8 1/9703
4 Cannabis Dependence 0.07375783 1/8 17/9703
5 Neoplastic Cell Transformation 0.07375783 2/8 139/9703
22 Neoplasm Invasiveness 0.07375783 2/8 184/9703
25 Peritoneal Neoplasms 0.07375783 1/8 10/9703
31 Gastric ulcer 0.07375783 1/8 18/9703
32 Aganglionosis, Colonic 0.07375783 1/8 15/9703
48 Carcinomatosis of peritoneal cavity 0.07375783 1/8 10/9703
******************************************
* *
* Welcome to WebGestaltR ! *
* *
******************************************
Loading the functional categories...
Loading the ID list...
Loading the reference list...
Performing the enrichment analysis...
Warning in oraEnrichment(interestGeneList, referenceGeneList, geneSet,
minNum = minNum, : No significant gene set is identified based on FDR 0.05!
NULL
sessionInfo()
R version 3.6.1 (2019-07-05)
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] WebGestaltR_0.4.4 disgenet2r_0.99.2 enrichR_3.0 cowplot_1.0.0
[5] ggplot2_3.3.3
loaded via a namespace (and not attached):
[1] bitops_1.0-6 matrixStats_0.57.0
[3] fs_1.3.1 bit64_4.0.5
[5] doParallel_1.0.16 progress_1.2.2
[7] httr_1.4.1 rprojroot_2.0.2
[9] GenomeInfoDb_1.20.0 doRNG_1.8.2
[11] tools_3.6.1 utf8_1.2.1
[13] R6_2.5.0 DBI_1.1.1
[15] BiocGenerics_0.30.0 colorspace_1.4-1
[17] withr_2.4.1 tidyselect_1.1.0
[19] prettyunits_1.0.2 bit_4.0.4
[21] curl_3.3 compiler_3.6.1
[23] git2r_0.26.1 Biobase_2.44.0
[25] DelayedArray_0.10.0 rtracklayer_1.44.0
[27] labeling_0.3 scales_1.1.0
[29] readr_1.4.0 apcluster_1.4.8
[31] stringr_1.4.0 digest_0.6.20
[33] Rsamtools_2.0.0 svglite_1.2.2
[35] rmarkdown_1.13 XVector_0.24.0
[37] pkgconfig_2.0.3 htmltools_0.3.6
[39] fastmap_1.1.0 BSgenome_1.52.0
[41] rlang_0.4.11 RSQLite_2.2.7
[43] generics_0.0.2 farver_2.1.0
[45] jsonlite_1.6 BiocParallel_1.18.0
[47] dplyr_1.0.7 VariantAnnotation_1.30.1
[49] RCurl_1.98-1.1 magrittr_2.0.1
[51] GenomeInfoDbData_1.2.1 Matrix_1.2-18
[53] Rcpp_1.0.6 munsell_0.5.0
[55] S4Vectors_0.22.1 fansi_0.5.0
[57] gdtools_0.1.9 lifecycle_1.0.0
[59] stringi_1.4.3 whisker_0.3-2
[61] yaml_2.2.0 SummarizedExperiment_1.14.1
[63] zlibbioc_1.30.0 plyr_1.8.4
[65] grid_3.6.1 blob_1.2.1
[67] parallel_3.6.1 promises_1.0.1
[69] crayon_1.4.1 lattice_0.20-38
[71] Biostrings_2.52.0 GenomicFeatures_1.36.3
[73] hms_1.1.0 knitr_1.23
[75] pillar_1.6.1 igraph_1.2.4.1
[77] GenomicRanges_1.36.0 rjson_0.2.20
[79] rngtools_1.5 codetools_0.2-16
[81] reshape2_1.4.3 biomaRt_2.40.1
[83] stats4_3.6.1 XML_3.98-1.20
[85] glue_1.4.2 evaluate_0.14
[87] data.table_1.14.0 foreach_1.5.1
[89] vctrs_0.3.8 httpuv_1.5.1
[91] gtable_0.3.0 purrr_0.3.4
[93] assertthat_0.2.1 cachem_1.0.5
[95] xfun_0.8 later_0.8.0
[97] tibble_3.1.2 iterators_1.0.13
[99] GenomicAlignments_1.20.1 AnnotationDbi_1.46.0
[101] memoise_2.0.0 IRanges_2.18.1
[103] workflowr_1.6.2 ellipsis_0.3.2