PWAS Schema

Plasma_Protein_EA_hg38.pos is a text file with proteins (ID) and locations on their encoding genes, (CHR, P0, P1), as well as pointers to their weights files (WGT).

The weights file for each protein is in Plasma_Protein_weights_EA, in TWAS/Fusion format. When loaded, each .RDat file contains snps (snp info), wgt.matrix (weights), and cv.performance (cross validation) data. The columns of the snps table are chromosome (V1), rsid (V2), position (V4), effect allele (V5) and reference allele (V6). In the wgt.matrix table, the rownames are the rsids, and the columns are the weights derived from elastic net and top1 methods for each snp.

Load Libraries

Run in R:

suppressPackageStartupMessages(library(RSQLite))
suppressPackageStartupMessages(library(tidyverse))
suppressPackageStartupMessages(library(data.table))
suppressPackageStartupMessages(library(biomaRt))

Initialize Extra Table

The file Plasma_Protein_EA_hg38.pos points to the TWAS/FUSION files that contains the weights for each gene. We create a dictionary file that joins gene name, TWAS/FUSION file name, and Ensembl ID.

ensembl = useEnsembl(biomart = "genes", dataset = "hsapiens_gene_ensembl")

Ensembl site unresponsive, trying useast mirror

gene_annot = getBM( attributes=
                    c("ensembl_gene_id",
                      "hgnc_symbol"),
                  values =TRUE,
                  mart = ensembl)
gene_annot = gene_annot[!duplicated(gene_annot$hgnc_symbol),]
gene_list = read.table("/Users/sabrinami/Github/ARIC/PWAS/PWAS_AA/Plasma_Protein_AA_hg38.pos", head=TRUE)
gene_list = gene_list[2:6]

Some of the genes in the PWAS do not have an Ensembl ID annotation, so we use HGNC symbol in its place.

extra = left_join(gene_list, gene_annot, by=c("ID"="hgnc_symbol"))
extra = extra %>% mutate(ensembl_gene_id = coalesce(ensembl_gene_id,ID))

Then add columns to match PrediXcan format.

extra$pred.perf.R2 <- NA
extra$pred.perf.pval <- NA
extra$pred.perf.qval <- NA

Convert File to Dataframe

make_df will load a file and store its data as a dataframe. This is only for a single gene, so later will be repeated for all genes. The input is the name of the .RDat file, and it returns returns dataframe with gene, position, chromosome, ref allele, eff allele, and non-zero enet weights.

make_df <- function(file) {
  if (file %in% extra$WGT) {
    load(file)  
    weights <- data.frame(wgt.matrix) 
    snps <- data.frame(snps) 
    
    index = which(extra$WGT == file)
    
    weights$gene <- extra$ensembl_gene_id[index]
    weights$rsid <- rownames(weights)
    weights$varID <- paste("chr",paste(snps$V1,snps$V4,snps$V6,snps$V5,"b38", sep="_"), sep="")
    weights$ref_allele <- snps$V6
    weights$eff_allele <- snps$V5
    weights = filter(weights, enet != 0)[,c(3,4,5,6,7,1)]
    rownames(weights) <- c()
    weights = rename(weights, c("weight"="enet"))
    
    rsq = cv.performance[1,1]
    pval = cv.performance[2,1]
    extra$pred.perf.R2[index] = rsq
    extra$pred.perf.pval[index] = pval
    assign('extra',extra,envir=.GlobalEnv)
    return(weights)
  }
}

Make Weights Table

First, combine .RDat file names in a vector. Then convert each of them to a dataframe in PrediXcan format, appending them in a weights table.

setwd("/Users/sabrinami/Github/ARIC/PWAS/PWAS_AA/Plasma_Protein_weights_AA")
files <- list.files(pattern = "\\.RDat")

weights = data.frame()
for(i in 1:length(files)) {
  weights <- rbind(weights, make_df(files[i]))
}

Make Extra Table

Generate number of snps for each gene from the weights table.

extra = rename(extra, c("genename"="ID", "gene"="ensembl_gene_id"))
n.snps = weights %>% group_by(gene) %>% summarise(n.snps.in.model = n())

`summarise()` ungrouping output (override with `.groups` argument)

extra = inner_join(extra, n.snps)

Joining, by = "gene"

extra <- extra[,c(6,2,10,7,8,9)]

Write to SQLite Database

Create database connection, and write the weights and extra tables to database.

model_db = "/Users/sabrinami/Github/ARIC/models/ARIC_AA_hg38.db"
conn <- dbConnect(RSQLite::SQLite(), model_db)
dbWriteTable(conn, "weights", weights, overwrite=TRUE)
dbWriteTable(conn, "extra", extra, overwrite=TRUE)

To double check, confirm there is a weights and extra table, and show their contents.

dbListTables(conn)

[1] "extra"   "weights"

dbGetQuery(conn, 'SELECT * FROM weights') %>% head

             gene       rsid                  varID ref_allele eff_allele
1 ENSG00000254521  rs3752135 chr19_51497370_T_G_b38          T          G
2 ENSG00000254521  rs3826667 chr19_51500820_C_T_b38          C          T
3 ENSG00000254521  rs3810114 chr19_51503097_T_C_b38          T          C
4 ENSG00000254521 rs10418495 chr19_51503534_G_A_b38          G          A
5 ENSG00000254521  rs8113048 chr19_51503888_T_G_b38          T          G
6 ENSG00000254521  rs8112579 chr19_51504162_C_T_b38          C          T
        weight
1 -0.083643830
2 -0.076623212
3 -0.031837335
4 -0.008567506
5 -0.031837782
6 -0.008003640

dbGetQuery(conn, 'SELECT * FROM extra') %>% head

             gene genename n.snps.in.model  pred.perf.R2 pred.perf.pval
1 ENSG00000254521 SIGLEC12               8  0.1066109691   6.020845e-48
2 ENSG00000111361   EIF2B1              20  0.0473916727   9.995404e-22
3 ENSG00000198931     APRT              54  0.1955522022   1.229389e-90
4 ENSG00000111674     ENO2             103  0.0253777708   2.520872e-12
5 ENSG00000168610    STAT3              23  0.0901864931   1.617896e-40
6 ENSG00000089127     OAS1              14 -0.0005060883   8.161129e-01
  pred.perf.qval
1             NA
2             NA
3             NA
4             NA
5             NA
6             NA

Lastly, disconnect from database connection

dbDisconnect(conn)

Session information

sessionInfo()

R version 4.0.3 (2020-10-10)
Platform: x86_64-apple-darwin17.0 (64-bit)
Running under: macOS Big Sur 10.16

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] biomaRt_2.44.4    data.table_1.13.2 forcats_0.5.0     stringr_1.4.0    
 [5] dplyr_1.0.2       purrr_0.3.4       readr_1.4.0       tidyr_1.1.2      
 [9] tibble_3.0.4      ggplot2_3.3.2     tidyverse_1.3.0   RSQLite_2.2.1    

loaded via a namespace (and not attached):
 [1] Biobase_2.48.0       httr_1.4.2           bit64_4.0.5         
 [4] jsonlite_1.7.1       modelr_0.1.8         assertthat_0.2.1    
 [7] askpass_1.1          BiocFileCache_1.12.1 highr_0.8           
[10] stats4_4.0.3         blob_1.2.1           cellranger_1.1.0    
[13] yaml_2.2.1           progress_1.2.2       pillar_1.4.6        
[16] backports_1.1.10     glue_1.4.2           digest_0.6.27       
[19] promises_1.1.1       rvest_0.3.6          colorspace_1.4-1    
[22] htmltools_0.5.0      httpuv_1.5.4         XML_3.99-0.5        
[25] pkgconfig_2.0.3      broom_0.7.2          haven_2.3.1         
[28] scales_1.1.1         later_1.1.0.1        git2r_0.27.1        
[31] openssl_1.4.3        generics_0.0.2       IRanges_2.22.2      
[34] ellipsis_0.3.1       withr_2.3.0          BiocGenerics_0.34.0 
[37] cli_2.1.0            magrittr_1.5         crayon_1.3.4        
[40] readxl_1.3.1         memoise_1.1.0        evaluate_0.14       
[43] fs_1.5.0             fansi_0.4.1          xml2_1.3.2          
[46] tools_4.0.3          prettyunits_1.1.1    hms_0.5.3           
[49] lifecycle_0.2.0      S4Vectors_0.26.1     munsell_0.5.0       
[52] reprex_0.3.0         AnnotationDbi_1.50.3 compiler_4.0.3      
[55] rlang_0.4.8          grid_4.0.3           rstudioapi_0.11     
[58] rappdirs_0.3.1       rmarkdown_2.5        gtable_0.3.0        
[61] curl_4.3             DBI_1.1.0            R6_2.4.1            
[64] lubridate_1.7.9      knitr_1.30           bit_4.0.4           
[67] workflowr_1.6.2      rprojroot_1.3-2      stringi_1.5.3       
[70] parallel_4.0.3       Rcpp_1.0.5           vctrs_0.3.4         
[73] dbplyr_1.4.4         tidyselect_1.1.0     xfun_0.18