Last updated: 2021-09-09
Checks: 6 1
Knit directory: ctwas_applied/
This reproducible R Markdown analysis was created with workflowr (version 1.6.2). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.
The R Markdown file has unstaged changes. To know which version of the R Markdown file created these results, you’ll want to first commit it to the Git repo. If you’re still working on the analysis, you can ignore this warning. When you’re finished, you can run wflow_publish
to commit the R Markdown file and build the HTML.
Great job! The global environment was empty. Objects defined in the global environment can affect the analysis in your R Markdown file in unknown ways. For reproduciblity it’s best to always run the code in an empty environment.
The command set.seed(20210726)
was run prior to running the code in the R Markdown file. Setting a seed ensures that any results that rely on randomness, e.g. subsampling or permutations, are reproducible.
Great job! Recording the operating system, R version, and package versions is critical for reproducibility.
Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.
Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.
Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility.
The results in this page were generated with repository version 59e5f4d. See the Past versions tab to see a history of the changes made to the R Markdown and HTML files.
Note that you need to be careful to ensure that all relevant files for the analysis have been committed to Git prior to generating the results (you can use wflow_publish
or wflow_git_commit
). workflowr only checks the R Markdown file, but you know if there are other scripts or data files that it depends on. Below is the status of the Git repository when the results were generated:
Unstaged changes:
Modified: analysis/ukb-d-30500_irnt_Liver.Rmd
Modified: analysis/ukb-d-30500_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30600_irnt_Liver.Rmd
Modified: analysis/ukb-d-30600_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30610_irnt_Liver.Rmd
Modified: analysis/ukb-d-30610_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30620_irnt_Liver.Rmd
Modified: analysis/ukb-d-30620_irnt_Whole_Blood.Rmd
Modified: analysis/ukb-d-30630_irnt_Liver.Rmd
Modified: analysis/ukb-d-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-30620_irnt_Whole_Blood.Rmd
) and HTML (docs/ukb-d-30620_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 |
These are the results of a ctwas
analysis of the UK Biobank trait Alanine aminotransferase (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-30620_irnt
. Results were obtained from from IEU rather than Neale Lab because they are in a standardard format (GWAS VCF). Note that 3 of the 34 biomarker traits were not available from IEU and were excluded from analysis.
The weights are mashr GTEx v8 models on Whole_Blood
eQTL obtained from PredictDB. We performed a full harmonization of the variants, including recovering strand ambiguous variants. This procedure is discussed in a separate document. (TO-DO: add report that describes harmonization)
LD matrices were computed from a 10% subset of Neale lab subjects. Subjects were matched using the plate and well information from genotyping. We included only biallelic variants with MAF>0.01 in the original Neale Lab GWAS. (TO-DO: add more details [number of subjects, variants, etc])
TO-DO: add enhanced QC reporting (total number of weights, why each variant was missing for all genes)
qclist_all <- list()
qc_files <- paste0(results_dir, "/", list.files(results_dir, pattern="exprqc.Rd"))
for (i in 1:length(qc_files)){
load(qc_files[i])
chr <- unlist(strsplit(rev(unlist(strsplit(qc_files[i], "_")))[1], "[.]"))[1]
qclist_all[[chr]] <- cbind(do.call(rbind, lapply(qclist,unlist)), as.numeric(substring(chr,4)))
}
qclist_all <- data.frame(do.call(rbind, qclist_all))
colnames(qclist_all)[ncol(qclist_all)] <- "chr"
rm(qclist, wgtlist, z_gene_chr)
#number of imputed weights
nrow(qclist_all)
[1] 11095
#number of imputed weights by chromosome
table(qclist_all$chr)
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15
1129 747 624 400 479 621 560 383 404 430 682 652 192 362 331
16 17 18 19 20 21 22
551 725 159 911 313 130 310
#proportion of imputed weights without missing variants
mean(qclist_all$nmiss==0)
[1] 0.762776
#load ctwas results
ctwas_res <- data.table::fread(paste0(results_dir, "/", analysis_id, "_ctwas.susieIrss.txt"))
#make unique identifier for regions
ctwas_res$region_tag <- paste(ctwas_res$region_tag1, ctwas_res$region_tag2, sep="_")
#compute PVE for each gene/SNP
ctwas_res$PVE = ctwas_res$susie_pip*ctwas_res$mu2/sample_size #check PVE calculation
#separate gene and SNP results
ctwas_gene_res <- ctwas_res[ctwas_res$type == "gene", ]
ctwas_gene_res <- data.frame(ctwas_gene_res)
ctwas_snp_res <- ctwas_res[ctwas_res$type == "SNP", ]
ctwas_snp_res <- data.frame(ctwas_snp_res)
#add gene information to results
sqlite <- RSQLite::dbDriver("SQLite")
db = RSQLite::dbConnect(sqlite, paste0("/project2/compbio/predictdb/mashr_models/mashr_", weight, ".db"))
query <- function(...) RSQLite::dbGetQuery(db, ...)
gene_info <- query("select gene, genename from extra")
gene_info <- query("select gene, genename, gene_type from extra")
RSQLite::dbDisconnect(db)
ctwas_gene_res <- cbind(ctwas_gene_res, gene_info[sapply(ctwas_gene_res$id, match, gene_info$gene), c("genename", "gene_type")])
#add z score to results
load(paste0(results_dir, "/", analysis_id, "_expr_z_gene.Rd"))
ctwas_gene_res$z <- z_gene[ctwas_gene_res$id,]$z
#load(paste0(results_dir, "/", analysis_id, "_expr_z_snp.Rd")) #for new version, stored after harmonization
z_snp <- readRDS(paste0(results_dir, "/", trait_id, ".RDS")) #for old version, unharmonized
z_snp <- z_snp[z_snp$id %in% ctwas_snp_res$id,] #subset snps to those included in analysis, note some are duplicated, need to match which allele was used
ctwas_snp_res$z <- z_snp$z[match(ctwas_snp_res$id, z_snp$id)] #for duplicated snps, this arbitrarily uses the first allele
ctwas_snp_res$z_flag <- as.numeric(ctwas_snp_res$id %in% z_snp$id[duplicated(z_snp$id)]) #mark the unclear z scores, flag=1
#formatting and rounding for tables
ctwas_gene_res$z <- round(ctwas_gene_res$z,2)
ctwas_snp_res$z <- round(ctwas_snp_res$z,2)
ctwas_gene_res$susie_pip <- round(ctwas_gene_res$susie_pip,3)
ctwas_snp_res$susie_pip <- round(ctwas_snp_res$susie_pip,3)
ctwas_gene_res$mu2 <- round(ctwas_gene_res$mu2,2)
ctwas_snp_res$mu2 <- round(ctwas_snp_res$mu2,2)
ctwas_gene_res$PVE <- signif(ctwas_gene_res$PVE, 2)
ctwas_snp_res$PVE <- signif(ctwas_snp_res$PVE, 2)
#merge gene and snp results with added information
ctwas_gene_res$z_flag=NA
ctwas_snp_res$genename=NA
ctwas_snp_res$gene_type=NA
ctwas_res <- rbind(ctwas_gene_res,
ctwas_snp_res[,colnames(ctwas_gene_res)])
#store columns to report
report_cols <- colnames(ctwas_gene_res)[!(colnames(ctwas_gene_res) %in% c("type", "region_tag1", "region_tag2", "cs_index", "gene_type", "z_flag", "id", "chrom", "pos"))]
first_cols <- c("genename", "region_tag")
report_cols <- c(first_cols, report_cols[!(report_cols %in% first_cols)])
report_cols_snps <- c("id", report_cols[-1])
#get number of SNPs from s1 results; adjust for thin
ctwas_res_s1 <- data.table::fread(paste0(results_dir, "/", analysis_id, "_ctwas.s1.susieIrss.txt"))
n_snps <- sum(ctwas_res_s1$type=="SNP")/thin
rm(ctwas_res_s1)
library(ggplot2)
library(cowplot)
********************************************************
Note: As of version 1.0.0, cowplot does not change the
default ggplot2 theme anymore. To recover the previous
behavior, execute:
theme_set(theme_cowplot())
********************************************************
load(paste0(results_dir, "/", analysis_id, "_ctwas.s2.susieIrssres.Rd"))
df <- data.frame(niter = rep(1:ncol(group_prior_rec), 2),
value = c(group_prior_rec[1,], group_prior_rec[2,]),
group = rep(c("Gene", "SNP"), each = ncol(group_prior_rec)))
df$group <- as.factor(df$group)
df$value[df$group=="SNP"] <- df$value[df$group=="SNP"]*thin #adjust parameter to account for thin argument
p_pi <- ggplot(df, aes(x=niter, y=value, group=group)) +
geom_line(aes(color=group)) +
geom_point(aes(color=group)) +
xlab("Iteration") + ylab(bquote(pi)) +
ggtitle("Prior mean") +
theme_cowplot()
df <- data.frame(niter = rep(1:ncol(group_prior_var_rec), 2),
value = c(group_prior_var_rec[1,], group_prior_var_rec[2,]),
group = rep(c("Gene", "SNP"), each = ncol(group_prior_var_rec)))
df$group <- as.factor(df$group)
p_sigma2 <- ggplot(df, aes(x=niter, y=value, group=group)) +
geom_line(aes(color=group)) +
geom_point(aes(color=group)) +
xlab("Iteration") + ylab(bquote(sigma^2)) +
ggtitle("Prior variance") +
theme_cowplot()
plot_grid(p_pi, p_sigma2)
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
#estimated group prior
estimated_group_prior <- group_prior_rec[,ncol(group_prior_rec)]
names(estimated_group_prior) <- c("gene", "snp")
estimated_group_prior["snp"] <- estimated_group_prior["snp"]*thin #adjust parameter to account for thin argument
print(estimated_group_prior)
gene snp
0.0186343157 0.0001879955
#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
11.96031 12.02071
#report sample size
print(sample_size)
[1] 344136
#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.007185434 0.057112804
#compare sum(PIP*mu2/sample_size) with above PVE calculation
c(sum(ctwas_gene_res$PVE),sum(ctwas_snp_res$PVE))
[1] 0.02928884 0.74343215
#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
9507 FES 15_43 0.997 31.75 9.2e-05 -5.57
2378 CAMTA2 17_5 0.996 235.82 6.8e-04 5.31
10004 FPR3 19_35 0.988 23.83 6.8e-05 -4.77
7955 MFSD3 8_94 0.985 53.08 1.5e-04 5.22
10841 MRPS18B 6_26 0.981 56.91 1.6e-04 -7.62
3058 PLEKHA3 2_108 0.975 34.25 9.7e-05 -5.69
5540 SYTL1 1_19 0.962 29.73 8.3e-05 -5.41
5690 WDPCP 2_41 0.960 26.90 7.5e-05 5.36
1778 KPNA3 13_21 0.955 21.70 6.0e-05 -4.39
7835 EVA1C 21_13 0.947 31.67 8.7e-05 5.55
5420 C18orf8 18_12 0.944 32.79 9.0e-05 -5.94
2540 NECTIN1 11_72 0.943 19.99 5.5e-05 -4.15
11104 STARD10 11_41 0.940 67.49 1.8e-04 9.15
5235 SUOX 12_35 0.934 49.65 1.3e-04 6.98
4360 TRIM5 11_4 0.928 39.30 1.1e-04 -5.20
5634 ADAM15 1_77 0.927 25.16 6.8e-05 -2.81
6533 WHAMM 15_38 0.927 27.27 7.3e-05 -4.47
325 BAK1 6_28 0.904 25.59 6.7e-05 4.90
8823 SH3PXD2B 5_103 0.886 22.44 5.8e-05 4.67
9223 ZBTB7A 19_4 0.877 24.95 6.4e-05 -4.84
11143 TMEM167B 1_67 0.868 21.37 5.4e-05 4.45
3758 ATXN1 6_13 0.868 20.25 5.1e-05 4.15
10313 ZFP62 5_109 0.853 20.36 5.0e-05 4.15
11411 HIST1H2BN 6_21 0.851 18.69 4.6e-05 3.77
3833 GPCPD1 20_5 0.846 22.00 5.4e-05 4.22
8124 DDX19A 16_37 0.844 20.03 4.9e-05 -4.22
8829 HRAS 11_1 0.829 23.21 5.6e-05 4.52
3101 MEF2D 1_77 0.828 46.09 1.1e-04 6.91
5100 OIT3 10_48 0.827 27.66 6.6e-05 5.56
208 PPP5C 19_32 0.823 21.51 5.1e-05 4.29
209 CEP68 2_42 0.813 29.21 6.9e-05 5.38
487 FOXC1 6_2 0.813 21.82 5.2e-05 -4.35
#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
4687 TMEM60 7_49 0.000 1540.88 0.0e+00 3.75
2528 PANX1 11_53 0.000 1461.44 0.0e+00 13.17
9608 PSMG1 21_19 0.000 1347.00 0.0e+00 -6.78
9834 BRWD1 21_19 0.000 636.42 0.0e+00 1.07
3657 GPR83 11_53 0.000 531.67 0.0e+00 -2.81
1366 CWF19L1 10_64 0.211 503.74 3.1e-04 -24.40
8411 TRMT61B 2_19 0.000 463.33 2.6e-16 2.53
11039 PPP1CB 2_19 0.000 455.93 2.5e-17 -2.33
10204 BLOC1S2 10_64 0.000 406.60 1.0e-08 -21.88
10931 HMGN1 21_19 0.000 397.10 0.0e+00 -1.59
11094 APTR 7_49 0.000 306.65 0.0e+00 1.10
9466 WRB 21_19 0.000 243.07 0.0e+00 -2.96
268 MRE11 11_53 0.000 240.72 0.0e+00 -0.49
2378 CAMTA2 17_5 0.996 235.82 6.8e-04 5.31
397 MED17 11_53 0.000 235.70 0.0e+00 2.20
4095 KIF1C 17_5 0.018 232.55 1.2e-05 4.18
10258 INCA1 17_5 0.030 221.19 2.0e-05 3.75
6626 LCA5L 21_19 0.000 205.24 0.0e+00 -2.86
5934 MFHAS1 8_11 0.000 204.76 7.0e-10 5.41
9806 SH3BGR 21_19 0.000 201.83 0.0e+00 1.20
#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
2378 CAMTA2 17_5 0.996 235.82 6.8e-04 5.31
1366 CWF19L1 10_64 0.211 503.74 3.1e-04 -24.40
11104 STARD10 11_41 0.940 67.49 1.8e-04 9.15
2849 SMC4 3_99 0.774 74.17 1.7e-04 8.76
10841 MRPS18B 6_26 0.981 56.91 1.6e-04 -7.62
7955 MFSD3 8_94 0.985 53.08 1.5e-04 5.22
5235 SUOX 12_35 0.934 49.65 1.3e-04 6.98
8813 MSL2 3_84 0.682 57.38 1.1e-04 10.26
3101 MEF2D 1_77 0.828 46.09 1.1e-04 6.91
4360 TRIM5 11_4 0.928 39.30 1.1e-04 -5.20
9656 KMT5A 12_75 0.733 47.04 1.0e-04 -8.04
3058 PLEKHA3 2_108 0.975 34.25 9.7e-05 -5.69
3378 GPAM 10_70 0.660 47.77 9.2e-05 7.40
8284 SPDYE5 7_48 0.720 44.07 9.2e-05 6.83
9507 FES 15_43 0.997 31.75 9.2e-05 -5.57
5420 C18orf8 18_12 0.944 32.79 9.0e-05 -5.94
7835 EVA1C 21_13 0.947 31.67 8.7e-05 5.55
10505 UGT2B17 4_48 0.787 36.95 8.5e-05 -6.01
5540 SYTL1 1_19 0.962 29.73 8.3e-05 -5.41
2607 SH2B3 12_67 0.626 44.03 8.0e-05 9.03
#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
1366 CWF19L1 10_64 0.211 503.74 3.1e-04 -24.40
10204 BLOC1S2 10_64 0.000 406.60 1.0e-08 -21.88
4137 MAU2 19_15 0.015 150.04 6.7e-06 -15.13
11446 RP11-9M16.2 9_59 0.006 201.16 3.4e-06 -14.56
2528 PANX1 11_53 0.000 1461.44 0.0e+00 13.17
6949 RPL8 8_94 0.000 137.02 2.2e-10 -12.54
6175 DLG5 10_50 0.032 69.77 6.5e-06 10.46
8813 MSL2 3_84 0.682 57.38 1.1e-04 10.26
10765 ZDHHC18 1_18 0.032 88.50 8.2e-06 9.92
2131 ATP13A1 19_15 0.002 78.73 5.6e-07 9.71
6551 SUPV3L1 10_46 0.022 76.39 4.9e-06 -9.19
11104 STARD10 11_41 0.940 67.49 1.8e-04 9.15
2289 DNMBP 10_64 0.000 110.27 4.3e-09 -9.08
10255 ZNF34 8_94 0.000 115.63 9.3e-11 9.04
2607 SH2B3 12_67 0.626 44.03 8.0e-05 9.03
2292 ERLIN1 10_64 0.000 89.62 2.8e-08 -8.84
2849 SMC4 3_99 0.774 74.17 1.7e-04 8.76
11726 YJEFN3 19_15 0.008 52.11 1.3e-06 8.76
4451 RSG1 1_11 0.014 66.90 2.8e-06 -8.56
9813 MUC1 1_77 0.002 80.65 3.8e-07 -8.55
#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.02045967
#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
1366 CWF19L1 10_64 0.211 503.74 3.1e-04 -24.40
10204 BLOC1S2 10_64 0.000 406.60 1.0e-08 -21.88
4137 MAU2 19_15 0.015 150.04 6.7e-06 -15.13
11446 RP11-9M16.2 9_59 0.006 201.16 3.4e-06 -14.56
2528 PANX1 11_53 0.000 1461.44 0.0e+00 13.17
6949 RPL8 8_94 0.000 137.02 2.2e-10 -12.54
6175 DLG5 10_50 0.032 69.77 6.5e-06 10.46
8813 MSL2 3_84 0.682 57.38 1.1e-04 10.26
10765 ZDHHC18 1_18 0.032 88.50 8.2e-06 9.92
2131 ATP13A1 19_15 0.002 78.73 5.6e-07 9.71
6551 SUPV3L1 10_46 0.022 76.39 4.9e-06 -9.19
11104 STARD10 11_41 0.940 67.49 1.8e-04 9.15
2289 DNMBP 10_64 0.000 110.27 4.3e-09 -9.08
10255 ZNF34 8_94 0.000 115.63 9.3e-11 9.04
2607 SH2B3 12_67 0.626 44.03 8.0e-05 9.03
2292 ERLIN1 10_64 0.000 89.62 2.8e-08 -8.84
2849 SMC4 3_99 0.774 74.17 1.7e-04 8.76
11726 YJEFN3 19_15 0.008 52.11 1.3e-06 8.76
4451 RSG1 1_11 0.014 66.90 2.8e-06 -8.56
9813 MUC1 1_77 0.002 80.65 3.8e-07 -8.55
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: 10_64"
genename region_tag susie_pip mu2 PVE z
3390 GOT1 10_64 0.000 21.02 4.3e-10 -4.11
11229 RP11-441O15.3 10_64 0.000 21.63 6.7e-10 -3.61
6475 SLC25A28 10_64 0.000 12.50 1.2e-10 -2.34
11988 RP11-85A1.3 10_64 0.000 12.33 1.2e-10 2.31
10532 ENTPD7 10_64 0.000 25.61 9.5e-10 3.63
3379 CUTC 10_64 0.000 17.85 1.7e-10 3.51
244 COX15 10_64 0.000 8.21 7.6e-11 0.89
290 ABCC2 10_64 0.000 37.61 1.3e-08 2.46
2289 DNMBP 10_64 0.000 110.27 4.3e-09 -9.08
2292 ERLIN1 10_64 0.000 89.62 2.8e-08 -8.84
1366 CWF19L1 10_64 0.211 503.74 3.1e-04 -24.40
10204 BLOC1S2 10_64 0.000 406.60 1.0e-08 -21.88
2294 PKD2L1 10_64 0.000 32.62 8.1e-10 4.96
11463 OLMALINC 10_64 0.000 17.20 7.3e-10 -2.20
891 SEC31B 10_64 0.000 9.30 1.6e-10 -0.86
7681 HIF1AN 10_64 0.000 7.67 1.1e-10 -0.38
7682 NDUFB8 10_64 0.000 5.67 6.0e-11 -0.10
3375 SLF2 10_64 0.796 29.49 6.8e-05 5.10
1367 SEMA4G 10_64 0.082 25.08 6.0e-06 4.01
2313 TWNK 10_64 0.011 20.86 6.8e-07 4.64
504 MRPL43 10_64 0.004 18.44 2.3e-07 -3.42
2314 LZTS2 10_64 0.033 22.86 2.2e-06 3.82
9967 PDZD7 10_64 0.000 4.97 5.1e-11 0.90
2315 SFXN3 10_64 0.000 6.99 8.8e-11 0.95
2316 KAZALD1 10_64 0.000 7.83 9.5e-11 0.19
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 19_15"
genename region_tag susie_pip mu2 PVE z
4199 LSM4 19_15 0.004 11.76 1.4e-07 -1.65
4197 PGPEP1 19_15 0.004 11.39 1.4e-07 -1.27
8907 LRRC25 19_15 0.010 22.03 6.4e-07 3.12
4196 SSBP4 19_15 0.002 5.07 3.1e-08 0.93
2112 ISYNA1 19_15 0.002 5.19 3.2e-08 -0.81
2113 ELL 19_15 0.014 24.63 9.8e-07 -3.21
2123 KXD1 19_15 0.002 9.60 6.5e-08 1.69
11192 UBA52 19_15 0.003 9.42 8.8e-08 1.54
7904 KLHL26 19_15 0.002 10.62 6.9e-08 2.24
52 UPF1 19_15 0.003 7.68 7.2e-08 -0.15
2115 COPE 19_15 0.002 5.10 3.2e-08 -0.40
2116 DDX49 19_15 0.003 8.91 8.7e-08 1.09
2118 ARMC6 19_15 0.002 5.24 3.7e-08 0.53
599 SUGP2 19_15 0.002 5.07 3.2e-08 -0.17
596 TMEM161A 19_15 0.004 12.39 1.4e-07 1.46
11075 MEF2B 19_15 0.064 53.86 1.0e-05 -6.43
11817 BORCS8 19_15 0.002 15.11 9.6e-08 -5.47
595 RFXANK 19_15 0.002 5.69 3.7e-08 0.57
4137 MAU2 19_15 0.015 150.04 6.7e-06 -15.13
7905 GATAD2A 19_15 0.008 45.36 1.0e-06 7.00
9879 NDUFA13 19_15 0.005 41.67 6.2e-07 6.71
9152 TSSK6 19_15 0.008 15.09 3.6e-07 -4.98
11726 YJEFN3 19_15 0.008 52.11 1.3e-06 8.76
6840 CILP2 19_15 0.008 14.84 3.4e-07 4.96
2128 PBX4 19_15 0.003 6.46 5.3e-08 2.93
597 LPAR2 19_15 0.037 27.91 3.0e-06 5.92
1235 GMIP 19_15 0.049 26.95 3.9e-06 5.34
2131 ATP13A1 19_15 0.002 78.73 5.6e-07 9.71
9450 ZNF101 19_15 0.004 9.29 1.0e-07 4.14
2126 ZNF14 19_15 0.002 8.48 5.5e-08 -0.67
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 9_59"
genename region_tag susie_pip mu2 PVE z
11333 ORM1 9_59 0.008 9.73 2.2e-07 2.67
11310 ORM2 9_59 0.013 13.15 4.9e-07 0.94
11446 RP11-9M16.2 9_59 0.006 201.16 3.4e-06 -14.56
4914 ATP6V1G1 9_59 0.006 6.84 1.3e-07 0.61
6633 TMEM268 9_59 0.010 20.88 5.9e-07 3.91
2258 TNFSF8 9_59 0.005 5.89 9.2e-08 1.59
9428 TNFSF15 9_59 0.008 8.54 1.9e-07 -0.80
392 TNC 9_59 0.007 7.82 1.6e-07 -1.15
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 11_53"
genename region_tag susie_pip mu2 PVE z
7666 CEP295 11_53 0 62.46 0 0.30
397 MED17 11_53 0 235.70 0 2.20
2528 PANX1 11_53 0 1461.44 0 13.17
3657 GPR83 11_53 0 531.67 0 -2.81
268 MRE11 11_53 0 240.72 0 -0.49
8126 ANKRD49 11_53 0 126.76 0 -3.89
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
[1] "Region: 8_94"
genename region_tag susie_pip mu2 PVE z
9429 ZFP41 8_94 0.000 7.53 3.3e-15 -0.94
11739 GLI4 8_94 0.000 7.09 2.9e-15 0.56
9844 ZNF696 8_94 0.000 8.16 3.8e-15 -0.98
9696 TOP1MT 8_94 0.000 41.62 5.7e-13 3.16
6670 RHPN1 8_94 0.000 9.70 5.2e-15 -1.13
1943 GSDMD 8_94 0.000 24.94 8.2e-14 1.88
239 ZC3H3 8_94 0.000 7.96 4.3e-15 -1.44
10868 MROH6 8_94 0.000 6.00 2.4e-15 0.98
5965 NAPRT 8_94 0.000 62.41 9.5e-12 -4.12
1944 EEF1D 8_94 0.000 62.95 1.0e-11 -3.46
9289 TIGD5 8_94 0.000 9.70 6.8e-15 -0.56
9398 ZNF707 8_94 0.000 14.16 1.1e-14 2.04
9378 FAM83H 8_94 0.000 18.58 3.2e-14 1.62
9302 PUF60 8_94 0.000 5.70 2.0e-15 -0.37
9164 PLEC 8_94 0.000 9.63 1.0e-14 -0.64
9194 PARP10 8_94 0.000 10.27 1.1e-14 -1.17
9198 GRINA 8_94 0.000 24.95 2.1e-13 -1.99
9209 OPLAH 8_94 0.000 16.38 9.4e-15 2.43
9240 CYC1 8_94 0.000 10.40 4.8e-15 1.08
9276 MAF1 8_94 0.000 6.89 2.5e-15 0.53
9278 WDR97 8_94 0.000 6.43 3.1e-15 -0.38
9268 SHARPIN 8_94 0.000 6.74 2.4e-15 -0.50
9766 HSF1 8_94 0.000 7.75 2.9e-15 -1.94
11985 SCX 8_94 0.000 6.28 2.4e-15 -0.88
9752 DGAT1 8_94 0.000 14.01 7.1e-15 1.22
12034 TMEM249 8_94 0.000 9.09 3.5e-15 -0.82
8663 ADCK5 8_94 0.000 9.98 4.0e-15 0.60
785 CPSF1 8_94 0.000 19.42 3.8e-14 0.30
6938 TONSL 8_94 0.000 20.96 3.9e-14 1.92
7957 KIFC2 8_94 0.000 33.78 6.8e-14 -4.79
6943 PPP1R16A 8_94 0.000 112.03 1.4e-13 6.05
6940 RECQL4 8_94 0.004 178.54 2.2e-06 7.95
7956 GPT 8_94 0.000 19.60 1.2e-14 4.98
7955 MFSD3 8_94 0.985 53.08 1.5e-04 5.22
6941 LRRC14 8_94 0.000 135.11 1.3e-13 7.72
5964 SLC39A4 8_94 0.000 16.37 2.8e-14 0.16
10559 ZNF251 8_94 0.000 16.16 6.6e-15 2.50
10255 ZNF34 8_94 0.000 115.63 9.3e-11 9.04
6949 RPL8 8_94 0.000 137.02 2.2e-10 -12.54
10422 ZNF517 8_94 0.000 57.30 1.2e-12 -6.87
8358 COMMD5 8_94 0.000 85.75 2.9e-10 -5.93
10214 ZNF250 8_94 0.000 36.82 6.4e-13 1.76
8360 ZNF16 8_94 0.000 28.12 6.8e-14 2.47
9490 C8orf33 8_94 0.000 47.38 3.2e-12 2.03
Version | Author | Date |
---|---|---|
dfd2b5f | wesleycrouse | 2021-09-07 |
#snps with PIP>0.8 or 20 highest PIPs
head(ctwas_snp_res[order(-ctwas_snp_res$susie_pip),report_cols_snps],
max(sum(ctwas_snp_res$susie_pip>0.8), 20))
id region_tag susie_pip mu2 PVE z
4517 rs4336844 1_11 1.000 182.38 5.3e-04 14.42
7755 rs79598313 1_18 1.000 111.95 3.3e-04 11.19
53029 rs2642420 1_112 1.000 43.08 1.3e-04 -6.69
62118 rs12239046 1_131 1.000 48.67 1.4e-04 -7.26
72145 rs569546056 2_19 1.000 832.92 2.4e-03 3.68
91172 rs894194 2_55 1.000 39.75 1.2e-04 -6.40
112383 rs1862069 2_102 1.000 50.73 1.5e-04 8.48
147330 rs2649750 3_28 1.000 34.60 1.0e-04 -6.03
171928 rs9870956 3_77 1.000 49.99 1.5e-04 7.14
181653 rs9817452 3_97 1.000 45.72 1.3e-04 6.81
186146 rs149368105 3_105 1.000 105.46 3.1e-04 -11.00
241496 rs11727676 4_94 1.000 43.34 1.3e-04 6.79
275686 rs536916238 5_33 1.000 42.84 1.2e-04 -0.43
290229 rs163895 5_63 1.000 35.05 1.0e-04 -5.66
392506 rs10277379 7_49 1.000 2675.85 7.8e-03 -4.68
392509 rs761767938 7_49 1.000 3442.53 1.0e-02 -3.42
392517 rs1544459 7_49 1.000 3379.03 9.8e-03 -3.87
422115 rs2428 8_11 1.000 592.65 1.7e-03 8.64
422120 rs758184196 8_11 1.000 642.54 1.9e-03 -2.30
428749 rs2293400 8_23 1.000 46.62 1.4e-04 5.84
437165 rs140753685 8_42 1.000 42.24 1.2e-04 6.62
493797 rs7040440 9_59 1.000 62.06 1.8e-04 -1.99
493805 rs2763193 9_59 1.000 279.81 8.1e-04 -17.31
493806 rs10739409 9_59 1.000 117.29 3.4e-04 -14.77
523481 rs9645500 10_46 1.000 149.90 4.4e-04 12.94
532240 rs139450722 10_64 1.000 92.15 2.7e-04 -4.95
562420 rs17157266 11_34 1.000 47.42 1.4e-04 -7.06
571107 rs144988974 11_52 1.000 37.01 1.1e-04 6.20
571291 rs74717621 11_54 1.000 37.52 1.1e-04 6.45
577665 rs1176746 11_67 1.000 1141.65 3.3e-03 -2.94
577667 rs2307599 11_67 1.000 1141.79 3.3e-03 -2.75
588457 rs12824533 12_11 1.000 64.00 1.9e-04 8.23
591300 rs66720652 12_15 1.000 49.16 1.4e-04 -6.83
659222 rs2332328 14_3 1.000 42.81 1.2e-04 6.77
667696 rs72681869 14_20 1.000 54.77 1.6e-04 -9.53
681621 rs1243165 14_49 1.000 50.17 1.5e-04 3.73
691494 rs511338 15_14 1.000 44.36 1.3e-04 7.09
696877 rs2070895 15_27 1.000 52.37 1.5e-04 -7.39
729354 rs11645522 16_46 1.000 63.73 1.9e-04 7.57
731241 rs2255451 16_49 1.000 50.06 1.5e-04 -7.19
750513 rs1801689 17_38 1.000 125.15 3.6e-04 11.78
752292 rs1477066 17_41 1.000 46.72 1.4e-04 6.87
755830 rs62076019 17_46 1.000 97.54 2.8e-04 10.38
772181 rs62094231 18_31 1.000 68.66 2.0e-04 -8.33
793518 rs814573 19_31 1.000 94.79 2.8e-04 -11.60
825489 rs62221472 21_23 1.000 48.03 1.4e-04 -6.44
834737 rs139050 22_19 1.000 158.48 4.6e-04 -12.57
834738 rs6006585 22_19 1.000 63.28 1.8e-04 -5.48
844023 rs4989532 1_6 1.000 498.41 1.4e-03 -3.70
844025 rs115843159 1_6 1.000 63.86 1.9e-04 0.87
848095 rs35130213 1_19 1.000 461.64 1.3e-03 2.49
855500 rs140584594 1_67 1.000 74.80 2.2e-04 8.43
951448 rs2844543 6_26 1.000 133.09 3.9e-04 10.88
951785 rs112436252 6_26 1.000 72.23 2.1e-04 6.91
987691 rs873884 8_94 1.000 213.09 6.2e-04 11.89
1011559 rs11601507 11_4 1.000 96.61 2.8e-04 9.71
1015846 rs2511241 11_41 1.000 39.24 1.1e-04 -6.84
1020885 rs148050219 11_53 1.000 37770.02 1.1e-01 -15.40
1020895 rs111443113 11_53 1.000 37728.14 1.1e-01 -0.61
1066876 rs373230966 17_5 1.000 412.37 1.2e-03 3.07
1135585 rs34079499 21_19 1.000 4321.11 1.3e-02 -5.26
77909 rs72800939 2_28 0.999 32.41 9.4e-05 5.73
532073 rs111286300 10_64 0.999 42.44 1.2e-04 8.46
709547 rs9788910 16_3 0.999 34.10 9.9e-05 7.01
834739 rs11090617 22_19 0.999 594.73 1.7e-03 30.16
1075621 rs117643180 17_6 0.999 79.28 2.3e-04 9.29
53007 rs337171 1_112 0.998 48.80 1.4e-04 7.55
783550 rs8113613 19_9 0.998 71.69 2.1e-04 -6.23
792417 rs11879413 19_30 0.998 32.95 9.6e-05 -5.23
301099 rs769204262 5_84 0.997 30.33 8.8e-05 5.47
321252 rs115740542 6_20 0.997 31.70 9.2e-05 5.67
392105 rs12539160 7_47 0.997 34.31 9.9e-05 5.46
561872 rs1593480 11_34 0.997 32.33 9.4e-05 -5.53
783547 rs10854115 19_9 0.997 44.04 1.3e-04 -2.96
223425 rs77094191 4_59 0.996 60.26 1.7e-04 -4.92
403677 rs10435378 7_70 0.996 46.37 1.3e-04 6.83
550446 rs7481951 11_15 0.996 48.04 1.4e-04 7.20
743238 rs3760511 17_22 0.996 28.29 8.2e-05 -5.16
832404 rs132642 22_14 0.996 164.71 4.8e-04 13.69
422095 rs11774455 8_11 0.995 113.11 3.3e-04 7.51
790090 rs889140 19_23 0.995 41.56 1.2e-04 -8.34
47962 rs4951163 1_104 0.994 32.28 9.3e-05 5.08
186167 rs234043 3_106 0.994 28.86 8.3e-05 -5.30
772474 rs12373325 18_31 0.993 167.58 4.8e-04 -14.96
322033 rs2235698 6_23 0.992 31.36 9.0e-05 5.96
457437 rs10956254 8_83 0.992 31.52 9.1e-05 7.08
987654 rs112574791 8_94 0.992 465.69 1.3e-03 -21.20
38163 rs9425587 1_84 0.991 29.94 8.6e-05 -5.39
274729 rs28499105 5_31 0.990 31.48 9.1e-05 5.56
790082 rs12610925 19_23 0.990 49.94 1.4e-04 -8.65
807380 rs6129631 20_24 0.989 45.84 1.3e-04 -5.34
241704 rs4552481 4_95 0.988 239.74 6.9e-04 16.59
656101 rs7318424 13_59 0.987 34.87 1.0e-04 -5.90
809285 rs6066141 20_29 0.987 32.35 9.3e-05 -5.84
16463 rs12140153 1_39 0.986 27.56 7.9e-05 -5.46
422859 rs11777976 8_13 0.985 58.40 1.7e-04 -8.60
464149 rs1016565 9_1 0.984 25.82 7.4e-05 4.91
844024 rs2072735 1_6 0.983 495.86 1.4e-03 -4.12
673581 rs12894822 14_33 0.981 43.05 1.2e-04 6.62
298750 rs72791146 5_79 0.979 29.77 8.5e-05 5.35
457396 rs13252684 8_83 0.979 53.78 1.5e-04 3.57
667779 rs2883893 14_20 0.979 29.33 8.3e-05 -5.00
701147 rs12594627 15_35 0.978 63.87 1.8e-04 8.23
722872 rs72784008 16_31 0.978 38.31 1.1e-04 -6.24
179322 rs6774253 3_92 0.973 37.74 1.1e-04 -7.57
486605 rs113609637 9_47 0.971 31.61 8.9e-05 -5.79
786198 rs11668601 19_14 0.965 86.42 2.4e-04 9.97
591656 rs67981690 12_16 0.964 34.80 9.7e-05 5.54
729353 rs13334801 16_46 0.964 33.67 9.4e-05 4.78
30332 rs1730862 1_66 0.961 28.62 8.0e-05 -5.25
782395 rs10401485 19_7 0.959 27.01 7.5e-05 4.98
477234 rs11557154 9_27 0.958 34.95 9.7e-05 6.11
815725 rs2823025 21_2 0.958 29.03 8.1e-05 5.29
543580 rs7939634 11_2 0.957 27.25 7.6e-05 -5.03
802595 rs11557577 20_13 0.957 27.04 7.5e-05 -4.99
177104 rs12695769 3_87 0.948 44.46 1.2e-04 -6.83
808536 rs4810422 20_28 0.945 29.58 8.1e-05 5.40
275668 rs173964 5_33 0.943 106.88 2.9e-04 8.09
188268 rs13089089 3_110 0.935 24.26 6.6e-05 -4.47
580906 rs11220136 11_77 0.934 36.66 9.9e-05 6.02
732626 rs75079463 16_51 0.934 24.53 6.7e-05 4.61
1093989 rs58542926 19_15 0.927 272.55 7.3e-04 19.31
792431 rs28602288 19_30 0.925 28.05 7.5e-05 4.57
154307 rs62247577 3_43 0.921 23.65 6.3e-05 4.37
681617 rs941594 14_49 0.920 56.88 1.5e-04 4.99
303634 rs35341726 5_88 0.918 25.06 6.7e-05 4.76
409 rs10910028 1_2 0.917 45.42 1.2e-04 7.16
171948 rs141809192 3_78 0.914 24.22 6.4e-05 4.15
691486 rs17723097 15_13 0.914 42.80 1.1e-04 7.20
774155 rs73963711 18_35 0.913 25.00 6.6e-05 -5.02
422385 rs13265179 8_12 0.911 88.78 2.3e-04 10.81
347271 rs7752846 6_75 0.908 24.42 6.4e-05 -4.77
138463 rs115532219 3_9 0.907 25.76 6.8e-05 3.99
268706 rs13172112 5_21 0.899 46.22 1.2e-04 9.04
322624 rs1126511 6_27 0.892 40.87 1.1e-04 7.43
73473 rs72787520 2_20 0.889 23.94 6.2e-05 -4.53
565578 rs11236797 11_42 0.888 28.01 7.2e-05 -4.85
138281 rs13085211 3_9 0.884 152.89 3.9e-04 -10.80
654691 rs61965347 13_56 0.884 27.44 7.0e-05 5.11
571267 rs72980431 11_54 0.880 24.09 6.2e-05 -4.58
30969 rs325937 1_69 0.873 26.12 6.6e-05 -4.74
422136 rs13265731 8_11 0.873 580.93 1.5e-03 8.08
279872 rs3010273 5_43 0.870 35.64 9.0e-05 -6.07
987989 rs113284751 8_94 0.870 58.28 1.5e-04 -4.33
667744 rs142004400 14_20 0.866 48.06 1.2e-04 -8.89
457385 rs79658059 8_83 0.858 176.75 4.4e-04 -15.98
465759 rs6915 9_5 0.854 26.32 6.5e-05 -4.88
45973 rs12144388 1_99 0.852 27.90 6.9e-05 4.99
297743 rs34576922 5_78 0.851 32.23 8.0e-05 -5.67
684894 rs17617994 14_54 0.847 31.02 7.6e-05 -5.37
982512 rs1050969 8_69 0.847 29.22 7.2e-05 -4.08
805790 rs4911358 20_20 0.844 54.60 1.3e-04 6.24
298906 rs6894249 5_79 0.841 28.92 7.1e-05 -5.04
697511 rs340029 15_27 0.841 50.15 1.2e-04 7.15
195658 rs13116176 4_4 0.839 30.42 7.4e-05 -6.04
180156 rs7610095 3_94 0.837 36.37 8.8e-05 -6.13
98359 rs192728998 2_70 0.835 24.26 5.9e-05 -4.50
72144 rs7562170 2_19 0.833 801.87 1.9e-03 3.52
629940 rs1778790 13_7 0.833 57.79 1.4e-04 -7.85
179270 rs7645585 3_92 0.823 54.89 1.3e-04 -7.96
705143 rs35408448 15_41 0.823 71.72 1.7e-04 -8.67
10338 rs260970 1_24 0.822 25.00 6.0e-05 4.55
71938 rs937813 2_16 0.819 27.47 6.5e-05 5.15
679984 rs243215 14_45 0.819 26.29 6.3e-05 4.77
72148 rs4580350 2_19 0.811 803.08 1.9e-03 -3.46
769143 rs9953845 18_26 0.811 30.88 7.3e-05 5.33
#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
1020885 rs148050219 11_53 1.000 37770.02 1.1e-01 -15.40
1020895 rs111443113 11_53 1.000 37728.14 1.1e-01 -0.61
1020894 rs60550219 11_53 0.140 37725.38 1.5e-02 -15.42
1020881 rs7105405 11_53 0.157 37715.61 1.7e-02 -15.44
1020919 rs67167563 11_53 0.067 37712.32 7.3e-03 -15.41
1020926 rs113426210 11_53 0.531 37699.41 5.8e-02 -15.46
1020877 rs9888156 11_53 0.000 37657.48 5.0e-06 -15.40
1020932 rs950878 11_53 0.007 37650.52 7.8e-04 -15.44
1020875 rs67232024 11_53 0.000 37586.96 1.1e-08 -15.34
1020854 rs7927828 11_53 0.000 37583.40 9.1e-11 -15.29
1020872 rs9888266 11_53 0.000 37532.49 2.3e-11 -15.31
1020860 rs67812366 11_53 0.000 37531.84 1.6e-11 -15.31
1020863 rs7109132 11_53 0.000 37531.70 1.3e-11 -15.31
1020855 rs57856352 11_53 0.000 37519.25 2.0e-13 -15.27
1020876 rs16919533 11_53 0.000 37518.96 4.1e-12 -15.37
1020874 rs67549397 11_53 0.000 37463.49 4.3e-15 -15.29
1020873 rs9888143 11_53 0.000 37403.19 1.2e-17 -15.27
1020864 rs60351354 11_53 0.000 37399.02 0.0e+00 -15.27
1020865 rs60546087 11_53 0.000 37398.98 0.0e+00 -15.27
1020869 rs1573567 11_53 0.000 37398.95 0.0e+00 -15.27
1020866 rs7109819 11_53 0.000 37398.87 0.0e+00 -15.27
1020832 rs7932290 11_53 0.000 37218.99 0.0e+00 -15.28
1020800 rs7934467 11_53 0.000 36913.01 0.0e+00 -15.03
1021195 rs72966603 11_53 0.000 31145.65 0.0e+00 -16.38
1021325 rs12419615 11_53 0.000 29498.09 0.0e+00 -16.55
1021376 rs58964858 11_53 0.000 25272.82 0.0e+00 -16.59
1021378 rs72968738 11_53 0.000 25224.88 0.0e+00 -16.54
1021402 rs138626734 11_53 0.000 24921.79 0.0e+00 -16.62
1021420 rs4408267 11_53 0.000 24909.30 0.0e+00 -16.60
1021388 rs72968745 11_53 0.000 24896.33 0.0e+00 -16.47
1021387 rs4491178 11_53 0.000 24894.79 0.0e+00 -16.46
1021448 rs11604580 11_53 0.000 24881.93 0.0e+00 -16.62
1021453 rs4342991 11_53 0.000 24879.20 0.0e+00 -16.62
1021392 rs4753124 11_53 0.000 24819.88 0.0e+00 -16.63
1021425 rs16919942 11_53 0.000 24802.25 0.0e+00 -16.64
1021306 rs7945841 11_53 0.000 24742.95 0.0e+00 -14.94
1021019 rs72962880 11_53 0.000 24719.19 0.0e+00 -12.36
1021008 rs55659547 11_53 0.000 24647.19 0.0e+00 -12.33
1021007 rs7950356 11_53 0.000 24642.75 0.0e+00 -12.32
1021018 rs56359140 11_53 0.000 24616.13 0.0e+00 -12.29
1021011 rs72962872 11_53 0.000 24613.98 0.0e+00 -12.30
1021013 rs140989262 11_53 0.000 24291.75 0.0e+00 -12.26
1021344 rs7119800 11_53 0.000 24259.32 0.0e+00 -15.03
1021033 rs72962891 11_53 0.000 24095.40 0.0e+00 -12.23
1021052 rs72964604 11_53 0.000 24046.27 0.0e+00 -12.26
1021348 rs2176565 11_53 0.000 23950.34 0.0e+00 -15.17
1021349 rs7949551 11_53 0.000 23243.87 0.0e+00 -15.51
1020815 rs1506657 11_53 0.000 22836.55 0.0e+00 12.23
1021352 rs72968710 11_53 0.000 22722.53 0.0e+00 -15.74
1021355 rs16919917 11_53 0.000 22529.43 0.0e+00 -15.60
#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
1020885 rs148050219 11_53 1.000 37770.02 0.11000 -15.40
1020895 rs111443113 11_53 1.000 37728.14 0.11000 -0.61
1020926 rs113426210 11_53 0.531 37699.41 0.05800 -15.46
1020881 rs7105405 11_53 0.157 37715.61 0.01700 -15.44
1020894 rs60550219 11_53 0.140 37725.38 0.01500 -15.42
1135585 rs34079499 21_19 1.000 4321.11 0.01300 -5.26
392509 rs761767938 7_49 1.000 3442.53 0.01000 -3.42
392517 rs1544459 7_49 1.000 3379.03 0.00980 -3.87
392506 rs10277379 7_49 1.000 2675.85 0.00780 -4.68
1020919 rs67167563 11_53 0.067 37712.32 0.00730 -15.41
392513 rs11972122 7_49 0.778 3154.12 0.00710 -4.21
1135549 rs2836974 21_19 0.438 4260.45 0.00540 -5.16
1135586 rs34578707 21_19 0.285 4259.41 0.00350 -5.13
1135599 rs77090950 21_19 0.278 4260.48 0.00340 -5.15
577665 rs1176746 11_67 1.000 1141.65 0.00330 -2.94
577667 rs2307599 11_67 1.000 1141.79 0.00330 -2.75
1135603 rs35560196 21_19 0.269 4260.46 0.00330 -5.14
72145 rs569546056 2_19 1.000 832.92 0.00240 3.68
1135620 rs8129147 21_19 0.195 4259.42 0.00240 -5.18
1135621 rs28360661 21_19 0.181 4257.33 0.00220 -5.17
392514 rs11406602 7_49 0.222 3149.98 0.00200 -4.15
72144 rs7562170 2_19 0.833 801.87 0.00190 3.52
72148 rs4580350 2_19 0.811 803.08 0.00190 -3.46
422120 rs758184196 8_11 1.000 642.54 0.00190 -2.30
422115 rs2428 8_11 1.000 592.65 0.00170 8.64
834739 rs11090617 22_19 0.999 594.73 0.00170 30.16
1135537 rs34672724 21_19 0.137 4253.23 0.00170 -5.16
422136 rs13265731 8_11 0.873 580.93 0.00150 8.08
844023 rs4989532 1_6 1.000 498.41 0.00140 -3.70
844024 rs2072735 1_6 0.983 495.86 0.00140 -4.12
848095 rs35130213 1_19 1.000 461.64 0.00130 2.49
987654 rs112574791 8_94 0.992 465.69 0.00130 -21.20
1135787 rs12482979 21_19 0.124 3639.38 0.00130 -5.91
1066876 rs373230966 17_5 1.000 412.37 0.00120 3.07
532213 rs17882431 10_64 0.752 511.94 0.00110 -24.55
844028 rs2072734 1_6 0.790 487.79 0.00110 -4.24
1135754 rs8128827 21_19 0.096 3655.88 0.00100 -5.86
1135759 rs2026267 21_19 0.092 3655.38 0.00098 -5.86
1135760 rs35123057 21_19 0.092 3655.36 0.00098 -5.86
1135768 rs2735306 21_19 0.090 3654.71 0.00096 -5.86
1135782 rs2735309 21_19 0.088 3645.36 0.00093 -5.88
1135784 rs2776309 21_19 0.084 3643.45 0.00089 -5.88
1135774 rs2776307 21_19 0.083 3649.26 0.00088 -5.87
1135781 rs2776308 21_19 0.082 3645.78 0.00087 -5.88
493805 rs2763193 9_59 1.000 279.81 0.00081 -17.31
1020932 rs950878 11_53 0.007 37650.52 0.00078 -15.44
1093989 rs58542926 19_15 0.927 272.55 0.00073 19.31
1135800 rs3216359 21_19 0.069 3604.73 0.00072 -5.96
241704 rs4552481 4_95 0.988 239.74 0.00069 16.59
848101 rs34563832 1_19 0.532 448.89 0.00069 2.81
#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
834739 rs11090617 22_19 0.999 594.73 1.7e-03 30.16
834742 rs1977081 22_19 0.001 584.16 1.1e-06 29.82
834749 rs13056555 22_19 0.000 573.44 4.0e-09 29.50
834746 rs2401512 22_19 0.000 572.04 2.9e-09 29.48
834747 rs4823176 22_19 0.000 572.12 2.9e-09 29.48
834748 rs4823178 22_19 0.000 572.21 3.0e-09 29.48
834745 rs2072905 22_19 0.000 571.94 2.8e-09 29.47
834744 rs1883348 22_19 0.000 562.66 6.1e-10 29.26
532213 rs17882431 10_64 0.752 511.94 1.1e-03 -24.55
532198 rs35614792 10_64 0.037 506.82 5.4e-05 -24.45
532232 rs34539738 10_64 0.000 460.22 1.4e-09 -23.32
532231 rs138052038 10_64 0.000 447.22 1.1e-09 -23.05
532224 rs568600628 10_64 0.000 403.80 6.6e-10 -21.98
987654 rs112574791 8_94 0.992 465.69 1.3e-03 -21.20
987602 rs35968570 8_94 0.008 457.58 1.1e-05 -20.97
532189 rs528153341 10_64 0.000 304.20 5.4e-10 -19.45
1093989 rs58542926 19_15 0.927 272.55 7.3e-04 19.31
532191 rs2104368 10_64 0.000 318.15 3.7e-10 -19.24
532192 rs11190401 10_64 0.000 316.50 3.4e-10 -19.20
532185 rs6584345 10_64 0.000 295.46 3.5e-10 -19.18
532187 rs2902371 10_64 0.000 294.80 3.5e-10 -19.16
532188 rs11595928 10_64 0.000 294.55 3.4e-10 -19.15
1094037 rs200210321 19_15 0.044 267.00 3.4e-05 19.14
532219 rs4917889 10_64 0.000 317.96 3.1e-10 -19.09
1094020 rs8107974 19_15 0.028 264.34 2.1e-05 19.09
1094070 rs10401969 19_15 0.004 262.10 2.7e-06 19.01
834754 rs738491 22_19 0.000 286.32 9.7e-12 18.97
1094088 rs739846 19_15 0.002 259.68 1.4e-06 18.93
1094182 rs73001065 19_15 0.004 263.31 3.0e-06 18.86
532177 rs11595257 10_64 0.000 277.63 3.2e-10 -18.62
1093951 rs72999033 19_15 0.004 261.11 2.7e-06 18.61
532167 rs34633631 10_64 0.000 276.29 3.2e-10 -18.59
532171 rs11599683 10_64 0.000 275.88 3.1e-10 -18.58
1094418 rs73002956 19_15 0.007 247.71 5.3e-06 18.36
1094235 rs56255430 19_15 0.003 246.01 2.4e-06 18.35
1094472 rs3794991 19_15 0.006 243.94 4.2e-06 18.25
532174 rs555100159 10_64 0.000 242.52 2.5e-10 -17.95
1094269 rs150268548 19_15 0.001 231.35 6.6e-07 17.93
457384 rs2980875 8_83 0.159 133.94 6.2e-05 -17.90
457380 rs2001844 8_83 0.312 134.16 1.2e-04 -17.85
457382 rs2980886 8_83 0.315 134.18 1.2e-04 -17.85
532184 rs2862991 10_64 0.000 250.52 2.4e-10 -17.82
532166 rs3829161 10_64 0.000 255.42 3.3e-10 -17.45
1094593 rs17216525 19_15 0.001 217.78 8.6e-07 17.42
532164 rs12782078 10_64 0.000 245.33 3.1e-10 -17.40
532165 rs11597528 10_64 0.000 245.47 3.1e-10 -17.40
532178 rs2862949 10_64 0.000 247.20 3.9e-10 -17.38
987776 rs189120668 8_94 0.000 326.65 1.6e-15 -17.33
1094589 rs16996148 19_15 0.001 215.08 7.3e-07 17.33
493805 rs2763193 9_59 1.000 279.81 8.1e-04 -17.31
#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] 32
if (length(genes)>0){
GO_enrichment <- enrichr(genes, dbs)
for (db in dbs){
print(db)
df <- GO_enrichment[[db]]
df <- df[df$Adjusted.P.value<0.05,c("Term", "Overlap", "Adjusted.P.value", "Genes")]
print(df)
}
#DisGeNET enrichment
# devtools::install_bitbucket("ibi_group/disgenet2r")
library(disgenet2r)
disgenet_api_key <- get_disgenet_api_key(
email = "wesleycrouse@gmail.com",
password = "uchicago1" )
Sys.setenv(DISGENET_API_KEY= disgenet_api_key)
res_enrich <-disease_enrichment(entities=genes, vocabulary = "HGNC",
database = "CURATED" )
df <- res_enrich@qresult[1:10, c("Description", "FDR", "Ratio", "BgRatio")]
print(df)
#WebGestalt enrichment
library(WebGestaltR)
background <- ctwas_gene_res$genename
#listGeneSet()
databases <- c("pathway_KEGG", "disease_GLAD4U", "disease_OMIM")
enrichResult <- WebGestaltR(enrichMethod="ORA", organism="hsapiens",
interestGene=genes, referenceGene=background,
enrichDatabase=databases, interestGeneType="genesymbol",
referenceGeneType="genesymbol", isOutput=F)
print(enrichResult[,c("description", "size", "overlap", "FDR", "database", "userId")])
}
Uploading data to Enrichr... Done.
Querying GO_Biological_Process_2021... Done.
Querying GO_Cellular_Component_2021... Done.
Querying GO_Molecular_Function_2021... Done.
Parsing results... Done.
[1] "GO_Biological_Process_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)
[1] "GO_Cellular_Component_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)
[1] "GO_Molecular_Function_2021"
[1] Term Overlap Adjusted.P.value Genes
<0 rows> (or 0-length row.names)
SYTL1 gene(s) from the input list not found in DisGeNET CURATEDZFP62 gene(s) from the input list not found in DisGeNET CURATEDMFSD3 gene(s) from the input list not found in DisGeNET CURATEDADAM15 gene(s) from the input list not found in DisGeNET CURATEDPPP5C gene(s) from the input list not found in DisGeNET CURATEDDDX19A gene(s) from the input list not found in DisGeNET CURATEDPLEKHA3 gene(s) from the input list not found in DisGeNET CURATEDFES gene(s) from the input list not found in DisGeNET CURATEDFPR3 gene(s) from the input list not found in DisGeNET CURATEDTMEM167B gene(s) from the input list not found in DisGeNET CURATEDWHAMM gene(s) from the input list not found in DisGeNET CURATEDC18orf8 gene(s) from the input list not found in DisGeNET CURATEDEVA1C gene(s) from the input list not found in DisGeNET CURATEDCAMTA2 gene(s) from the input list not found in DisGeNET CURATEDHIST1H2BN gene(s) from the input list not found in DisGeNET CURATED
Description FDR Ratio BgRatio
64 Splenic Neoplasms 0.01616818 1/17 1/9703
74 Ataxia, Spinocerebellar 0.01616818 2/17 34/9703
80 Malignant neoplasm of spleen 0.01616818 1/17 1/9703
105 Sulfite oxidase deficiency 0.01616818 1/17 1/9703
115 Woolly hair nevus 0.01616818 1/17 1/9703
137 Spinocerebellar Ataxia Type 1 0.01616818 2/17 34/9703
138 Spinocerebellar Ataxia Type 2 0.01616818 2/17 34/9703
139 Spinocerebellar Ataxia Type 4 0.01616818 2/17 35/9703
140 Spinocerebellar Ataxia Type 5 0.01616818 2/17 34/9703
141 Spinocerebellar Ataxia Type 6 (disorder) 0.01616818 2/17 34/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