Albumin (quantile) - Whole_Blood
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-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-30620_irnt_Liver.Rmd
Modified: analysis/ukb-d-30630_irnt_Liver.Rmd
Modified: analysis/ukb-d-30640_irnt_Liver.Rmd
Modified: analysis/ukb-d-30650_irnt_Liver.Rmd
Modified: analysis/ukb-d-30660_irnt_Liver.Rmd
Modified: analysis/ukb-d-30670_irnt_Liver.Rmd
Modified: analysis/ukb-d-30680_irnt_Liver.Rmd
Modified: analysis/ukb-d-30690_irnt_Liver.Rmd
Modified: analysis/ukb-d-30700_irnt_Liver.Rmd
Modified: analysis/ukb-d-30710_irnt_Liver.Rmd
Modified: analysis/ukb-d-30720_irnt_Liver.Rmd
Modified: analysis/ukb-d-30730_irnt_Liver.Rmd
Modified: analysis/ukb-d-30740_irnt_Liver.Rmd
Modified: analysis/ukb-d-30750_irnt_Liver.Rmd
Modified: analysis/ukb-d-30760_irnt_Liver.Rmd
Modified: analysis/ukb-d-30770_irnt_Liver.Rmd
Modified: analysis/ukb-d-30780_irnt_Liver.Rmd
Modified: analysis/ukb-d-30790_irnt_Liver.Rmd
Modified: analysis/ukb-d-30800_irnt_Liver.Rmd
Modified: analysis/ukb-d-30810_irnt_Liver.Rmd
Modified: analysis/ukb-d-30820_irnt_Liver.Rmd
Modified: analysis/ukb-d-30830_irnt_Liver.Rmd
Modified: analysis/ukb-d-30840_irnt_Liver.Rmd
Modified: analysis/ukb-d-30850_irnt_Liver.Rmd
Modified: analysis/ukb-d-30860_irnt_Liver.Rmd
Modified: analysis/ukb-d-30870_irnt_Liver.Rmd
Modified: analysis/ukb-d-30880_irnt_Liver.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-30600_irnt_Whole_Blood.Rmd) and HTML (docs/ukb-d-30600_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 | 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 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-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 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])
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] 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.762776Load 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.0138851681 0.0001968373 #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
25.09011 12.93447 #report sample size
print(sample_size)[1] 315268#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.01226030 0.07023637 #compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))[1] 0.658995 1.298569Genes 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
8332 MGMT 10_81 1.000 73.91 2.3e-04 4.36
1980 FCGRT 19_34 1.000 197238.15 6.3e-01 17.88
6593 FCHO2 5_43 0.999 260.76 8.3e-04 16.28
7960 SERPINF2 17_2 0.997 265.74 8.4e-04 -15.96
4462 MPRIP 17_15 0.994 29.95 9.4e-05 -6.01
1483 TRIOBP 22_16 0.993 50.17 1.6e-04 -8.81
267 RUNX3 1_17 0.989 36.87 1.2e-04 -5.71
6206 ZNF827 4_95 0.988 30.43 9.5e-05 5.13
9863 LAMP1 13_62 0.988 32.17 1.0e-04 5.53
2342 CCDC6 10_39 0.976 45.62 1.4e-04 6.76
9181 BEND3 6_71 0.974 39.16 1.2e-04 -6.21
1870 XYLT1 16_16 0.972 28.03 8.6e-05 -4.69
5025 THBS1 15_13 0.971 27.98 8.6e-05 5.15
6299 PELO 5_31 0.969 36.08 1.1e-04 -5.83
9820 TLCD2 17_2 0.966 270.43 8.3e-04 -11.75
7786 CATSPER2 15_16 0.962 306.36 9.3e-04 18.82
1477 SH3BP1 22_16 0.957 49.27 1.5e-04 8.07
4458 ALOX5AP 13_9 0.953 25.47 7.7e-05 4.73
5095 DNAJC13 3_82 0.949 71.97 2.2e-04 5.54
1179 ASAP3 1_16 0.911 56.50 1.6e-04 5.94
236 TACC3 4_3 0.880 25.10 7.0e-05 4.67
3762 RAB17 2_140 0.869 22.12 6.1e-05 4.54
3347 IRF2BPL 14_36 0.850 26.31 7.1e-05 5.01
2072 TYK2 19_9 0.850 21.35 5.8e-05 3.63
2539 CBL 11_71 0.841 39.84 1.1e-04 -6.57
8815 ANGEL2 1_108 0.838 26.63 7.1e-05 4.85
2215 MEST 7_79 0.834 24.90 6.6e-05 4.54
7233 EOMES 3_20 0.830 27.32 7.2e-05 -4.85
4283 TRAF3 14_54 0.830 126.72 3.3e-04 -6.09
6495 NUP205 7_82 0.828 25.21 6.6e-05 4.46
7060 RAVER2 1_41 0.813 42.32 1.1e-04 6.31
10392 SND1 7_79 0.804 30.06 7.7e-05 5.53Genes 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
1980 FCGRT 19_34 1 197238.15 6.3e-01 17.88
5520 RCN3 19_34 0 62633.19 0.0e+00 14.88
8165 CPT1C 19_34 0 13536.39 0.0e+00 -5.30
4687 TMEM60 7_49 0 3832.78 0.0e+00 4.67
571 SLC6A16 19_34 0 3394.69 0.0e+00 -2.99
10492 CTC-301O7.4 19_34 0 3214.65 0.0e+00 -1.36
168 SPRTN 1_118 0 3209.34 1.8e-14 -2.83
3138 EXOC8 1_118 0 2324.78 0.0e+00 -3.19
11220 ADM5 19_34 0 2068.00 0.0e+00 -0.56
6980 ALDH16A1 19_34 0 2040.87 0.0e+00 5.05
846 TEAD2 19_34 0 1922.10 0.0e+00 -0.30
2551 PTPMT1 11_29 0 851.92 1.1e-14 5.51
11094 APTR 7_49 0 608.56 0.0e+00 -0.25
98 PHTF2 7_49 0 587.94 0.0e+00 0.66
4609 MYBPC3 11_29 0 567.40 1.2e-18 4.10
1989 CD37 19_34 0 504.64 0.0e+00 2.28
3140 TSNAX 1_118 0 455.05 0.0e+00 -1.73
5519 NOSIP 19_34 0 440.85 0.0e+00 -5.61
1230 GRAMD1A 19_24 0 388.63 0.0e+00 11.75
9074 PPFIA3 19_34 0 386.88 0.0e+00 -1.73Genes 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
1980 FCGRT 19_34 1.000 197238.15 0.63000 17.88
7786 CATSPER2 15_16 0.962 306.36 0.00093 18.82
7960 SERPINF2 17_2 0.997 265.74 0.00084 -15.96
6593 FCHO2 5_43 0.999 260.76 0.00083 16.28
9820 TLCD2 17_2 0.966 270.43 0.00083 -11.75
4283 TRAF3 14_54 0.830 126.72 0.00033 -6.09
8332 MGMT 10_81 1.000 73.91 0.00023 4.36
5095 DNAJC13 3_82 0.949 71.97 0.00022 5.54
5778 UCN2 3_34 0.788 83.32 0.00021 -9.21
1267 PABPC4 1_24 0.778 72.21 0.00018 8.75
9796 C14orf80 14_55 0.550 93.54 0.00016 11.28
1179 ASAP3 1_16 0.911 56.50 0.00016 5.94
1483 TRIOBP 22_16 0.993 50.17 0.00016 -8.81
1477 SH3BP1 22_16 0.957 49.27 0.00015 8.07
2342 CCDC6 10_39 0.976 45.62 0.00014 6.76
8293 CHRNB1 17_7 0.690 58.49 0.00013 9.29
267 RUNX3 1_17 0.989 36.87 0.00012 -5.71
9181 BEND3 6_71 0.974 39.16 0.00012 -6.21
7060 RAVER2 1_41 0.813 42.32 0.00011 6.31
6299 PELO 5_31 0.969 36.08 0.00011 -5.83Genes 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
7786 CATSPER2 15_16 0.962 306.36 9.3e-04 18.82
1980 FCGRT 19_34 1.000 197238.15 6.3e-01 17.88
6593 FCHO2 5_43 0.999 260.76 8.3e-04 16.28
7960 SERPINF2 17_2 0.997 265.74 8.4e-04 -15.96
5520 RCN3 19_34 0.000 62633.19 0.0e+00 14.88
7782 CASC4 15_17 0.031 202.63 2.0e-05 -14.86
2953 NRBP1 2_16 0.005 186.83 3.2e-06 14.14
2956 SNX17 2_16 0.008 173.74 4.3e-06 13.63
9820 TLCD2 17_2 0.966 270.43 8.3e-04 -11.75
1230 GRAMD1A 19_24 0.000 388.63 0.0e+00 11.75
9796 C14orf80 14_55 0.550 93.54 1.6e-04 11.28
9748 TMEM121 14_55 0.217 92.30 6.4e-05 11.24
8285 SERPINA6 14_49 0.000 187.41 9.0e-17 -10.93
2223 TMEM176B 7_94 0.014 84.59 3.7e-06 -10.21
5291 MFAP1 15_16 0.011 69.98 2.5e-06 10.15
2406 KAT2A 17_25 0.016 101.41 5.2e-06 9.94
2112 ISYNA1 19_15 0.059 82.25 1.5e-05 -9.36
8293 CHRNB1 17_7 0.690 58.49 1.3e-04 9.29
5778 UCN2 3_34 0.788 83.32 2.1e-04 -9.21
9945 MIR22HG 17_2 0.000 124.47 7.4e-11 -9.09Comparing 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.01919784#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
7786 CATSPER2 15_16 0.962 306.36 9.3e-04 18.82
1980 FCGRT 19_34 1.000 197238.15 6.3e-01 17.88
6593 FCHO2 5_43 0.999 260.76 8.3e-04 16.28
7960 SERPINF2 17_2 0.997 265.74 8.4e-04 -15.96
5520 RCN3 19_34 0.000 62633.19 0.0e+00 14.88
7782 CASC4 15_17 0.031 202.63 2.0e-05 -14.86
2953 NRBP1 2_16 0.005 186.83 3.2e-06 14.14
2956 SNX17 2_16 0.008 173.74 4.3e-06 13.63
9820 TLCD2 17_2 0.966 270.43 8.3e-04 -11.75
1230 GRAMD1A 19_24 0.000 388.63 0.0e+00 11.75
9796 C14orf80 14_55 0.550 93.54 1.6e-04 11.28
9748 TMEM121 14_55 0.217 92.30 6.4e-05 11.24
8285 SERPINA6 14_49 0.000 187.41 9.0e-17 -10.93
2223 TMEM176B 7_94 0.014 84.59 3.7e-06 -10.21
5291 MFAP1 15_16 0.011 69.98 2.5e-06 10.15
2406 KAT2A 17_25 0.016 101.41 5.2e-06 9.94
2112 ISYNA1 19_15 0.059 82.25 1.5e-05 -9.36
8293 CHRNB1 17_7 0.690 58.49 1.3e-04 9.29
5778 UCN2 3_34 0.788 83.32 2.1e-04 -9.21
9945 MIR22HG 17_2 0.000 124.47 7.4e-11 -9.09Locus 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: 15_16"
genename region_tag susie_pip mu2 PVE z
1325 SNAP23 15_16 0.015 6.98 3.4e-07 0.34
9382 LRRC57 15_16 0.014 6.42 2.9e-07 0.11
5030 HAUS2 15_16 0.051 33.74 5.5e-06 4.62
6785 STARD9 15_16 0.015 7.11 3.4e-07 0.04
5300 CDAN1 15_16 0.015 7.27 3.6e-07 -0.27
4064 TTBK2 15_16 0.076 21.81 5.2e-06 -1.56
7829 CCNDBP1 15_16 0.031 18.72 1.8e-06 1.91
1905 TGM5 15_16 0.020 36.47 2.4e-06 7.52
8115 ADAL 15_16 0.012 44.66 1.7e-06 -7.91
8116 LCMT2 15_16 0.012 44.66 1.7e-06 -7.91
5034 TUBGCP4 15_16 0.012 6.05 2.2e-07 -0.56
5295 ZSCAN29 15_16 0.011 12.16 4.3e-07 -4.18
7831 MAP1A 15_16 0.026 34.44 2.8e-06 8.03
7786 CATSPER2 15_16 0.962 306.36 9.3e-04 18.82
7839 PDIA3 15_16 0.045 23.78 3.4e-06 -5.66
4065 ELL3 15_16 0.068 21.22 4.6e-06 -2.96
5294 SERF2 15_16 0.068 21.22 4.6e-06 -2.96
5291 MFAP1 15_16 0.011 69.98 2.5e-06 10.15
1323 WDR76 15_16 0.139 31.74 1.4e-05 -3.62
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 19_34"
genename region_tag susie_pip mu2 PVE z
2098 BCAT2 19_34 0 41.77 0.00 1.77
1143 HSD17B14 19_34 0 16.94 0.00 2.43
2100 PLEKHA4 19_34 0 101.72 0.00 2.70
1142 PPP1R15A 19_34 0 96.58 0.00 2.70
1967 NUCB1 19_34 0 12.38 0.00 1.49
1968 DHDH 19_34 0 32.99 0.00 0.70
1146 FTL 19_34 0 161.82 0.00 -2.42
1969 GYS1 19_34 0 32.92 0.00 0.33
9576 RUVBL2 19_34 0 24.81 0.00 -0.46
12166 CTB-60B18.10 19_34 0 54.64 0.00 -0.35
1978 LIN7B 19_34 0 10.68 0.00 -0.60
1976 SNRNP70 19_34 0 365.08 0.00 -2.04
11184 C19orf73 19_34 0 148.02 0.00 -1.37
9074 PPFIA3 19_34 0 386.88 0.00 -1.73
4200 TRPM4 19_34 0 109.20 0.00 4.74
571 SLC6A16 19_34 0 3394.69 0.00 -2.99
10492 CTC-301O7.4 19_34 0 3214.65 0.00 -1.36
1989 CD37 19_34 0 504.64 0.00 2.28
846 TEAD2 19_34 0 1922.10 0.00 -0.30
6980 ALDH16A1 19_34 0 2040.87 0.00 5.05
1980 FCGRT 19_34 1 197238.15 0.63 17.88
5520 RCN3 19_34 0 62633.19 0.00 14.88
5519 NOSIP 19_34 0 440.85 0.00 -5.61
11220 ADM5 19_34 0 2068.00 0.00 -0.56
8165 CPT1C 19_34 0 13536.39 0.00 -5.30
3903 TSKS 19_34 0 108.35 0.00 3.13
10357 AP2A1 19_34 0 36.13 0.00 1.09
181 FUZ 19_34 0 15.34 0.00 -1.86
2006 MED25 19_34 0 8.18 0.00 0.61
2001 PTOV1 19_34 0 5.42 0.00 -0.51
381 PNKP 19_34 0 56.31 0.00 -0.17
1998 TBC1D17 19_34 0 24.42 0.00 1.56
8162 ATF5 19_34 0 165.78 0.00 -3.38
10990 NUP62 19_34 0 10.51 0.00 -1.34
6983 SIGLEC11 19_34 0 52.33 0.00 2.39
2015 VRK3 19_34 0 294.59 0.00 2.56
5516 ZNF473 19_34 0 55.00 0.00 3.46
2061 MYH14 19_34 0 203.89 0.00 -0.06
4295 NR1H2 19_34 0 11.74 0.00 1.04
4294 NAPSA 19_34 0 5.73 0.00 -0.43
6988 EMC10 19_34 0 17.49 0.00 -0.41
6989 JOSD2 19_34 0 21.40 0.00 0.22
10859 ASPDH 19_34 0 42.87 0.00 -0.11
2083 CLEC11A 19_34 0 6.82 0.00 -0.59
7965 C19orf48 19_34 0 20.86 0.00 -0.93
9327 LINC01869 19_34 0 6.57 0.00 -1.09
7966 KLK1 19_34 0 10.25 0.00 0.91
8152 KLK7 19_34 0 14.30 0.00 1.34
7968 KLK11 19_34 0 16.07 0.00 0.29
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 5_43"
genename region_tag susie_pip mu2 PVE z
4316 MAP1B 5_43 0.002 5.83 4.3e-08 -0.56
2789 MRPS27 5_43 0.006 12.37 2.3e-07 0.07
1055 TNPO1 5_43 0.002 79.37 5.6e-07 -8.86
6593 FCHO2 5_43 0.999 260.76 8.3e-04 16.28
6595 TMEM171 5_43 0.002 25.27 1.9e-07 -4.24
5832 BTF3 5_43 0.003 7.26 5.8e-08 0.94
7436 UTP15 5_43 0.004 10.12 1.3e-07 0.29
7434 ANKRA2 5_43 0.002 5.36 4.1e-08 -0.50
11122 ARHGEF28 5_43 0.006 23.87 4.3e-07 -4.07
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 17_2"
genename region_tag susie_pip mu2 PVE z
7864 CRK 17_2 0.000 8.96 1.8e-12 2.01
10502 MYO1C 17_2 0.000 7.23 3.8e-13 -0.46
4379 INPP5K 17_2 0.000 12.18 2.7e-12 -2.45
8783 PITPNA 17_2 0.000 6.25 5.5e-13 -1.33
7958 SLC43A2 17_2 0.000 8.76 5.3e-13 0.37
7959 RILP 17_2 0.000 13.46 1.5e-12 -0.33
8781 PRPF8 17_2 0.000 16.75 1.3e-12 -3.42
9820 TLCD2 17_2 0.966 270.43 8.3e-04 -11.75
9945 MIR22HG 17_2 0.000 124.47 7.4e-11 -9.09
7961 WDR81 17_2 0.000 91.45 2.9e-11 -0.55
7960 SERPINF2 17_2 0.997 265.74 8.4e-04 -15.96
4382 SERPINF1 17_2 0.000 42.62 9.2e-11 2.86
9939 SMYD4 17_2 0.000 31.31 1.6e-11 2.67
4381 RPA1 17_2 0.000 6.27 3.4e-13 -0.32
| Version | Author | Date |
|---|---|---|
| dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 15_17"
genename region_tag susie_pip mu2 PVE z
7782 CASC4 15_17 0.031 202.63 2.0e-05 -14.86
11335 PATL2 15_17 0.018 7.35 4.2e-07 1.54
7780 B2M 15_17 0.017 5.69 3.0e-07 -0.67
9861 TRIM69 15_17 0.020 7.70 4.9e-07 1.07
5293 SORD 15_17 0.025 9.65 7.8e-07 1.28
5042 DUOX1 15_17 0.017 6.01 3.3e-07 -0.57
8498 GATM 15_17 0.027 10.21 8.9e-07 -1.07
8497 SPATA5L1 15_17 0.016 5.05 2.6e-07 0.05
12481 RP11-96O20.5 15_17 0.028 10.44 9.2e-07 1.28
5023 SQRDL 15_17 0.074 20.61 4.8e-06 -2.17
12408 CTD-2306A12.1 15_17 0.019 6.20 3.6e-07 0.28
| 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
5089 rs4336844 1_11 1.000 78.68 2.5e-04 9.46
34107 rs72692783 1_74 1.000 44.42 1.4e-04 6.56
57124 rs766167074 1_118 1.000 3177.23 1.0e-02 2.76
63077 rs12239046 1_131 1.000 38.10 1.2e-04 -6.29
72803 rs780093 2_16 1.000 418.22 1.3e-03 -22.27
135890 rs12619647 2_144 1.000 66.77 2.1e-04 -8.25
169048 rs630850 3_67 1.000 43.55 1.4e-04 6.83
183675 rs9817452 3_97 1.000 73.90 2.3e-04 8.79
188773 rs234043 3_106 1.000 38.58 1.2e-04 6.19
197818 rs1406674 4_4 1.000 79.04 2.5e-04 -0.65
197829 rs3752442 4_4 1.000 146.75 4.7e-04 -13.82
197843 rs36205397 4_4 1.000 134.82 4.3e-04 15.48
221826 rs150783681 4_52 1.000 104.84 3.3e-04 -13.63
221904 rs72649167 4_52 1.000 96.13 3.0e-04 -13.42
230905 rs35518360 4_67 1.000 99.38 3.2e-04 -10.63
230971 rs13140033 4_68 1.000 63.88 2.0e-04 -8.35
291889 rs35676551 5_67 1.000 42.23 1.3e-04 -6.84
314542 rs70995354 6_2 1.000 1778.14 5.6e-03 3.41
314545 rs3929752 6_2 1.000 1776.42 5.6e-03 3.11
364229 rs56393506 6_104 1.000 97.82 3.1e-04 -5.06
369004 rs73041368 7_4 1.000 93.09 3.0e-04 -4.45
395054 rs761767938 7_49 1.000 7846.85 2.5e-02 -5.65
395062 rs1544459 7_49 1.000 7913.30 2.5e-02 -5.51
409168 rs125124 7_80 1.000 41.92 1.3e-04 -6.83
428715 rs118162691 8_23 1.000 56.40 1.8e-04 -7.72
558174 rs2785172 11_23 1.000 52.48 1.7e-04 -7.58
565349 rs75592015 11_37 1.000 39.07 1.2e-04 6.31
566848 rs78488993 11_41 1.000 46.87 1.5e-04 6.03
604473 rs7397189 12_36 1.000 69.29 2.2e-04 -8.80
610537 rs11337960 12_47 1.000 2287.39 7.3e-03 -3.02
610541 rs1716526 12_47 1.000 2315.91 7.3e-03 -2.98
621134 rs141105880 12_67 1.000 62.88 2.0e-04 -6.69
669780 rs2883893 14_20 1.000 101.35 3.2e-04 -8.20
682826 rs1998057 14_49 1.000 437.16 1.4e-03 3.38
682827 rs12893029 14_49 1.000 662.62 2.1e-03 -2.06
682843 rs1243165 14_49 1.000 299.19 9.5e-04 7.87
748410 rs755736 17_29 1.000 84.33 2.7e-04 9.84
753462 rs113408695 17_39 1.000 204.06 6.5e-04 -17.59
793369 rs56361048 19_23 1.000 60.12 1.9e-04 -6.47
793371 rs889140 19_23 1.000 85.89 2.7e-04 -9.12
859645 rs11589479 1_77 1.000 149.09 4.7e-04 13.73
951239 rs10456852 6_71 1.000 49.58 1.6e-04 7.11
984004 rs201524046 10_81 1.000 201.07 6.4e-04 -2.97
988384 rs3072639 11_29 1.000 4586.40 1.5e-02 2.85
1030751 rs56282047 14_54 1.000 4467.07 1.4e-02 3.56
1030761 rs942021 14_54 1.000 4466.86 1.4e-02 3.76
1051010 rs11078597 17_2 1.000 383.38 1.2e-03 20.74
1052776 rs9901675 17_7 1.000 62.09 2.0e-04 -8.15
1075698 rs41523449 19_24 1.000 1748.62 5.5e-03 23.52
1075702 rs749726391 19_24 1.000 2479.61 7.9e-03 -5.48
1075703 rs12461158 19_24 1.000 2393.53 7.6e-03 -7.85
1081522 rs374141296 19_34 1.000 184631.64 5.9e-01 -16.04
112398 rs11689257 2_97 0.999 32.23 1.0e-04 5.62
562744 rs487579 11_32 0.999 40.07 1.3e-04 -5.02
731104 rs78839491 16_44 0.999 46.51 1.5e-04 -6.87
762132 rs958079 18_7 0.999 34.19 1.1e-04 3.55
762142 rs2061135 18_7 0.999 41.73 1.3e-04 4.89
937143 rs2734335 6_26 0.999 57.24 1.8e-04 6.60
36051 rs2446630 1_79 0.998 43.90 1.4e-04 7.77
593437 rs66720652 12_15 0.998 32.20 1.0e-04 -5.38
197811 rs115019205 4_4 0.997 41.08 1.3e-04 5.88
260015 rs2736100 5_2 0.997 29.13 9.2e-05 -5.23
727984 rs57186116 16_38 0.997 36.36 1.1e-04 4.94
54006 rs34179415 1_112 0.996 69.66 2.2e-04 7.90
776624 rs7237414 18_34 0.996 92.84 2.9e-04 10.12
36047 rs4657041 1_79 0.995 37.91 1.2e-04 6.91
694678 rs72743115 15_22 0.995 27.85 8.8e-05 -5.04
747102 rs2074292 17_27 0.994 41.49 1.3e-04 -7.63
176716 rs7649873 3_83 0.993 28.43 9.0e-05 -3.89
355539 rs212776 6_88 0.993 32.18 1.0e-04 5.63
682264 rs11624512 14_46 0.993 94.88 3.0e-04 10.43
682935 rs2069987 14_49 0.993 66.15 2.1e-04 -8.24
701952 rs72732561 15_36 0.993 39.44 1.2e-04 -3.71
565727 rs3740643 11_38 0.992 31.52 9.9e-05 5.63
291388 rs35552666 5_66 0.991 27.61 8.7e-05 -5.04
751649 rs12452590 17_36 0.991 27.14 8.5e-05 -4.98
753428 rs4147967 17_39 0.991 41.69 1.3e-04 8.13
733436 rs2255451 16_49 0.989 33.57 1.1e-04 -5.73
327236 rs78470916 6_32 0.988 31.16 9.8e-05 -5.45
395058 rs11972122 7_49 0.988 7311.48 2.3e-02 -5.14
789100 rs11668950 19_14 0.986 53.35 1.7e-04 7.10
561803 rs11040598 11_30 0.985 40.78 1.3e-04 6.64
421901 rs12543422 8_10 0.981 26.80 8.3e-05 4.96
701926 rs56357772 15_36 0.981 91.46 2.8e-04 -11.74
564757 rs11227230 11_36 0.980 35.29 1.1e-04 -5.57
188752 rs149368105 3_105 0.979 25.51 7.9e-05 4.81
375831 rs542176135 7_17 0.979 38.65 1.2e-04 4.75
458524 rs6982940 8_83 0.979 27.15 8.4e-05 -2.99
715332 rs45545642 16_12 0.979 27.82 8.6e-05 5.04
549199 rs3750952 11_7 0.976 26.65 8.2e-05 -4.98
348124 rs519573 6_73 0.975 36.71 1.1e-04 -5.97
813305 rs11698812 20_28 0.974 27.93 8.6e-05 -5.55
396458 rs3839804 7_51 0.973 30.48 9.4e-05 -5.81
738972 rs402614 17_6 0.972 26.85 8.3e-05 4.91
502099 rs1886296 9_73 0.970 25.45 7.8e-05 -4.79
955803 rs62621812 7_79 0.970 33.05 1.0e-04 6.00
461776 rs2315839 8_88 0.969 30.82 9.5e-05 5.61
364148 rs60425481 6_104 0.968 39.47 1.2e-04 -5.09
508658 rs10906857 10_13 0.968 25.69 7.9e-05 -4.87
53600 rs4846567 1_112 0.966 76.15 2.3e-04 -8.89
245697 rs72727873 4_98 0.965 28.98 8.9e-05 5.31
551090 rs67757744 11_10 0.965 25.44 7.8e-05 -4.80
807788 rs73605562 20_14 0.965 25.18 7.7e-05 4.64
793351 rs1423066 19_23 0.962 38.62 1.2e-04 -0.18
753272 rs3072826 17_39 0.961 64.78 2.0e-04 8.06
159187 rs7625774 3_48 0.960 42.03 1.3e-04 6.75
148233 rs186945223 3_25 0.959 24.95 7.6e-05 4.70
509615 rs58262664 10_14 0.959 29.83 9.1e-05 -5.45
747805 rs8073834 17_27 0.958 25.55 7.8e-05 -3.54
1050103 rs9905106 17_2 0.958 46.42 1.4e-04 -6.52
246930 rs6816767 4_101 0.957 43.14 1.3e-04 -6.73
211197 rs58932203 4_32 0.956 28.65 8.7e-05 -5.21
253880 rs309704 4_114 0.949 93.20 2.8e-04 7.53
623463 rs11064707 12_73 0.949 24.25 7.3e-05 4.62
280871 rs250722 5_45 0.948 30.61 9.2e-05 6.35
785016 rs191715972 19_5 0.948 25.03 7.5e-05 -4.72
566796 rs7944225 11_41 0.947 41.05 1.2e-04 5.41
1081632 rs36013629 19_34 0.947 35230.23 1.1e-01 -26.15
938579 rs368213214 6_26 0.945 81.93 2.5e-04 -8.41
214641 rs4864786 4_39 0.941 25.18 7.5e-05 -3.08
1051014 rs62090014 17_2 0.938 113.09 3.4e-04 14.17
311981 rs190940335 5_106 0.937 28.64 8.5e-05 5.38
691990 rs598899 15_15 0.935 24.63 7.3e-05 -5.18
324442 rs2246856 6_23 0.934 35.64 1.1e-04 -4.30
493122 rs10759266 9_54 0.933 29.41 8.7e-05 4.90
582039 rs1945396 11_75 0.931 45.69 1.3e-04 6.89
858009 rs4970834 1_67 0.931 68.90 2.0e-04 -8.01
346012 rs9496567 6_67 0.930 24.54 7.2e-05 -4.69
99864 rs13027072 2_69 0.928 25.04 7.4e-05 3.74
401607 rs3197597 7_61 0.927 24.65 7.2e-05 -2.73
141191 rs13085211 3_9 0.917 156.27 4.5e-04 -13.23
112953 rs7607980 2_100 0.915 58.36 1.7e-04 -10.35
591044 rs3782567 12_12 0.914 45.23 1.3e-04 7.29
197828 rs3748034 4_4 0.912 62.64 1.8e-04 10.54
112947 rs13389219 2_100 0.904 80.59 2.3e-04 -11.66
401620 rs6974537 7_62 0.903 24.62 7.0e-05 -4.77
785385 rs146992497 19_6 0.893 23.61 6.7e-05 -4.43
313196 rs10073754 5_108 0.892 26.76 7.6e-05 -5.05
364202 rs141383002 6_104 0.882 28.22 7.9e-05 2.99
353584 rs7763784 6_84 0.881 40.48 1.1e-04 3.46
487054 rs1360200 9_45 0.877 26.07 7.3e-05 4.78
711498 rs2601781 16_4 0.872 27.45 7.6e-05 -5.05
152353 rs116643069 3_35 0.868 55.20 1.5e-04 7.55
300176 rs145769851 5_82 0.865 25.08 6.9e-05 -4.68
394326 rs1179610 7_48 0.865 25.94 7.1e-05 4.62
697288 rs340019 15_27 0.856 86.73 2.4e-04 8.28
669806 rs12587842 14_20 0.849 33.42 9.0e-05 3.36
466813 rs10970967 9_4 0.848 23.88 6.4e-05 -4.37
272444 rs10214273 5_24 0.846 31.03 8.3e-05 5.51
313509 rs307812 5_109 0.846 26.91 7.2e-05 -4.95
99852 rs1562256 2_69 0.845 31.89 8.5e-05 4.51
36090 rs114602456 1_79 0.843 27.83 7.4e-05 -4.57
330453 rs62406009 6_39 0.842 28.53 7.6e-05 4.49
797380 rs12459419 19_35 0.838 29.67 7.9e-05 5.31
69911 rs6531051 2_10 0.833 25.73 6.8e-05 4.54
702185 rs78190829 15_37 0.833 34.24 9.0e-05 5.88
497382 rs56141363 9_62 0.828 24.19 6.4e-05 4.33
669772 rs145991799 14_20 0.823 32.26 8.4e-05 0.91
72792 rs373142409 2_16 0.821 29.63 7.7e-05 -5.96
36092 rs10917685 1_79 0.815 35.39 9.1e-05 -3.97
13029 rs977747 1_29 0.809 31.57 8.1e-05 4.90
369029 rs10232355 7_4 0.808 54.50 1.4e-04 -8.03SNPs 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
1081522 rs374141296 19_34 1 184631.6 0.59 -16.04
1081519 rs113176985 19_34 0 183770.3 0.00 -17.91
1081526 rs2946865 19_34 0 183761.4 0.00 -17.88
1081512 rs35295508 19_34 0 183352.3 0.00 -17.67
1081510 rs61371437 19_34 0 182840.1 0.00 -17.72
1081517 rs73056069 19_34 0 182798.0 0.00 -17.68
1081514 rs2878354 19_34 0 182236.2 0.00 -17.55
1081500 rs739349 19_34 0 182082.3 0.00 -17.69
1081501 rs756628 19_34 0 182082.3 0.00 -17.69
1081497 rs739347 19_34 0 181714.2 0.00 -17.68
1081498 rs2073614 19_34 0 181485.5 0.00 -17.65
1081503 rs2077300 19_34 0 180992.5 0.00 -17.56
1081493 rs4802613 19_34 0 180679.1 0.00 -17.66
1081507 rs73056059 19_34 0 180679.1 0.00 -17.60
1081527 rs60815603 19_34 0 180375.9 0.00 -17.88
1081530 rs1316885 19_34 0 180275.9 0.00 -18.01
1081535 rs2946863 19_34 0 179975.8 0.00 -18.03
1081528 rs35443645 19_34 0 179669.7 0.00 -17.98
1081532 rs60746284 19_34 0 179262.8 0.00 -17.79
1081491 rs10403394 19_34 0 178229.6 0.00 -17.58
1081492 rs17555056 19_34 0 178102.7 0.00 -17.56
1081508 rs73056062 19_34 0 175964.6 0.00 -16.98
1081538 rs553431297 19_34 0 173973.4 0.00 -17.22
1081521 rs112283514 19_34 0 173297.6 0.00 -16.69
1081523 rs11270139 19_34 0 172305.5 0.00 -16.94
1081488 rs10421294 19_34 0 161149.4 0.00 -17.39
1081487 rs8108175 19_34 0 161127.3 0.00 -17.39
1081480 rs59192944 19_34 0 160818.2 0.00 -17.37
1081486 rs1858742 19_34 0 160759.0 0.00 -17.35
1081477 rs55991145 19_34 0 160679.7 0.00 -17.36
1081472 rs3786567 19_34 0 160615.5 0.00 -17.35
1081468 rs2271952 19_34 0 160550.5 0.00 -17.35
1081471 rs4801801 19_34 0 160531.1 0.00 -17.36
1081467 rs2271953 19_34 0 160352.3 0.00 -17.37
1081469 rs2271951 19_34 0 160344.0 0.00 -17.36
1081458 rs60365978 19_34 0 160231.1 0.00 -17.40
1081464 rs4802612 19_34 0 159574.9 0.00 -17.26
1081474 rs2517977 19_34 0 159343.6 0.00 -17.35
1081461 rs55893003 19_34 0 159132.0 0.00 -17.07
1081453 rs55992104 19_34 0 155356.0 0.00 -16.71
1081447 rs60403475 19_34 0 155339.0 0.00 -16.72
1081450 rs4352151 19_34 0 155288.7 0.00 -16.72
1081444 rs11878448 19_34 0 155175.9 0.00 -16.71
1081438 rs9653100 19_34 0 155143.7 0.00 -16.70
1081434 rs4802611 19_34 0 155041.6 0.00 -16.69
1081426 rs7251338 19_34 0 154801.9 0.00 -16.67
1081425 rs59269605 19_34 0 154783.5 0.00 -16.67
1081446 rs1042120 19_34 0 154372.0 0.00 -16.60
1081442 rs113220577 19_34 0 154234.9 0.00 -16.59
1081436 rs9653118 19_34 0 153995.2 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
1081522 rs374141296 19_34 1.000 184631.64 0.59000 -16.04
1081632 rs36013629 19_34 0.947 35230.23 0.11000 -26.15
395054 rs761767938 7_49 1.000 7846.85 0.02500 -5.65
395062 rs1544459 7_49 1.000 7913.30 0.02500 -5.51
395058 rs11972122 7_49 0.988 7311.48 0.02300 -5.14
988384 rs3072639 11_29 1.000 4586.40 0.01500 2.85
1030751 rs56282047 14_54 1.000 4467.07 0.01400 3.56
1030761 rs942021 14_54 1.000 4466.86 0.01400 3.76
57124 rs766167074 1_118 1.000 3177.23 0.01000 2.76
1075702 rs749726391 19_24 1.000 2479.61 0.00790 -5.48
1075703 rs12461158 19_24 1.000 2393.53 0.00760 -7.85
610537 rs11337960 12_47 1.000 2287.39 0.00730 -3.02
610541 rs1716526 12_47 1.000 2315.91 0.00730 -2.98
1081607 rs10419198 19_34 0.053 35845.69 0.00610 -26.01
314542 rs70995354 6_2 1.000 1778.14 0.00560 3.41
314545 rs3929752 6_2 1.000 1776.42 0.00560 3.11
1075698 rs41523449 19_24 1.000 1748.62 0.00550 23.52
1075737 rs45512696 19_24 0.795 2003.90 0.00510 24.77
1030777 rs11622216 14_54 0.251 4382.67 0.00350 4.07
57121 rs10489611 1_118 0.277 3207.33 0.00280 3.14
57115 rs2256908 1_118 0.257 3207.15 0.00260 3.14
57122 rs2486737 1_118 0.232 3207.19 0.00240 3.13
57111 rs1076804 1_118 0.206 3203.40 0.00210 3.16
682827 rs12893029 14_49 1.000 662.62 0.00210 -2.06
57123 rs971534 1_118 0.192 3207.05 0.00200 3.13
57118 rs2790891 1_118 0.167 3206.83 0.00170 3.12
57119 rs2491405 1_118 0.167 3206.83 0.00170 3.12
988367 rs1905285 11_29 0.114 4583.89 0.00170 2.98
682835 rs760335 14_49 0.638 801.92 0.00160 14.28
682826 rs1998057 14_49 1.000 437.16 0.00140 3.38
988386 rs7949513 11_29 0.096 4585.71 0.00140 2.95
72803 rs780093 2_16 1.000 418.22 0.00130 -22.27
988440 rs7119161 11_29 0.084 4585.13 0.00120 2.95
988448 rs7946068 11_29 0.083 4585.08 0.00120 2.95
1051010 rs11078597 17_2 1.000 383.38 0.00120 20.74
57133 rs1416913 1_118 0.096 3201.98 0.00098 3.13
988390 rs11039670 11_29 0.066 4585.67 0.00096 2.94
988441 rs1872023 11_29 0.066 4577.96 0.00096 2.99
682843 rs1243165 14_49 1.000 299.19 0.00095 7.87
57136 rs2790874 1_118 0.090 3201.66 0.00092 3.13
682834 rs2005945 14_49 0.362 800.94 0.00092 14.26
988462 rs12276188 11_29 0.063 4559.66 0.00092 3.11
1075732 rs149131600 19_24 0.138 2001.00 0.00087 24.74
988430 rs12295434 11_29 0.056 4584.19 0.00082 2.93
988439 rs11039677 11_29 0.056 4584.17 0.00081 2.93
988344 rs71045559 11_29 0.052 4583.64 0.00076 2.96
753462 rs113408695 17_39 1.000 204.06 0.00065 -17.59
984004 rs201524046 10_81 1.000 201.07 0.00064 -2.97
988422 rs7124318 11_29 0.044 4585.36 0.00064 2.93
57131 rs2211176 1_118 0.062 3204.98 0.00063 3.10SNPs 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
1081632 rs36013629 19_34 0.947 35230.23 1.1e-01 -26.15
1081607 rs10419198 19_34 0.053 35845.69 6.1e-03 -26.01
1075737 rs45512696 19_24 0.795 2003.90 5.1e-03 24.77
1075732 rs149131600 19_24 0.138 2001.00 8.7e-04 24.74
1075739 rs58895965 19_24 0.062 1998.47 3.9e-04 24.70
1075756 rs28616221 19_24 0.005 1949.95 3.3e-05 24.58
1075782 rs11671010 19_24 0.000 1928.68 7.3e-08 24.40
1075785 rs2018519 19_24 0.000 1930.46 1.3e-07 24.40
1081690 rs111476047 19_34 0.000 33389.94 0.0e+00 -23.87
1075698 rs41523449 19_24 1.000 1748.62 5.5e-03 23.52
1081564 rs10469298 19_34 0.000 56200.98 0.0e+00 -22.88
1081577 rs1132990 19_34 0.000 56295.28 0.0e+00 -22.80
72803 rs780093 2_16 1.000 418.22 1.3e-03 -22.27
1081540 rs2335534 19_34 0.000 92919.17 0.0e+00 -22.14
1075768 rs1688031 19_24 0.000 232.28 0.0e+00 22.11
1075766 rs2305746 19_24 0.000 467.92 0.0e+00 22.10
1075769 rs1672991 19_24 0.000 467.56 0.0e+00 22.10
1075724 rs1672981 19_24 0.000 471.08 0.0e+00 22.09
1075746 rs1688042 19_24 0.000 459.85 0.0e+00 22.07
1075749 rs1350292 19_24 0.000 459.68 0.0e+00 22.07
1075750 rs1350291 19_24 0.000 459.68 0.0e+00 22.07
1075747 rs1672979 19_24 0.000 459.31 0.0e+00 22.06
1075752 rs1350290 19_24 0.000 455.68 0.0e+00 22.05
1075744 rs1688043 19_24 0.000 466.08 0.0e+00 22.04
1075754 rs2452000 19_24 0.000 455.78 0.0e+00 22.04
1075755 rs2445819 19_24 0.000 455.82 0.0e+00 22.04
1075759 rs4806073 19_24 0.000 455.80 0.0e+00 22.04
1075743 rs55726903 19_24 0.000 453.57 0.0e+00 22.03
1075751 rs1688040 19_24 0.000 454.68 0.0e+00 22.03
1075758 rs1615767 19_24 0.000 455.01 0.0e+00 22.03
1075770 rs1672992 19_24 0.000 228.85 0.0e+00 22.03
1075741 rs1688044 19_24 0.000 453.19 0.0e+00 21.98
1075745 rs2445818 19_24 0.000 453.62 0.0e+00 21.98
1075771 rs1688030 19_24 0.000 458.01 0.0e+00 21.93
1075720 rs55714647 19_24 0.000 448.90 0.0e+00 21.81
1075784 rs897764 19_24 0.000 437.98 0.0e+00 21.51
1075800 rs10419703 19_24 0.000 416.43 0.0e+00 21.14
1081758 rs2379087 19_34 0.000 26871.15 0.0e+00 -20.86
1081814 rs10414643 19_34 0.000 27809.95 0.0e+00 -20.85
1081785 rs10426059 19_34 0.000 25545.88 0.0e+00 -20.82
1081798 rs7249925 19_34 0.000 26461.72 0.0e+00 -20.82
1081788 rs112727702 19_34 0.000 25927.47 0.0e+00 -20.80
1081762 rs10416310 19_34 0.000 26710.89 0.0e+00 -20.79
1081793 rs2288921 19_34 0.000 25966.37 0.0e+00 -20.79
1081807 rs11878568 19_34 0.000 26484.92 0.0e+00 -20.79
1081829 rs2116922 19_34 0.000 28474.55 0.0e+00 -20.79
1081782 rs8113357 19_34 0.000 25957.85 0.0e+00 -20.78
1081827 rs10417980 19_34 0.000 28346.97 0.0e+00 -20.78
1081779 rs2890072 19_34 0.000 26064.32 0.0e+00 -20.76
1081790 rs10406941 19_34 0.000 25992.98 0.0e+00 -20.76Gene 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] 32if (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"
Term
1 negative regulation of plasminogen activation (GO:0010757)
2 negative regulation of fibrinolysis (GO:0051918)
3 positive regulation of transforming growth factor beta production (GO:0071636)
4 regulation of plasminogen activation (GO:0010755)
5 regulation of fibrinolysis (GO:0051917)
6 negative regulation of protein processing (GO:0010955)
7 positive regulation of blood coagulation (GO:0030194)
8 response to progesterone (GO:0032570)
9 response to ketone (GO:1901654)
Overlap Adjusted.P.value Genes
1 2/5 0.01216499 SERPINF2;THBS1
2 2/9 0.01882651 SERPINF2;THBS1
3 2/11 0.01882651 SERPINF2;THBS1
4 2/12 0.01882651 SERPINF2;THBS1
5 2/13 0.01882651 SERPINF2;THBS1
6 2/15 0.02107735 SERPINF2;THBS1
7 2/17 0.02335348 SERPINF2;THBS1
8 2/21 0.03142720 CATSPER2;THBS1
9 2/28 0.04993343 CATSPER2;THBS1
[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 CURATEDLAMP1 gene(s) from the input list not found in DisGeNET CURATEDPELO gene(s) from the input list not found in DisGeNET CURATEDFCHO2 gene(s) from the input list not found in DisGeNET CURATEDFCGRT gene(s) from the input list not found in DisGeNET CURATEDSH3BP1 gene(s) from the input list not found in DisGeNET CURATEDTLCD2 gene(s) from the input list not found in DisGeNET CURATEDBEND3 gene(s) from the input list not found in DisGeNET CURATEDANGEL2 gene(s) from the input list not found in DisGeNET CURATEDRAVER2 gene(s) from the input list not found in DisGeNET CURATEDRAB17 gene(s) from the input list not found in DisGeNET CURATED Description
86 gliosarcoma
100 Giant Cell Glioblastoma
152 Deafness, Autosomal Recessive 28
157 Tyrosine Kinase 2 Deficiency
163 ALPHA-2-PLASMIN INHIBITOR DEFICIENCY
164 Gastro-enteropancreatic neuroendocrine tumor
165 NOONAN SYNDROME-LIKE DISORDER WITH OR WITHOUT JUVENILE MYELOMONOCYTIC LEUKEMIA
169 Anti-plasmin deficiency, congenital
174 ENCEPHALOPATHY, ACUTE, INFECTION-INDUCED (HERPES-SPECIFIC), SUSCEPTIBILITY TO, 5
177 DESBUQUOIS DYSPLASIA 2
FDR Ratio BgRatio
86 0.02891156 2/21 16/9703
100 0.02891156 3/21 84/9703
152 0.02891156 1/21 1/9703
157 0.02891156 1/21 1/9703
163 0.02891156 1/21 1/9703
164 0.02891156 2/21 15/9703
165 0.02891156 1/21 1/9703
169 0.02891156 1/21 1/9703
174 0.02891156 1/21 1/9703
177 0.02891156 1/21 1/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