mashR: tissue-specific effects

gene-SNP pair

Previous pipeline: https://github.com/RajLabMSSM/mashR

GTEx paper: https://www.nature.com/articles/s41588-018-0268-8

https://stephenslab.github.io/gtexresults/Tspecific.html

mash_rds = readRDS(paste0(fastqtl_to_mash_output, "QTLSumStats.mash.rds"))
# names(mash_rds)

maxb <- mash_rds$strong.b
maxz <- mash_rds$strong.z

mash_posterior = readRDS(paste0(mashr_flashr_workflow_output, "QTLSumStats.mash.EZ.FL_PC3.V_simple.posterior.rds"))
# names(mash_posterior)
# head(mash_posterior$lfsr)
# dim(mash_posterior$lfsr)

pm.mash <- mash_posterior$PosteriorMean
lfsr.all <- mash_posterior$lfsr
standard.error <- maxb/maxz
pm.mash.beta <- pm.mash*standard.error

#pm.mash.beta <- pm.mash.beta[rowSums(lfsr.all<0.10)>0,]

Number of TS

Threshold = lfsr < 0.05

# Compute number of “tissue-specific” effects for each tissue
# “tissue-specific” to mean that the effect is at least 2-fold larger in one tissue than in any other tissue. 

source(paste0(code_folder, "normfuncs.R"))
thresh <- 0.05

lfsr.mash = lfsr.all
nsig                <- rowSums(lfsr.mash < thresh) # ensembl | number of columns its significative
pm.mash.beta.norm   <- het.norm(effectsize = pm.mash.beta)
pm.mash.beta.norm   <- pm.mash.beta.norm[nsig > 0,]
lfsr.mash           <- as.matrix(lfsr.mash[nsig > 0,])
colnames(lfsr.mash) <- colnames(maxz)
a         <- which(rowSums(pm.mash.beta.norm > 0.5) == 1)
lfsr.fold <- as.matrix(lfsr.mash[a,])
pm        <- as.matrix(pm.mash.beta.norm[a,])
tspec     <- NULL
for(i in 1:ncol(pm))
  tspec[i] <- sum(pm[,i] > 0.5)
tspec           <- as.matrix(tspec)
rownames(tspec) <- colnames(maxz)
rownames(tspec) <- gsub("(.*?)_(.*)", "\\1", rownames(tspec))
colnames(tspec) <- c("number_TS")
createDT(tspec)

Barplot of TS

Plot number of “tissue-specific” effects for each tissue. “Tissue-specific” to mean that the effect is at least 2-fold larger in one tissue than in any other tissue.

# Plot number of “tissue-specific” effects for each tissue 
colnames(lfsr.fold) <- gsub("(.*?)_(.*)", "\\1", colnames(lfsr.fold))
tissue.colors = pal_futurama()(4)

barplot(as.numeric(t(tspec)),las = 2,cex.names = 0.75,col = tissue.colors,
        names = colnames(lfsr.fold))

TS effect total

Number of total “tissue-specific” effects:

# Unique gene-SNP pair 
dim(lfsr.fold)[1]

[1] 1791

TS effects lists

What are the gene-SNP pairs from the barplot?

colnames(pm) <- gsub("(.*?)_(.*)", "\\1", colnames(pm))
gene_snp_list = list()

for(i in 1:ncol(pm))
  gene_snp_list[[colnames(pm)[i]]] <- rownames(pm)[pm[,i] > 0.5]

## Get conversion table for Gencode 30
gencode_30 = read.table("~/pd-omics/katia/ens.geneid.gencode.v30", header = T)
gencode_30$ensembl = gsub("(.*?)\\.(.*)", "\\1", gencode_30$gene_id)

MFG TS list

MFG = as.data.frame(gene_snp_list$MFG)
colnames(MFG) = c("gene_snp")
MFG$ensembl = gsub("(.*?)_(.*)", "\\1", MFG$gene_snp)
MFG$gene_snp = gsub("(.*?)_(.*)", "\\2", MFG$gene_snp)
MFG$snp = gsub("s", "rs", MFG$gene_snp)
MFG$gene_snp = NULL
MFG_symbol = merge(MFG, gencode_30, by = "ensembl")
MFG_symbol$gene_id = NULL 

write.table(MFG_symbol, file = paste0(work_dir, "ts_effect_5per/ts_effect_2fold_MFG.txt"), quote = F, row.names = F, sep = "\t")
createDT(MFG_symbol)

STG TS list

STG = as.data.frame(gene_snp_list$STG)
colnames(STG) = c("gene_snp")
STG$ensembl = gsub("(.*?)_(.*)", "\\1", STG$gene_snp)
STG$gene_snp = gsub("(.*?)_(.*)", "\\2", STG$gene_snp)
STG$snp = gsub("s", "rs", STG$gene_snp)
STG$gene_snp = NULL
STG_symbol = merge(STG, gencode_30, by = "ensembl")
STG_symbol$gene_id = NULL 

write.table(STG_symbol, file = paste0(work_dir, "ts_effect_5per/ts_effect_2fold_STG.txt"), quote = F, row.names = F, sep = "\t")
createDT(STG_symbol)

THA TS list

THA = as.data.frame(gene_snp_list$THA)
colnames(THA) = c("gene_snp")
THA$ensembl = gsub("(.*?)_(.*)", "\\1", THA$gene_snp)
THA$gene_snp = gsub("(.*?)_(.*)", "\\2", THA$gene_snp)
THA$snp = gsub("s", "rs", THA$gene_snp)
THA$gene_snp = NULL
THA_symbol = merge(THA, gencode_30, by = "ensembl")
THA_symbol$gene_id = NULL 

write.table(THA_symbol, file = paste0(work_dir, "ts_effect_5per/ts_effect_2fold_THA.txt"), quote = F, row.names = F, sep = "\t")
createDT(THA_symbol)

SVZ TS list

SVZ = as.data.frame(gene_snp_list$SVZ)
colnames(SVZ) = c("gene_snp")
SVZ$ensembl = gsub("(.*?)_(.*)", "\\1", SVZ$gene_snp)
SVZ$gene_snp = gsub("(.*?)_(.*)", "\\2", SVZ$gene_snp)
SVZ$snp = gsub("s", "rs", SVZ$gene_snp)
SVZ$gene_snp = NULL
SVZ_symbol = merge(SVZ, gencode_30, by = "ensembl")
SVZ_symbol$gene_id = NULL 

write.table(SVZ_symbol, file = paste0(work_dir, "ts_effect_5per/ts_effect_2fold_SVZ.txt"), quote = F, row.names = F, sep = "\t")
createDT(SVZ_symbol)

sessionInfo()

R version 3.6.1 (2019-07-05) Platform: x86_64-apple-darwin15.6.0 (64-bit) Running under: macOS Catalina 10.15.4

Matrix products: default BLAS: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib LAPACK: /Library/Frameworks/R.framework/Versions/3.6/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] ggsci_2.9 rmeta_3.0 dplyr_0.8.5 lattice_0.20-38

loaded via a namespace (and not attached): [1] Rcpp_1.0.4 later_1.0.0 pillar_1.4.3 compiler_3.6.1
[5] tools_3.6.1 digest_0.6.25 jsonlite_1.6.1 evaluate_0.14
[9] tibble_2.1.3 lifecycle_0.2.0 gtable_0.3.0 pkgconfig_2.0.3
[13] rlang_0.4.5 shiny_1.4.0 crosstalk_1.0.0 yaml_2.2.0
[17] xfun_0.11 fastmap_1.0.1 stringr_1.4.0 knitr_1.26
[21] htmlwidgets_1.5.1 grid_3.6.1 DT_0.13 tidyselect_0.2.5 [25] glue_1.4.0 R6_2.4.1 rmarkdown_2.0 ggplot2_3.3.0
[29] purrr_0.3.3 magrittr_1.5 promises_1.1.0 scales_1.1.0
[33] htmltools_0.4.0 assertthat_0.2.1 xtable_1.8-4 mime_0.8
[37] colorspace_1.4-1 httpuv_1.5.2 stringi_1.4.6 munsell_0.5.0
[41] crayon_1.3.4