Albumin (quantile) - Liver
wesleycrouse
2021-09-09
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-30600_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-30600_irnt_Liver.Rmd) and HTML (docs/ukb-d-30600_irnt_Liver.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 | 627a4e1 | wesleycrouse | 2021-09-07 | adding heritability |
| 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 | 05a98b7 | wesleycrouse | 2021-08-07 | adding additional results |
| html | 05a98b7 | wesleycrouse | 2021-08-07 | adding additional results |
| html | 03e541c | wesleycrouse | 2021-07-29 | Cleaning up report generation |
| Rmd | 276893d | wesleycrouse | 2021-07-29 | Updating reports |
| html | 276893d | wesleycrouse | 2021-07-29 | Updating reports |
Overview
These are the results of a ctwas analysis of the UK Biobank trait Albumin (quantile) using Liver 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-30600_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 Liver 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])
Weight QC
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] 10901#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
1070 768 652 417 494 611 548 408 405 434 634 629 195 365 354
16 17 18 19 20 21 22
526 663 160 859 306 114 289 #proportion of imputed weights without missing variants
mean(qclist_all$nmiss==0)[1] 0.8366205Load ctwas results
#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)Check convergence of parameters
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.0201865738 0.0001964596 #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
29.28187 12.24032 #report sample size
print(sample_size)[1] 315268#report group size
group_size <- c(nrow(ctwas_gene_res), n_snps)
print(group_size)[1] 10901 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.02043845 0.06633948 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.5564539 1.5878886Genes with highest PIPs
#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
1144 ASAP3 1_16 1.000 55.42 1.8e-04 7.28
5789 ABRACL 6_92 1.000 180.93 5.7e-04 6.64
2173 TMEM176B 7_93 1.000 218.91 6.9e-04 -14.87
12008 HPR 16_38 1.000 130.98 4.2e-04 -9.93
5389 RPS11 19_34 1.000 163204.55 5.2e-01 17.58
3212 CCND2 12_4 0.996 64.82 2.0e-04 -7.54
11885 SLC35G6 17_7 0.996 79.32 2.5e-04 -9.53
1429 SH3BP1 22_15 0.995 53.59 1.7e-04 8.26
6936 RAVER2 1_41 0.994 48.81 1.5e-04 6.79
257 RUNX3 1_17 0.993 36.98 1.2e-04 -5.71
7656 CATSPER2 15_16 0.993 303.50 9.6e-04 18.82
12467 RP11-219B17.3 15_27 0.993 85.74 2.7e-04 -8.67
6121 ZNF827 4_95 0.992 30.70 9.7e-05 -5.13
6220 PELO 5_31 0.991 37.66 1.2e-04 -5.94
6494 PHKG2 16_24 0.988 38.39 1.2e-04 5.58
1231 PABPC4 1_24 0.987 79.71 2.5e-04 9.17
9613 C14orf80 14_55 0.985 46.21 1.4e-04 10.26
12074 RP11-131K5.2 17_12 0.983 40.49 1.3e-04 -6.20
2809 COL7A1 3_35 0.982 81.43 2.5e-04 9.04
7084 EIF4E3 3_48 0.973 41.36 1.3e-04 6.47
8746 LRRC25 19_15 0.972 46.49 1.4e-04 9.38
2341 DDX5 17_37 0.971 23.50 7.2e-05 4.51
8531 TNKS 8_12 0.963 43.46 1.3e-04 8.95
10602 RNF5 6_26 0.958 46.40 1.4e-04 7.83
8107 PCDH7 4_26 0.957 26.70 8.1e-05 4.97
10763 NYNRIN 14_3 0.947 52.70 1.6e-04 7.15
5799 SLC22A3 6_104 0.945 162.77 4.9e-04 -6.89
6093 CSNK1G3 5_75 0.942 78.18 2.3e-04 9.49
11237 TNF 6_25 0.937 56.02 1.7e-04 5.24
11478 HLA-DMB 6_27 0.934 38.72 1.1e-04 -6.75
3591 SDC4 20_28 0.933 31.83 9.4e-05 -5.48
2142 EIF3B 7_4 0.932 35.00 1.0e-04 -6.46
7915 GLYCTK 3_36 0.927 32.92 9.7e-05 -4.96
11389 GATSL3 22_10 0.923 28.18 8.2e-05 -5.03
4638 GLUL 1_90 0.914 23.29 6.8e-05 -4.29
708 PTPN3 9_55 0.907 34.65 1.0e-04 5.74
8614 FAM53A 4_3 0.906 24.43 7.0e-05 4.57
906 UBE2K 4_32 0.903 30.42 8.7e-05 -5.23
11558 LINC01184 5_78 0.903 57.96 1.7e-04 7.46
7682 VPS39 15_15 0.894 23.52 6.7e-05 -5.06
2147 MOSPD3 7_62 0.893 25.77 7.3e-05 -4.73
2474 CBL 11_71 0.892 42.81 1.2e-04 -6.57
9855 PALM3 19_11 0.886 20.25 5.7e-05 3.92
1455 JOSD1 22_15 0.885 37.81 1.1e-04 6.20
7098 EOMES 3_20 0.878 27.57 7.7e-05 -4.85
6637 NPM2 8_23 0.853 33.70 9.1e-05 4.54
4842 CPEB2 4_14 0.851 44.12 1.2e-04 6.37
6702 ACE 17_37 0.844 30.58 8.2e-05 -5.09
5089 SCAF11 12_29 0.840 22.03 5.9e-05 -4.17
11306 LINC00094 9_70 0.837 34.20 9.1e-05 5.52
6919 TAL1 1_29 0.831 33.94 8.9e-05 -4.79
2645 ASCC3 6_68 0.831 22.00 5.8e-05 4.39
10880 MEF2B 19_15 0.831 25.59 6.7e-05 -5.61
2998 RALGPS2 1_87 0.817 90.60 2.3e-04 9.36
10442 GLMP 1_76 0.805 35.20 9.0e-05 6.67Genes with largest effect sizes
#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
5389 RPS11 19_34 1 163204.55 5.2e-01 17.58
1227 FLT3LG 19_34 0 141042.21 0.0e+00 -16.55
5393 RCN3 19_34 0 52574.94 0.0e+00 -14.77
1931 FCGRT 19_34 0 47948.16 0.0e+00 -6.60
3804 PRRG2 19_34 0 23602.17 0.0e+00 -20.72
3803 PRMT1 19_34 0 16030.76 0.0e+00 -9.25
3805 SCAF1 19_34 0 15951.61 0.0e+00 -11.38
3802 IRF3 19_34 0 15477.49 0.0e+00 -10.77
1940 SLC17A7 19_34 0 11229.71 0.0e+00 -4.42
4556 TMEM60 7_49 0 7106.13 0.0e+00 5.13
1932 PIH1D1 19_34 0 4923.97 0.0e+00 3.78
10291 CTC-301O7.4 19_34 0 3848.50 0.0e+00 2.57
4634 EGLN1 1_118 0 2875.12 1.3e-08 2.92
3058 EXOC8 1_118 0 2433.35 2.9e-12 -3.91
8030 CPT1C 19_34 0 2092.69 0.0e+00 0.75
545 SLC6A16 19_34 0 1462.87 0.0e+00 0.40
10903 APTR 7_49 0 1353.58 0.0e+00 1.63
5390 RPL13A 19_34 0 1127.45 0.0e+00 -7.77
2486 PTPMT1 11_29 0 772.81 1.3e-14 5.67
9811 RSBN1L 7_49 0 754.57 0.0e+00 1.29Genes with highest PVE
#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
5389 RPS11 19_34 1.000 163204.55 0.52000 17.58
7656 CATSPER2 15_16 0.993 303.50 0.00096 18.82
2173 TMEM176B 7_93 1.000 218.91 0.00069 -14.87
7523 SLC39A13 11_29 0.708 294.12 0.00066 -9.24
5789 ABRACL 6_92 1.000 180.93 0.00057 6.64
5799 SLC22A3 6_104 0.945 162.77 0.00049 -6.89
12008 HPR 16_38 1.000 130.98 0.00042 -9.93
12467 RP11-219B17.3 15_27 0.993 85.74 0.00027 -8.67
1231 PABPC4 1_24 0.987 79.71 0.00025 9.17
2809 COL7A1 3_35 0.982 81.43 0.00025 9.04
11885 SLC35G6 17_7 0.996 79.32 0.00025 -9.53
2998 RALGPS2 1_87 0.817 90.60 0.00023 9.36
6093 CSNK1G3 5_75 0.942 78.18 0.00023 9.49
2891 SNX17 2_16 0.523 125.51 0.00021 -12.86
5400 EPHA2 1_11 0.722 86.00 0.00020 -9.56
3212 CCND2 12_4 0.996 64.82 0.00020 -7.54
1144 ASAP3 1_16 1.000 55.42 0.00018 7.28
10085 RFX8 2_59 0.795 66.88 0.00017 -8.12
11558 LINC01184 5_78 0.903 57.96 0.00017 7.46
11237 TNF 6_25 0.937 56.02 0.00017 5.24Genes with largest z scores
#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
3804 PRRG2 19_34 0.000 23602.17 0.0e+00 -20.72
7656 CATSPER2 15_16 0.993 303.50 9.6e-04 18.82
5389 RPS11 19_34 1.000 163204.55 5.2e-01 17.58
1227 FLT3LG 19_34 0.000 141042.21 0.0e+00 -16.55
2173 TMEM176B 7_93 1.000 218.91 6.9e-04 -14.87
5393 RCN3 19_34 0.000 52574.94 0.0e+00 -14.77
7985 LCMT2 15_16 0.032 177.81 1.8e-05 -14.66
2887 NRBP1 2_16 0.015 183.23 8.9e-06 14.10
2891 SNX17 2_16 0.523 125.51 2.1e-04 -12.86
1058 GCKR 2_16 0.495 81.58 1.3e-04 -12.35
10987 C2orf16 2_16 0.495 81.58 1.3e-04 -12.35
7702 MAP1A 15_16 0.022 111.17 7.8e-06 12.19
9635 TLCD2 17_2 0.000 108.89 9.4e-12 11.75
8284 RBKS 2_16 0.020 108.93 6.9e-06 11.43
18 TMEM176A 7_93 0.004 147.13 2.0e-06 -11.43
3805 SCAF1 19_34 0.000 15951.61 0.0e+00 -11.38
9570 TMEM121 14_55 0.727 54.68 1.3e-04 11.31
11738 RP11-115J16.2 8_12 0.020 102.97 6.6e-06 11.06
8149 SERPINA6 14_48 0.000 186.68 1.2e-16 -10.93
3802 IRF3 19_34 0.000 15477.49 0.0e+00 -10.77Comparing z scores and PIPs
#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.02045684#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
3804 PRRG2 19_34 0.000 23602.17 0.0e+00 -20.72
7656 CATSPER2 15_16 0.993 303.50 9.6e-04 18.82
5389 RPS11 19_34 1.000 163204.55 5.2e-01 17.58
1227 FLT3LG 19_34 0.000 141042.21 0.0e+00 -16.55
2173 TMEM176B 7_93 1.000 218.91 6.9e-04 -14.87
5393 RCN3 19_34 0.000 52574.94 0.0e+00 -14.77
7985 LCMT2 15_16 0.032 177.81 1.8e-05 -14.66
2887 NRBP1 2_16 0.015 183.23 8.9e-06 14.10
2891 SNX17 2_16 0.523 125.51 2.1e-04 -12.86
1058 GCKR 2_16 0.495 81.58 1.3e-04 -12.35
10987 C2orf16 2_16 0.495 81.58 1.3e-04 -12.35
7702 MAP1A 15_16 0.022 111.17 7.8e-06 12.19
9635 TLCD2 17_2 0.000 108.89 9.4e-12 11.75
8284 RBKS 2_16 0.020 108.93 6.9e-06 11.43
18 TMEM176A 7_93 0.004 147.13 2.0e-06 -11.43
3805 SCAF1 19_34 0.000 15951.61 0.0e+00 -11.38
9570 TMEM121 14_55 0.727 54.68 1.3e-04 11.31
11738 RP11-115J16.2 8_12 0.020 102.97 6.6e-06 11.06
8149 SERPINA6 14_48 0.000 186.68 1.2e-16 -10.93
3802 IRF3 19_34 0.000 15477.49 0.0e+00 -10.77Locus plots for genes and SNPs
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: 19_34"
genename region_tag susie_pip mu2 PVE z
2042 BCAT2 19_34 0 6.81 0.00 -0.60
1110 HSD17B14 19_34 0 25.57 0.00 2.40
2044 PLEKHA4 19_34 0 35.00 0.00 -2.66
1921 NUCB1 19_34 0 18.47 0.00 2.34
1920 TULP2 19_34 0 28.50 0.00 -0.97
1922 DHDH 19_34 0 13.38 0.00 0.36
1113 FTL 19_34 0 151.04 0.00 2.42
9401 RUVBL2 19_34 0 31.02 0.00 -0.40
1928 SNRNP70 19_34 0 263.67 0.00 -1.94
1929 LIN7B 19_34 0 36.13 0.00 1.61
10994 C19orf73 19_34 0 137.74 0.00 -1.37
8899 PPFIA3 19_34 0 330.30 0.00 -1.73
4086 TRPM4 19_34 0 79.09 0.00 5.11
545 SLC6A16 19_34 0 1462.87 0.00 0.40
10291 CTC-301O7.4 19_34 0 3848.50 0.00 2.57
1940 SLC17A7 19_34 0 11229.71 0.00 -4.42
1932 PIH1D1 19_34 0 4923.97 0.00 3.78
6859 ALDH16A1 19_34 0 559.68 0.00 -0.18
1227 FLT3LG 19_34 0 141042.21 0.00 -16.55
5390 RPL13A 19_34 0 1127.45 0.00 -7.77
5389 RPS11 19_34 1 163204.55 0.52 17.58
1931 FCGRT 19_34 0 47948.16 0.00 -6.60
5393 RCN3 19_34 0 52574.94 0.00 -14.77
3804 PRRG2 19_34 0 23602.17 0.00 -20.72
5392 NOSIP 19_34 0 655.36 0.00 -8.69
3805 SCAF1 19_34 0 15951.61 0.00 -11.38
3802 IRF3 19_34 0 15477.49 0.00 -10.77
3803 PRMT1 19_34 0 16030.76 0.00 -9.25
8030 CPT1C 19_34 0 2092.69 0.00 0.75
3807 TSKS 19_34 0 71.11 0.00 2.31
10164 AP2A1 19_34 0 19.19 0.00 -0.39
162 FUZ 19_34 0 34.57 0.00 -1.30
1958 MED25 19_34 0 10.81 0.00 0.61
365 PNKP 19_34 0 79.91 0.00 0.28
1951 TBC1D17 19_34 0 40.20 0.00 -0.28
10797 NUP62 19_34 0 59.82 0.00 -1.12
8028 ATF5 19_34 0 296.42 0.00 -1.84
6860 SIGLEC11 19_34 0 161.38 0.00 0.09
5388 ZNF473 19_34 0 53.74 0.00 -3.32
1967 VRK3 19_34 0 72.87 0.00 -0.21
2009 MYH14 19_34 0 5.93 0.00 -0.24
4176 NR1H2 19_34 0 9.96 0.00 1.04
4174 KCNC3 19_34 0 5.38 0.00 0.08
4175 NAPSA 19_34 0 7.03 0.00 -0.43
543 POLD1 19_34 0 97.57 0.00 0.80
12177 SPIB 19_34 0 150.33 0.00 -2.18
1108 MYBPC2 19_34 0 14.87 0.00 0.32
10671 ASPDH 19_34 0 98.33 0.00 2.25
2030 CLEC11A 19_34 0 7.01 0.00 -0.58
7829 C19orf48 19_34 0 7.23 0.00 -1.12
9147 LINC01869 19_34 0 10.67 0.00 -1.51
7830 KLK1 19_34 0 54.43 0.00 0.50
4005 KLK10 19_34 0 20.96 0.00 -1.67
7831 KLK11 19_34 0 19.09 0.00 -1.48
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 15_16"
genename region_tag susie_pip mu2 PVE z
1853 ZNF106 15_16 0.071 41.35 9.3e-06 -5.45
9202 LRRC57 15_16 0.021 6.60 4.5e-07 0.10
6691 STARD9 15_16 0.017 5.03 2.7e-07 0.02
5189 CDAN1 15_16 0.023 7.47 5.4e-07 -0.27
3962 TTBK2 15_16 0.097 21.00 6.4e-06 -1.56
4903 TMEM62 15_16 0.081 23.95 6.2e-06 -2.15
7984 ADAL 15_16 0.020 42.74 2.7e-06 -7.91
7985 LCMT2 15_16 0.032 177.81 1.8e-05 -14.66
4898 TUBGCP4 15_16 0.016 59.49 3.1e-06 9.42
5180 ZSCAN29 15_16 0.016 13.45 6.9e-07 -4.45
7702 MAP1A 15_16 0.022 111.17 7.8e-06 12.19
7656 CATSPER2 15_16 0.993 303.50 9.6e-04 18.82
7709 PDIA3 15_16 0.196 40.59 2.5e-05 8.42
5178 MFAP1 15_16 0.017 65.76 3.5e-06 10.15
1286 WDR76 15_16 0.016 59.62 3.1e-06 9.48
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 7_93"
genename region_tag susie_pip mu2 PVE z
2166 ACTR3C 7_93 0.011 11.50 3.9e-07 -1.19
3854 LRRC61 7_93 0.020 16.47 1.0e-06 1.58
9945 ZBED6CL 7_93 0.008 9.44 2.3e-07 -1.25
2168 RARRES2 7_93 0.027 19.38 1.7e-06 1.83
10083 ZNF775 7_93 0.027 19.12 1.6e-06 -1.74
11468 LINC00996 7_93 0.004 5.54 7.3e-08 -0.09
8273 GIMAP8 7_93 0.008 8.86 2.3e-07 0.34
4368 GIMAP4 7_93 0.004 6.05 8.1e-08 -0.99
2172 GIMAP2 7_93 0.004 8.48 1.2e-07 1.94
10812 GIMAP1 7_93 0.007 8.60 1.8e-07 0.03
18 TMEM176A 7_93 0.004 147.13 2.0e-06 -11.43
2173 TMEM176B 7_93 1.000 218.91 6.9e-04 -14.87
14 AOC1 7_93 0.004 33.21 4.7e-07 -5.18
7377 NOS3 7_93 0.004 6.34 8.4e-08 1.24
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 2_16"
genename region_tag susie_pip mu2 PVE z
2881 CENPA 2_16 0.011 16.19 5.8e-07 -1.94
11149 OST4 2_16 0.008 11.12 2.8e-07 -1.76
4939 EMILIN1 2_16 0.006 14.99 2.9e-07 -3.45
4927 KHK 2_16 0.016 18.20 9.4e-07 -2.10
4935 PREB 2_16 0.007 14.60 3.3e-07 2.75
4941 ATRAID 2_16 0.012 46.26 1.7e-06 -5.81
4936 SLC5A6 2_16 0.010 47.30 1.6e-06 5.93
1060 CAD 2_16 0.010 33.08 1.0e-06 4.61
2885 SLC30A3 2_16 0.019 28.43 1.7e-06 5.26
7169 UCN 2_16 0.022 19.49 1.4e-06 0.18
2891 SNX17 2_16 0.523 125.51 2.1e-04 -12.86
7170 ZNF513 2_16 0.005 32.25 5.2e-07 4.44
2887 NRBP1 2_16 0.015 183.23 8.9e-06 14.10
4925 IFT172 2_16 0.005 14.70 2.4e-07 2.16
1058 GCKR 2_16 0.495 81.58 1.3e-04 -12.35
10987 C2orf16 2_16 0.495 81.58 1.3e-04 -12.35
10407 GPN1 2_16 0.062 71.57 1.4e-05 8.16
8847 CCDC121 2_16 0.005 10.55 1.8e-07 -1.94
6575 BRE 2_16 0.008 11.90 3.0e-07 -1.30
8284 RBKS 2_16 0.020 108.93 6.9e-06 11.43
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 17_2"
genename region_tag susie_pip mu2 PVE z
2368 YWHAE 17_2 0 8.42 1.2e-12 0.97
4256 INPP5K 17_2 0 12.05 3.2e-12 -2.41
8623 PITPNA 17_2 0 8.12 5.3e-13 0.16
7821 SLC43A2 17_2 0 9.87 8.4e-13 -0.29
7822 RILP 17_2 0 11.96 1.8e-12 -0.33
8621 PRPF8 17_2 0 10.08 1.1e-12 0.31
9635 TLCD2 17_2 0 108.89 9.4e-12 11.75
7824 WDR81 17_2 0 35.12 8.9e-12 5.78
9749 MIR22HG 17_2 0 90.33 5.2e-11 9.10
7823 SERPINF2 17_2 0 81.10 9.8e-12 0.47
4258 RPA1 17_2 0 6.08 3.5e-13 -1.76
9678 RTN4RL1 17_2 0 5.28 3.0e-13 1.13
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
SNPs with highest PIPs
#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
32178 rs72692783 1_74 1.000 45.16 1.4e-04 6.56
49581 rs79687284 1_108 1.000 50.87 1.6e-04 7.25
54976 rs766167074 1_118 1.000 2924.00 9.3e-03 2.76
60929 rs12239046 1_131 1.000 37.87 1.2e-04 -6.29
133630 rs12619647 2_144 1.000 66.08 2.1e-04 -8.25
165046 rs630850 3_67 1.000 43.23 1.4e-04 6.83
179791 rs9817452 3_97 1.000 73.38 2.3e-04 8.79
184889 rs234043 3_106 1.000 39.33 1.2e-04 6.19
193472 rs1406674 4_4 1.000 77.53 2.5e-04 -0.65
193483 rs3752442 4_4 1.000 145.55 4.6e-04 -13.82
193497 rs36205397 4_4 1.000 134.36 4.3e-04 15.48
216571 rs150783681 4_50 1.000 104.03 3.3e-04 -13.63
216649 rs72649167 4_50 1.000 95.80 3.0e-04 -13.42
225650 rs35518360 4_67 1.000 98.80 3.1e-04 -10.63
225716 rs13140033 4_68 1.000 63.57 2.0e-04 -8.35
271109 rs202048682 5_33 1.000 36.67 1.2e-04 5.75
287792 rs35676551 5_67 1.000 42.01 1.3e-04 -6.84
308911 rs70995354 6_2 1.000 1693.57 5.4e-03 3.41
308914 rs3929752 6_2 1.000 1661.91 5.3e-03 3.11
388146 rs761767938 7_49 1.000 7242.47 2.3e-02 -5.65
388154 rs1544459 7_49 1.000 7307.59 2.3e-02 -5.51
402905 rs125124 7_80 1.000 41.70 1.3e-04 -6.83
550217 rs2785172 11_23 1.000 52.23 1.7e-04 -7.58
557392 rs75592015 11_37 1.000 40.12 1.3e-04 6.31
558891 rs78488993 11_41 1.000 45.14 1.4e-04 6.03
596231 rs7397189 12_36 1.000 67.65 2.1e-04 -8.80
602295 rs11337960 12_47 1.000 2121.59 6.7e-03 -3.02
602299 rs1716526 12_47 1.000 2149.01 6.8e-03 -2.98
612892 rs141105880 12_67 1.000 62.38 2.0e-04 -6.69
662131 rs2883893 14_20 1.000 99.10 3.1e-04 -8.20
675956 rs1998057 14_48 1.000 420.58 1.3e-03 3.38
675957 rs12893029 14_48 1.000 638.89 2.0e-03 -2.06
675973 rs1243165 14_48 1.000 292.29 9.3e-04 7.87
739800 rs755736 17_29 1.000 84.20 2.7e-04 9.84
744552 rs113408695 17_39 1.000 202.99 6.4e-04 -17.59
784835 rs56361048 19_23 1.000 59.02 1.9e-04 -6.47
784837 rs889140 19_23 1.000 84.63 2.7e-04 -9.12
785429 rs1688031 19_24 1.000 338.19 1.1e-03 22.11
785430 rs4806075 19_24 1.000 436.24 1.4e-03 -9.46
785623 rs2251125 19_24 1.000 60.88 1.9e-04 -8.07
874706 rs11589479 1_76 1.000 138.44 4.4e-04 13.73
887479 rs1260326 2_16 1.000 366.57 1.2e-03 -23.45
1020877 rs2308005 6_92 1.000 428.84 1.4e-03 1.95
1023636 rs60425481 6_104 1.000 754.47 2.4e-03 -5.09
1074411 rs3072639 11_29 1.000 4088.30 1.3e-02 2.85
1140490 rs11078597 17_2 1.000 161.15 5.1e-04 20.74
1140604 rs118103201 17_2 1.000 98.95 3.1e-04 -7.20
1167319 rs113176985 19_34 1.000 150929.09 4.8e-01 -17.91
1167322 rs374141296 19_34 1.000 151976.27 4.8e-01 -16.04
34265 rs2446630 1_79 0.999 43.70 1.4e-04 7.77
109667 rs11689257 2_97 0.999 31.78 1.0e-04 5.62
554787 rs487579 11_32 0.999 38.41 1.2e-04 -5.02
722895 rs78839491 16_44 0.999 46.08 1.5e-04 -6.87
738492 rs2074292 17_27 0.999 41.32 1.3e-04 -7.63
753222 rs958079 18_7 0.999 33.84 1.1e-04 3.55
753232 rs2061135 18_7 0.999 41.34 1.3e-04 4.89
1142256 rs9901675 17_7 0.999 62.01 2.0e-04 -8.15
585641 rs66720652 12_15 0.998 31.96 1.0e-04 -5.38
34261 rs4657041 1_79 0.997 37.22 1.2e-04 6.91
193465 rs115019205 4_4 0.997 40.75 1.3e-04 5.88
254760 rs2736100 5_2 0.997 28.92 9.1e-05 -5.23
1037879 rs17256042 7_93 0.997 43.44 1.4e-04 0.42
51858 rs34179415 1_112 0.996 69.69 2.2e-04 7.90
767714 rs7237414 18_34 0.996 91.87 2.9e-04 10.12
687730 rs72743115 15_22 0.994 27.63 8.7e-05 -5.04
675394 rs11624512 14_46 0.993 95.57 3.0e-04 10.43
676065 rs2069987 14_48 0.993 65.73 2.1e-04 -8.24
287291 rs35552666 5_66 0.992 27.62 8.7e-05 -5.04
349861 rs212776 6_88 0.992 31.78 1.0e-04 5.63
557770 rs3740643 11_38 0.992 31.69 1.0e-04 5.63
173238 rs7649873 3_83 0.991 28.07 8.8e-05 -3.89
744518 rs4147967 17_39 0.991 41.54 1.3e-04 8.13
743039 rs12452590 17_36 0.990 26.95 8.5e-05 -4.98
388150 rs11972122 7_49 0.987 6754.21 2.1e-02 -5.14
707485 rs45545642 16_11 0.987 27.27 8.5e-05 5.04
780566 rs11668950 19_14 0.987 53.27 1.7e-04 7.10
307878 rs307812 5_109 0.985 26.74 8.4e-05 -4.95
553846 rs11040598 11_30 0.984 40.35 1.3e-04 6.64
541242 rs3750952 11_6 0.983 26.40 8.2e-05 -4.98
322605 rs78470916 6_32 0.982 30.45 9.5e-05 -5.45
416503 rs12543422 8_10 0.981 26.65 8.3e-05 4.96
694079 rs56357772 15_36 0.980 89.72 2.8e-04 -11.74
730763 rs402614 17_6 0.980 27.02 8.4e-05 4.91
184868 rs149368105 3_105 0.979 25.43 7.9e-05 4.81
368923 rs542176135 7_17 0.978 38.20 1.2e-04 4.75
276774 rs250722 5_45 0.976 32.94 1.0e-04 6.35
343440 rs519573 6_73 0.974 36.54 1.1e-04 -5.97
389550 rs3839804 7_51 0.972 30.34 9.4e-05 -5.81
451832 rs6982940 8_83 0.972 26.52 8.2e-05 -2.99
455084 rs2315839 8_88 0.970 30.87 9.5e-05 5.61
493386 rs1886296 9_73 0.969 25.35 7.8e-05 -4.79
240442 rs72727873 4_98 0.967 28.24 8.7e-05 5.31
51452 rs4846567 1_112 0.966 75.55 2.3e-04 -8.89
499945 rs10906857 10_13 0.966 25.53 7.8e-05 -4.87
556800 rs11227230 11_36 0.966 32.98 1.0e-04 -5.57
543133 rs67757744 11_10 0.964 25.36 7.8e-05 -4.80
800014 rs73605562 20_14 0.962 25.08 7.7e-05 4.64
145973 rs186945223 3_25 0.958 24.89 7.6e-05 4.70
785418 rs575494485 19_24 0.958 40.89 1.2e-04 -5.75
500902 rs58262664 10_14 0.957 29.72 9.0e-05 -5.45
241675 rs6816767 4_101 0.956 42.86 1.3e-04 -6.73
1139583 rs9905106 17_2 0.956 46.41 1.4e-04 -6.52
307565 rs10073754 5_108 0.954 26.12 7.9e-05 -5.05
784817 rs1423066 19_23 0.953 37.20 1.1e-04 -0.18
739195 rs8073834 17_27 0.950 25.35 7.6e-05 -3.54
248625 rs309704 4_114 0.947 92.03 2.8e-04 7.53
343036 rs2354558 6_71 0.939 25.12 7.5e-05 4.80
209386 rs4864786 4_39 0.936 24.86 7.4e-05 -3.08
615221 rs11064707 12_73 0.934 24.23 7.2e-05 4.62
306350 rs190940335 5_106 0.933 28.66 8.5e-05 5.38
318811 rs2246856 6_23 0.933 35.24 1.0e-04 -4.30
485902 rs10759266 9_54 0.931 29.20 8.6e-05 4.90
574489 rs1945396 11_75 0.931 45.45 1.3e-04 6.89
776106 rs191715972 19_5 0.929 24.95 7.4e-05 -4.72
1140494 rs62090014 17_2 0.928 82.51 2.4e-04 14.17
694105 rs72732561 15_36 0.926 38.49 1.1e-04 -3.71
341381 rs9496567 6_67 0.924 24.55 7.2e-05 -4.69
402347 rs3757387 7_78 0.923 24.93 7.3e-05 -4.65
97133 rs13027072 2_69 0.922 24.93 7.3e-05 3.74
97121 rs1562256 2_69 0.919 31.03 9.0e-05 4.51
138931 rs13085211 3_9 0.916 153.66 4.5e-04 -13.23
1167432 rs36013629 19_34 0.916 28916.88 8.4e-02 -26.15
479834 rs1360200 9_45 0.915 24.88 7.2e-05 4.78
193482 rs3748034 4_4 0.911 62.55 1.8e-04 10.54
110222 rs7607980 2_100 0.910 57.13 1.6e-04 -10.35
776475 rs146992497 19_6 0.897 23.49 6.7e-05 -4.43
110216 rs13389219 2_100 0.896 77.90 2.2e-04 -11.66
34304 rs114602456 1_79 0.889 25.71 7.3e-05 -4.57
394699 rs3197597 7_61 0.886 24.44 6.9e-05 -2.73
583248 rs3782567 12_12 0.882 47.65 1.3e-04 7.29
1109184 rs11621145 14_55 0.876 79.16 2.2e-04 -11.39
703651 rs2601781 16_4 0.871 27.44 7.6e-05 -5.05
387418 rs1179610 7_48 0.867 25.89 7.1e-05 4.62
489345 rs56141363 9_62 0.863 23.53 6.4e-05 4.33
325822 rs62406009 6_38 0.862 28.18 7.7e-05 4.49
372568 rs3757640 7_23 0.861 25.08 6.8e-05 4.55
694338 rs78190829 15_37 0.847 34.55 9.3e-05 5.88
267189 rs10214273 5_24 0.843 30.77 8.2e-05 5.51
460121 rs10970967 9_4 0.841 23.96 6.4e-05 -4.37
789606 rs12459419 19_35 0.837 29.46 7.8e-05 5.31
295012 rs145769851 5_82 0.835 25.45 6.7e-05 -4.68
869217 rs4970834 1_67 0.834 69.65 1.8e-04 -8.01
513606 rs1171830 10_39 0.832 42.76 1.1e-04 6.76
1119535 rs28364531 15_15 0.816 109.98 2.8e-04 11.09
666387 rs117940590 14_29 0.804 24.48 6.2e-05 -4.41
728851 rs12449600 17_3 0.801 61.19 1.6e-04 -7.97SNPs with largest effect sizes
#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
1167322 rs374141296 19_34 1 151976.3 0.48 -16.04
1167310 rs61371437 19_34 0 151122.8 0.00 -17.72
1167319 rs113176985 19_34 1 150929.1 0.48 -17.91
1167300 rs739349 19_34 0 150801.1 0.00 -17.69
1167301 rs756628 19_34 0 150798.2 0.00 -17.69
1167312 rs35295508 19_34 0 150679.4 0.00 -17.67
1167297 rs739347 19_34 0 150559.8 0.00 -17.68
1167298 rs2073614 19_34 0 150381.8 0.00 -17.65
1167326 rs2946865 19_34 0 150153.5 0.00 -17.88
1167293 rs4802613 19_34 0 150087.3 0.00 -17.66
1167317 rs73056069 19_34 0 150023.2 0.00 -17.68
1167303 rs2077300 19_34 0 149894.2 0.00 -17.56
1167314 rs2878354 19_34 0 149761.8 0.00 -17.55
1167307 rs73056059 19_34 0 149600.4 0.00 -17.60
1167291 rs10403394 19_34 0 149172.6 0.00 -17.58
1167292 rs17555056 19_34 0 149035.6 0.00 -17.56
1167327 rs60815603 19_34 0 147951.1 0.00 -17.88
1167330 rs1316885 19_34 0 147296.7 0.00 -18.01
1167335 rs2946863 19_34 0 147000.9 0.00 -18.03
1167332 rs60746284 19_34 0 146990.7 0.00 -17.79
1167328 rs35443645 19_34 0 146890.5 0.00 -17.98
1167308 rs73056062 19_34 0 145641.5 0.00 -16.98
1167338 rs553431297 19_34 0 143129.9 0.00 -17.22
1167321 rs112283514 19_34 0 142761.3 0.00 -16.69
1167323 rs11270139 19_34 0 141873.9 0.00 -16.94
1167288 rs10421294 19_34 0 134879.3 0.00 -17.39
1167287 rs8108175 19_34 0 134862.8 0.00 -17.39
1167280 rs59192944 19_34 0 134619.6 0.00 -17.37
1167286 rs1858742 19_34 0 134597.2 0.00 -17.35
1167277 rs55991145 19_34 0 134520.7 0.00 -17.36
1167272 rs3786567 19_34 0 134472.6 0.00 -17.35
1167268 rs2271952 19_34 0 134423.0 0.00 -17.35
1167271 rs4801801 19_34 0 134416.7 0.00 -17.36
1167267 rs2271953 19_34 0 134269.9 0.00 -17.37
1167269 rs2271951 19_34 0 134262.1 0.00 -17.36
1167258 rs60365978 19_34 0 134164.1 0.00 -17.40
1167264 rs4802612 19_34 0 133602.1 0.00 -17.26
1167274 rs2517977 19_34 0 133290.2 0.00 -17.35
1167261 rs55893003 19_34 0 133108.7 0.00 -17.07
1167253 rs55992104 19_34 0 130105.1 0.00 -16.71
1167247 rs60403475 19_34 0 130083.7 0.00 -16.72
1167250 rs4352151 19_34 0 130064.1 0.00 -16.72
1167244 rs11878448 19_34 0 129976.1 0.00 -16.71
1167238 rs9653100 19_34 0 129941.2 0.00 -16.70
1167234 rs4802611 19_34 0 129857.9 0.00 -16.69
1167226 rs7251338 19_34 0 129675.2 0.00 -16.67
1167225 rs59269605 19_34 0 129660.6 0.00 -16.67
1167246 rs1042120 19_34 0 129291.8 0.00 -16.60
1167242 rs113220577 19_34 0 129181.7 0.00 -16.59
1167236 rs9653118 19_34 0 128984.7 0.00 -16.60SNPs with highest PVE
#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
1167319 rs113176985 19_34 1.000 150929.09 0.48000 -17.91
1167322 rs374141296 19_34 1.000 151976.27 0.48000 -16.04
1167432 rs36013629 19_34 0.916 28916.88 0.08400 -26.15
388146 rs761767938 7_49 1.000 7242.47 0.02300 -5.65
388154 rs1544459 7_49 1.000 7307.59 0.02300 -5.51
388150 rs11972122 7_49 0.987 6754.21 0.02100 -5.14
1074411 rs3072639 11_29 1.000 4088.30 0.01300 2.85
54976 rs766167074 1_118 1.000 2924.00 0.00930 2.76
1167407 rs10419198 19_34 0.084 29418.97 0.00780 -26.01
602299 rs1716526 12_47 1.000 2149.01 0.00680 -2.98
602295 rs11337960 12_47 1.000 2121.59 0.00670 -3.02
308911 rs70995354 6_2 1.000 1693.57 0.00540 3.41
308914 rs3929752 6_2 1.000 1661.91 0.00530 3.11
54973 rs10489611 1_118 0.302 2899.15 0.00280 3.14
54967 rs2256908 1_118 0.290 2898.95 0.00270 3.14
54963 rs1076804 1_118 0.279 2895.37 0.00260 3.16
54974 rs2486737 1_118 0.271 2899.06 0.00250 3.13
1023636 rs60425481 6_104 1.000 754.47 0.00240 -5.09
54975 rs971534 1_118 0.241 2898.96 0.00220 3.13
54970 rs2790891 1_118 0.222 2898.71 0.00200 3.12
54971 rs2491405 1_118 0.222 2898.71 0.00200 3.12
675957 rs12893029 14_48 1.000 638.89 0.00200 -2.06
54985 rs1416913 1_118 0.180 2894.45 0.00160 3.13
54988 rs2790874 1_118 0.172 2894.16 0.00160 3.13
675965 rs760335 14_48 0.636 781.68 0.00160 14.28
785430 rs4806075 19_24 1.000 436.24 0.00140 -9.46
1020877 rs2308005 6_92 1.000 428.84 0.00140 1.95
1074394 rs1905285 11_29 0.107 4086.05 0.00140 2.98
675956 rs1998057 14_48 1.000 420.58 0.00130 3.38
887479 rs1260326 2_16 1.000 366.57 0.00120 -23.45
54982 rs2248646 1_118 0.115 2897.19 0.00110 3.09
54983 rs2211176 1_118 0.119 2897.33 0.00110 3.10
54984 rs2790882 1_118 0.119 2897.33 0.00110 3.10
785429 rs1688031 19_24 1.000 338.19 0.00110 22.11
1074413 rs7949513 11_29 0.087 4087.57 0.00110 2.95
1074489 rs12276188 11_29 0.078 4064.90 0.00100 3.11
1074467 rs7119161 11_29 0.077 4087.05 0.00099 2.95
1074475 rs7946068 11_29 0.076 4087.01 0.00098 2.95
54962 rs2474635 1_118 0.105 2880.23 0.00096 3.19
675973 rs1243165 14_48 1.000 292.29 0.00093 7.87
675964 rs2005945 14_48 0.364 780.71 0.00090 14.26
785427 rs2445819 19_24 0.469 594.14 0.00088 22.04
1074468 rs1872023 11_29 0.067 4080.81 0.00086 2.99
1074417 rs11039670 11_29 0.061 4087.52 0.00079 2.94
1074457 rs12295434 11_29 0.053 4086.19 0.00068 2.93
1074466 rs11039677 11_29 0.052 4086.17 0.00068 2.93
1074371 rs71045559 11_29 0.050 4085.76 0.00065 2.96
744552 rs113408695 17_39 1.000 202.99 0.00064 -17.59
785426 rs1688040 19_24 0.328 593.29 0.00062 22.03
1074449 rs7124318 11_29 0.041 4087.21 0.00053 2.93SNPs with largest z scores
#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
1167432 rs36013629 19_34 0.916 28916.88 8.4e-02 -26.15
1167407 rs10419198 19_34 0.084 29418.97 7.8e-03 -26.01
1167490 rs111476047 19_34 0.000 27404.90 0.0e+00 -23.87
887479 rs1260326 2_16 1.000 366.57 1.2e-03 -23.45
1167364 rs10469298 19_34 0.000 45799.20 0.0e+00 -22.88
1167377 rs1132990 19_34 0.000 46097.71 0.0e+00 -22.80
887501 rs780093 2_16 0.000 325.06 3.4e-07 -22.27
887266 rs4665972 2_16 0.000 321.63 2.6e-07 -22.26
887499 rs780094 2_16 0.000 323.66 3.1e-07 -22.25
1167340 rs2335534 19_34 0.000 75915.21 0.0e+00 -22.14
785429 rs1688031 19_24 1.000 338.19 1.1e-03 22.11
785427 rs2445819 19_24 0.469 594.14 8.8e-04 22.04
785426 rs1688040 19_24 0.328 593.29 6.2e-04 22.03
785425 rs2445818 19_24 0.203 592.41 3.8e-04 21.98
887530 rs11127048 2_16 0.000 313.58 4.3e-07 -21.79
887486 rs6547692 2_16 0.000 291.52 4.6e-07 -21.55
1167558 rs2379087 19_34 0.000 22013.70 0.0e+00 -20.86
1167614 rs10414643 19_34 0.000 22778.95 0.0e+00 -20.85
887516 rs1260333 2_16 0.000 272.51 3.6e-07 -20.83
1167585 rs10426059 19_34 0.000 20925.18 0.0e+00 -20.82
1167598 rs7249925 19_34 0.000 21679.71 0.0e+00 -20.82
887497 rs780096 2_16 0.000 270.58 3.1e-07 -20.80
887515 rs1260334 2_16 0.000 271.27 3.4e-07 -20.80
887520 rs1313566 2_16 0.000 271.20 3.4e-07 -20.80
1167588 rs112727702 19_34 0.000 21237.00 0.0e+00 -20.80
887498 rs780095 2_16 0.000 269.98 3.0e-07 -20.79
1167562 rs10416310 19_34 0.000 21884.49 0.0e+00 -20.79
1167593 rs2288921 19_34 0.000 21269.26 0.0e+00 -20.79
1167607 rs11878568 19_34 0.000 21697.94 0.0e+00 -20.79
1167629 rs2116922 19_34 0.000 23323.05 0.0e+00 -20.79
1167582 rs8113357 19_34 0.000 21262.00 0.0e+00 -20.78
1167627 rs10417980 19_34 0.000 23216.66 0.0e+00 -20.78
1167579 rs2890072 19_34 0.000 21348.57 0.0e+00 -20.76
1167590 rs10406941 19_34 0.000 21291.04 0.0e+00 -20.76
1167552 rs111310942 19_34 0.000 22098.58 0.0e+00 -20.75
1167578 rs2379088 19_34 0.000 21316.69 0.0e+00 -20.75
1140490 rs11078597 17_2 1.000 161.15 5.1e-04 20.74
1167572 rs10421333 19_34 0.000 21406.06 0.0e+00 -20.73
1167587 rs56873913 19_34 0.000 21279.90 0.0e+00 -20.73
1167568 rs3760708 19_34 0.000 21573.09 0.0e+00 -20.72
1167609 rs10404887 19_34 0.000 21242.53 0.0e+00 -20.72
1167610 rs7251295 19_34 0.000 21260.61 0.0e+00 -20.71
1167612 rs7251877 19_34 0.000 21146.83 0.0e+00 -20.69
1167584 rs762866768 19_34 0.000 21238.24 0.0e+00 -20.67
1167608 rs11083979 19_34 0.000 21271.89 0.0e+00 -20.67
1167615 rs16981329 19_34 0.000 21205.85 0.0e+00 -20.65
887524 rs2950835 2_16 0.000 266.67 3.2e-07 -20.64
887525 rs2911711 2_16 0.000 266.67 3.2e-07 -20.64
1167531 rs10412446 19_34 0.000 22412.28 0.0e+00 -20.63
1167550 rs7254718 19_34 0.000 22347.32 0.0e+00 -20.60Gene set enrichment for genes with PIP>0.8
#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] 55if (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)ZNF827 gene(s) from the input list not found in DisGeNET CURATEDRP11-219B17.3 gene(s) from the input list not found in DisGeNET CURATEDLRRC25 gene(s) from the input list not found in DisGeNET CURATEDRNF5 gene(s) from the input list not found in DisGeNET CURATEDPELO gene(s) from the input list not found in DisGeNET CURATEDRPS11 gene(s) from the input list not found in DisGeNET CURATEDGLMP gene(s) from the input list not found in DisGeNET CURATEDRALGPS2 gene(s) from the input list not found in DisGeNET CURATEDFAM53A gene(s) from the input list not found in DisGeNET CURATEDSH3BP1 gene(s) from the input list not found in DisGeNET CURATEDTMEM176B gene(s) from the input list not found in DisGeNET CURATEDCPEB2 gene(s) from the input list not found in DisGeNET CURATEDUBE2K gene(s) from the input list not found in DisGeNET CURATEDLINC00094 gene(s) from the input list not found in DisGeNET CURATEDRP11-131K5.2 gene(s) from the input list not found in DisGeNET CURATEDPCDH7 gene(s) from the input list not found in DisGeNET CURATEDCSNK1G3 gene(s) from the input list not found in DisGeNET CURATEDMOSPD3 gene(s) from the input list not found in DisGeNET CURATEDABRACL gene(s) from the input list not found in DisGeNET CURATEDJOSD1 gene(s) from the input list not found in DisGeNET CURATEDNPM2 gene(s) from the input list not found in DisGeNET CURATEDTNKS gene(s) from the input list not found in DisGeNET CURATEDMEF2B gene(s) from the input list not found in DisGeNET CURATEDHLA-DMB gene(s) from the input list not found in DisGeNET CURATEDGATSL3 gene(s) from the input list not found in DisGeNET CURATEDSCAF11 gene(s) from the input list not found in DisGeNET CURATEDSLC35G6 gene(s) from the input list not found in DisGeNET CURATEDHPR gene(s) from the input list not found in DisGeNET CURATEDLINC01184 gene(s) from the input list not found in DisGeNET CURATEDPALM3 gene(s) from the input list not found in DisGeNET CURATEDC14orf80 gene(s) from the input list not found in DisGeNET CURATEDRAVER2 gene(s) from the input list not found in DisGeNET CURATEDNYNRIN gene(s) from the input list not found in DisGeNET CURATED Description FDR Ratio
51 Coxsackievirus Infections 0.007324122 2/22
148 Myocarditis 0.007324122 2/22
188 Respiratory Tract Diseases 0.007324122 2/22
385 Carditis 0.007324122 2/22
395 Coronary Restenosis 0.007324122 2/22
23 Berylliosis 0.017145768 2/22
91 Hepatic Coma 0.017145768 2/22
92 Hepatic Encephalopathy 0.017145768 2/22
360 Fulminant Hepatic Failure with Cerebral Edema 0.017145768 2/22
361 Hepatic Stupor 0.017145768 2/22
BgRatio
51 6/9703
148 6/9703
188 6/9703
385 6/9703
395 5/9703
23 11/9703
91 13/9703
92 13/9703
360 13/9703
361 13/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