GTEx V6 using pipeline

Last updated: 2019-10-02

Checks: 7 0

Knit directory: GTEx/

This reproducible R Markdown analysis was created with workflowr (version 1.4.0). The Checks tab describes the reproducibility checks that were applied when the results were created. The Past versions tab lists the development history.

R Markdown file: up-to-date

Great! Since the R Markdown file has been committed to the Git repository, you know the exact version of the code that produced these results.

Environment: empty

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.

Seed: set.seed(20181220)

The command set.seed(20181220) 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.

Session information: recorded

Great job! Recording the operating system, R version, and package versions is critical for reproducibility.

Cache: none

Nice! There were no cached chunks for this analysis, so you can be confident that you successfully produced the results during this run.

File paths: relative

Great job! Using relative paths to the files within your workflowr project makes it easier to run your code on other machines.

Repository version: 0526bcb

Great! You are using Git for version control. Tracking code development and connecting the code version to the results is critical for reproducibility. The version displayed above was the version of the Git repository at the time these results were generated.

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:


Ignored files:
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    dsc-mash-gtex/

Untracked files:
    Untracked:  .DS_Store
    Untracked:  analysis/.DS_Store
    Untracked:  code/Demo_SumstatQuery.R
    Untracked:  code/SexDE.TAU.R
    Untracked:  code/SexDE_get_genes.R
    Untracked:  data/.DS_Store
    Untracked:  data/cor_tissues_non_ash_voom_pearson.rda
    Untracked:  data/gene_names_GTEX_V6.txt
    Untracked:  data/genewide_ash_out_tissue_mat_halfuniform_non_mode.rda
    Untracked:  data/order_index.rda
    Untracked:  data/samples_id.txt
    Untracked:  data/sexde/
    Untracked:  data/tissuewide_pearson_halfuniform_tissuewide_non_mode.rda
    Untracked:  output/.DS_Store
    Untracked:  output/GTExV6/
    Untracked:  output/GTExV6pipeline/
    Untracked:  output/corshrink_noise_gene_1.rds
    Untracked:  output/sexde/

Unstaged changes:
    Modified:   analysis/SexDEPlots.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 R Markdown and HTML files. If you’ve configured a remote Git repository (see ?wflow_git_remote), click on the hyperlinks in the table below to view them.

File	Version	Author	Date	Message
Rmd	0526bcb	zouyuxin	2019-10-02	wflow_publish(“analysis/GTExV6pipeline.Rmd”)
html	933fc27	zouyuxin	2019-02-08	Build site.
Rmd	77b789d	zouyuxin	2019-02-08	wflow_publish(“analysis/GTExV6pipeline.Rmd”)
html	58e443e	zouyuxin	2019-01-29	Build site.
Rmd	8607dea	zouyuxin	2019-01-29	wflow_publish(“analysis/GTExV6pipeline.Rmd”)
html	5d17b16	zouyuxin	2019-01-29	Build site.
Rmd	cadbb28	zouyuxin	2019-01-29	wflow_publish(“analysis/GTExV6pipeline.Rmd”)
html	02da57c	zouyuxin	2019-01-27	Build site.
Rmd	378989c	zouyuxin	2019-01-27	wflow_publish(“analysis/GTExV6pipeline.Rmd”)

library(mashr)

Loading required package: ashr

library(knitr)
library(kableExtra)
library(ggplot2)
library(gridExtra)

gtex <- readRDS(gzcon(url("https://github.com/stephenslab/gtexresults/blob/master/data/MatrixEQTLSumStats.Portable.Z.rds?raw=TRUE")))
missing.tissues <- c(7, 8, 19, 20, 24, 25, 31, 34, 37)
gtex.colors <- read.table("https://github.com/stephenslab/gtexresults/blob/master/data/GTExColors.txt?raw=TRUE", sep = '\t', comment.char = '')[-missing.tissues, 2]
gtex.colors <- as.character(gtex.colors)
gene.names = as.character(read.table('data/gene_names.txt')[,1])

The results are from mashr_flashr_pipeline. We include the data driven covariance matrices based on the first three principal components and factors from flash.

Flash model based on effects:

factors = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.flash.model.rds')$factors
par(mfrow = c(2, 3))
for(k in 1:13){
  barplot(factors[,k], col=gtex.colors, names.arg = FALSE, axes = FALSE, main=paste0("Factor ", k))
}

Version	Author	Date
02da57c	zouyuxin	2019-01-27

Version	Author	Date
02da57c	zouyuxin	2019-01-27

Version	Author	Date
02da57c	zouyuxin	2019-01-27

Flash model based on z scores:

factors = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.flash.model.rds')$factors
par(mfrow = c(2, 3))
for(k in 1:18){
  barplot(factors[,k], col=gtex.colors, names.arg = FALSE, axes = FALSE, main=paste0("Factor ", k))
}

Version	Author	Date
02da57c	zouyuxin	2019-01-27

Version	Author	Date
02da57c	zouyuxin	2019-01-27

Version	Author	Date
02da57c	zouyuxin	2019-01-27

# read model
m_simple_EE = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.mash_model_V_simple.rds')
m_simple_EE$result = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.mash_model_V_simple.posterior.rds')
m_simple_EZ = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.mash_model_V_simple.rds')
m_simple_EZ$result = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.mash_model_V_simple.posterior.rds')

m_mle_EE = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.mash_model_V_mle.rds')
m_mle_EE$result = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.mash_model_V_mle.posterior.rds')
m_mle_EZ = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.mash_model_V_mle.rds')
m_mle_EZ$result = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.mash_model_V_mle.posterior.rds')

m_Vgene_EE_kushal = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.V_corshrink_xcondition_kushal.mash_model.rds')
m_Vgene_EE_kushal$result = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.V_corshrink_xcondition_kushal.posterior.rds')
m_Vgene_EZ_kushal = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.V_corshrink_xcondition_kushal.mash_model.rds')
m_Vgene_EZ_kushal$result = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.V_corshrink_xcondition_kushal.posterior.rds')

m_Vgene_EE_simple_corshrink = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.V_corshrink_xcondition_nullz.mash_model.rds')
m_Vgene_EE_simple_corshrink$result = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.V_corshrink_xcondition_nullz.posterior.rds')
m_Vgene_EZ_simple_corshrink = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.V_corshrink_xcondition_nullz.mash_model.rds')
m_Vgene_EZ_simple_corshrink$result = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.V_corshrink_xcondition_nullz.posterior.rds')

m_Vgene_EE_simple = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.V_simple_specific_nullz_step_1.mash_model.rds')
m_Vgene_EE_simple$result = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.V_simple_specific_nullz_step_2.posterior.rds')

m_Vgene_EZ_simple = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.V_simple_specific_nullz_step_1.mash_model.rds')
m_Vgene_EZ_simple$result = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.V_simple_specific_nullz_step_2.posterior.rds')

Estimated null cor V

V.simple = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.V_simple.rds')
corrplot::corrplot(V.simple, method='color', type='upper', tl.col="black", tl.srt=45, tl.cex = 0.5, diag = FALSE, col=colorRampPalette(c("blue", "white", "red"))(200), cl.lim = c(-1,1), title = 'Simple', mar=c(0,0,5,0))

Version	Author	Date
02da57c	zouyuxin	2019-01-27

# dev.off()

V.mle.EE = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EE.FL_PC3.V_mle.rds')
corrplot::corrplot(V.mle.EE, method='color', type='upper', tl.col="black", tl.srt=45, tl.cex = 0.5, diag = FALSE, col=colorRampPalette(c("blue", "white", "red"))(200), cl.lim = c(-1,1), title = 'MLE EE', mar=c(0,0,5,0))

Version	Author	Date
02da57c	zouyuxin	2019-01-27

V.mle.EZ = readRDS('output/GTExV6pipeline/MatrixEQTLSumStats.Portable.Z.EZ.FL_PC3.V_mle.rds')
corrplot::corrplot(V.mle.EZ, method='color', type='upper', tl.col="black", tl.srt=45, tl.cex = 0.5, diag = FALSE, col=colorRampPalette(c("blue", "white", "red"))(200), cl.lim = c(-1,1), title = 'MLE EZ', mar=c(0,0,5,0))

Version	Author	Date
02da57c	zouyuxin	2019-01-27

Results

Vgene Kushal uses gene-specific correlation matrices, the results are from Kushal’s paper (CorShrink on z score matrix for each gene, which include signals).
Vgene simple uses gene-specific correlation matrices, the correlation matrix for each gene is the correlation of z scores close to null.
Vgene simple corshrink uses gene-specific correlation matrix, the correlation matrix for each gene is the CorShrink result using Vgene simple.

Log likelihoods:

logliks = c(get_loglik(m_simple_EE), get_loglik(m_mle_EE), get_loglik(m_Vgene_EE_kushal), get_loglik(m_Vgene_EE_simple), get_loglik(m_Vgene_EE_simple_corshrink))
logliks_EZ = c(get_loglik(m_simple_EZ), get_loglik(m_mle_EZ), get_loglik(m_Vgene_EZ_kushal), get_loglik(m_Vgene_EZ_simple), get_loglik(m_Vgene_EZ_simple_corshrink))
tmp = cbind(logliks, logliks_EZ)
row.names(tmp) = c('Simple', 'MLE', 'Vgene Kushal', 'Vgene simple', 'Vgene simple corshrink')
colnames(tmp) = c('EE', 'EZ')
tmp %>% kable() %>% kable_styling()

	EE	EZ
Simple	936478.4	937254.7
MLE	940058.8	940457.4
Vgene Kushal	886368.9	907004.4
Vgene simple	1001931.7	1011877.2
Vgene simple corshrink	1021859.8	1036059.9

Estimated prior weights:ß

par(mfrow=c(1,2))
barplot(get_estimated_pi(m_simple_EE), las=2, cex.names = 0.7, main = 'Simple EE')
barplot(get_estimated_pi(m_mle_EE), las=2, cex.names = 0.7, main = 'MLE EE')

Version	Author	Date
02da57c	zouyuxin	2019-01-27

barplot(get_estimated_pi(m_Vgene_EE_kushal), las=2, cex.names = 0.7, main = 'V gene specific EE Kushal')
barplot(get_estimated_pi(m_Vgene_EE_simple), las=2, cex.names = 0.7, main = 'V gene specific EE simple')

Version	Author	Date
933fc27	zouyuxin	2019-02-08
5d17b16	zouyuxin	2019-01-29
02da57c	zouyuxin	2019-01-27

barplot(get_estimated_pi(m_Vgene_EE_simple_corshrink), las=2, cex.names = 0.7, main = 'V gene specific EE simple corshrink')

barplot(get_estimated_pi(m_simple_EZ), las=2, cex.names = 0.7, main = 'Simple EZ')

Version	Author	Date
933fc27	zouyuxin	2019-02-08
5d17b16	zouyuxin	2019-01-29

barplot(get_estimated_pi(m_mle_EZ), las=2, cex.names = 0.7, main = 'MLE EZ')
barplot(get_estimated_pi(m_Vgene_EZ_kushal), las=2, cex.names = 0.7, main = 'V gene specific EZ Kushal')

Version	Author	Date
933fc27	zouyuxin	2019-02-08
5d17b16	zouyuxin	2019-01-29

barplot(get_estimated_pi(m_Vgene_EZ_simple), las=2, cex.names = 0.7, main = 'V gene specific EZ simple')
barplot(get_estimated_pi(m_Vgene_EZ_simple_corshrink), las=2, cex.names = 0.7, main = 'V gene specific EZ simple corshrink')

Version	Author	Date
933fc27	zouyuxin	2019-02-08

Number of significant:

numsig_EE = c(length(get_significant_results(m_simple_EE)), 
              length(get_significant_results(m_mle_EE)), 
              length(get_significant_results(m_Vgene_EE_kushal)), 
              length(get_significant_results(m_Vgene_EE_simple)),
              length(get_significant_results(m_Vgene_EE_simple_corshrink)))
numsig_EZ = c(length(get_significant_results(m_simple_EZ)), 
              length(get_significant_results(m_mle_EZ)),
              length(get_significant_results(m_Vgene_EZ_kushal)),
              length(get_significant_results(m_Vgene_EZ_simple)),
              length(get_significant_results(m_Vgene_EZ_simple_corshrink)))
tmp = cbind(numsig_EE, numsig_EZ)
row.names(tmp) = c('Simple', 'MLE', 'Vgene Kushal', 'Vgene simple', 'Vgene simple corshrink')
colnames(tmp) = c('EE', 'EZ')
tmp %>% kable() %>% kable_styling()

	EE	EZ
Simple	13068	13519
MLE	12654	12986
Vgene Kushal	15767	16066
Vgene simple	15684	15838
Vgene simple corshrink	15916	15967

The gene significant in Simple EZ, not in MLE EZ:

ind = setdiff(get_significant_results(m_simple_EZ), get_significant_results(m_mle_EZ))[9]
stronggene = data.frame(gtex$strong.b[ind,])
colnames(stronggene) = 'EffectSize'
stronggene$Group = row.names(stronggene)
stronggene$se = gtex$strong.s[ind,]
p1 = ggplot(stronggene, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + coord_flip() + ggtitle(paste0(gene.names[ind], ' raw')) + ylim(c(-1,1)) + geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneSimple = data.frame(m_simple_EZ$result$PosteriorMean[ind,])
colnames(stronggeneSimple) = 'EffectSize'
stronggeneSimple$Group = row.names(stronggeneSimple)
stronggeneSimple$se = m_simple_EZ$result$PosteriorSD[ind,]
p2 = ggplot(stronggeneSimple, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + coord_flip() + ggtitle(paste0(gene.names[ind],' Simple EZ')) + ylim(c(-1,1)) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneMLE = data.frame(m_mle_EZ$result$PosteriorMean[ind,])
colnames(stronggeneMLE) = 'EffectSize'
stronggeneMLE$Group = row.names(stronggeneMLE)
stronggeneMLE$se = m_mle_EZ$result$PosteriorSD[ind,]
p3 = ggplot(stronggeneMLE, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1,1)) + coord_flip() + ggtitle(paste0(gene.names[ind],' MLE EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneVgeneKushal = data.frame(m_Vgene_EZ_kushal$result$PosteriorMean[ind,])
colnames(stronggeneVgeneKushal) = 'EffectSize'
stronggeneVgeneKushal$Group = row.names(stronggeneVgeneKushal)
stronggeneVgeneKushal$se = m_Vgene_EZ_kushal$result$PosteriorSD[ind,]
p4 = ggplot(stronggeneVgeneKushal, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1,1)) + coord_flip() + ggtitle(paste0(gene.names[ind],' V gene specific Kushal EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneVgeneSimple = data.frame(m_Vgene_EZ_simple$result$PosteriorMean[ind,])
colnames(stronggeneVgeneSimple) = 'EffectSize'
stronggeneVgeneSimple$Group = row.names(stronggeneVgeneSimple)
stronggeneVgeneSimple$se = m_Vgene_EZ_simple$result$PosteriorSD[ind,]
p5 = ggplot(stronggeneVgeneSimple, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1,1)) + coord_flip() + ggtitle(paste0(gene.names[ind],' V gene specific simple EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneVgeneSimpleCor = data.frame(m_Vgene_EZ_simple_corshrink$result$PosteriorMean[ind,])
colnames(stronggeneVgeneSimpleCor) = 'EffectSize'
stronggeneVgeneSimpleCor$Group = row.names(stronggeneVgeneSimpleCor)
stronggeneVgeneSimpleCor$se = m_Vgene_EZ_simple_corshrink$result$PosteriorSD[ind,]
p6 = ggplot(stronggeneVgeneSimpleCor, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1,1)) + coord_flip() + ggtitle(paste0(gene.names[ind],' V gene specific simple corshrink EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 2)

Version	Author	Date
933fc27	zouyuxin	2019-02-08
5d17b16	zouyuxin	2019-01-29
02da57c	zouyuxin	2019-01-27

The gene MCPH1:

ind = 13837
stronggene = data.frame(gtex$strong.b[13837,])
colnames(stronggene) = 'EffectSize'
stronggene$Group = row.names(stronggene)
stronggene$se = gtex$strong.s[13837,]
p1 = ggplot(stronggene, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + coord_flip() + ggtitle('ENSG00000249898 row') + ylim(c(-1.3,1.1)) + geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneSimple = data.frame(m_simple_EZ$result$PosteriorMean[13837,])
colnames(stronggeneSimple) = 'EffectSize'
stronggeneSimple$Group = row.names(stronggeneSimple)
stronggeneSimple$se = m_simple_EZ$result$PosteriorSD[13837,]
p2 = ggplot(stronggeneSimple, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1.3,1.1)) + coord_flip() + ggtitle('ENSG00000249898 Simple EZ') + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneMLE = data.frame(m_mle_EZ$result$PosteriorMean[13837,])
colnames(stronggeneMLE) = 'EffectSize'
stronggeneMLE$Group = row.names(stronggeneMLE)
stronggeneMLE$se = m_mle_EZ$result$PosteriorSD[13837,]
p3 = ggplot(stronggeneMLE, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + coord_flip() + ggtitle('ENSG00000249898 MLE EZ') + ylim(c(-1.3,1.1)) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneVgeneKushal = data.frame(m_Vgene_EZ_kushal$result$PosteriorMean[ind,])
colnames(stronggeneVgeneKushal) = 'EffectSize'
stronggeneVgeneKushal$Group = row.names(stronggeneVgeneKushal)
stronggeneVgeneKushal$se = m_Vgene_EZ_kushal$result$PosteriorSD[ind,]
p4 = ggplot(stronggeneVgeneKushal, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1.3,1.1)) + coord_flip() + ggtitle(paste0(gene.names[ind],' V gene specific Kushal EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneVgeneSimple = data.frame(m_Vgene_EZ_simple$result$PosteriorMean[ind,])
colnames(stronggeneVgeneSimple) = 'EffectSize'
stronggeneVgeneSimple$Group = row.names(stronggeneVgeneSimple)
stronggeneVgeneSimple$se = m_Vgene_EZ_simple$result$PosteriorSD[ind,]
p5 = ggplot(stronggeneVgeneSimple, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1,1)) + coord_flip() + ggtitle(paste0(gene.names[ind],' V gene specific simple EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneVgeneSimpleCor = data.frame(m_Vgene_EZ_simple_corshrink$result$PosteriorMean[ind,])
colnames(stronggeneVgeneSimpleCor) = 'EffectSize'
stronggeneVgeneSimpleCor$Group = row.names(stronggeneVgeneSimpleCor)
stronggeneVgeneSimpleCor$se = m_Vgene_EZ_simple_corshrink$result$PosteriorSD[ind,]
p6 = ggplot(stronggeneVgeneSimpleCor, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1,1)) + coord_flip() + ggtitle(paste0(gene.names[ind],' V gene specific simple corshrink EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))


grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 2)

Version	Author	Date
933fc27	zouyuxin	2019-02-08
5d17b16	zouyuxin	2019-01-29
02da57c	zouyuxin	2019-01-27

The gene significant in V gene specific EZ (simple corshrink) tissuewide, not in mle EZ:

ind = setdiff(get_significant_results(m_Vgene_EZ_simple_corshrink), get_significant_results(m_mle_EZ))[10]
stronggene = data.frame(gtex$strong.b[ind,])
colnames(stronggene) = 'EffectSize'
stronggene$Group = row.names(stronggene)
stronggene$se = gtex$strong.s[ind,]
p1 = ggplot(stronggene, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + coord_flip() + ggtitle(paste0(gene.names[ind],' row')) + ylim(c(-1.3,1.4)) + geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneSimple = data.frame(m_simple_EZ$result$PosteriorMean[ind,])
colnames(stronggeneSimple) = 'EffectSize'
stronggeneSimple$Group = row.names(stronggeneSimple)
stronggeneSimple$se = m_simple_EZ$result$PosteriorSD[ind,]
p2 = ggplot(stronggeneSimple, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1.3,1.4)) + coord_flip() + ggtitle(paste0(gene.names[ind],' Simple EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneMLE = data.frame(m_mle_EZ$result$PosteriorMean[ind,])
colnames(stronggeneMLE) = 'EffectSize'
stronggeneMLE$Group = row.names(stronggeneMLE)
stronggeneMLE$se = m_mle_EZ$result$PosteriorSD[ind,]
p3 = ggplot(stronggeneMLE, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + coord_flip() + ggtitle(paste0(gene.names[ind],' MLE EZ')) + ylim(c(-1.3,1.4)) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneVgeneKushal = data.frame(m_Vgene_EZ_kushal$result$PosteriorMean[ind,])
colnames(stronggeneVgeneKushal) = 'EffectSize'
stronggeneVgeneKushal$Group = row.names(stronggeneVgeneKushal)
stronggeneVgeneKushal$se = m_Vgene_EZ_kushal$result$PosteriorSD[ind,]
p4 = ggplot(stronggeneVgeneKushal, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1,1)) + coord_flip() + ggtitle(paste0(gene.names[ind],' V gene specific Kushal EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneVgeneSimple = data.frame(m_Vgene_EZ_simple$result$PosteriorMean[ind,])
colnames(stronggeneVgeneSimple) = 'EffectSize'
stronggeneVgeneSimple$Group = row.names(stronggeneVgeneSimple)
stronggeneVgeneSimple$se = m_Vgene_EZ_simple$result$PosteriorSD[ind,]
p5 = ggplot(stronggeneVgeneSimple, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1,1)) + coord_flip() + ggtitle(paste0(gene.names[ind],' V gene specific simple EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))

stronggeneVgeneSimpleCor = data.frame(m_Vgene_EZ_simple_corshrink$result$PosteriorMean[ind,])
colnames(stronggeneVgeneSimpleCor) = 'EffectSize'
stronggeneVgeneSimpleCor$Group = row.names(stronggeneVgeneSimpleCor)
stronggeneVgeneSimpleCor$se = m_Vgene_EZ_simple_corshrink$result$PosteriorSD[ind,]
p6 = ggplot(stronggeneVgeneSimpleCor, aes(y = EffectSize, x = Group)) + 
  geom_point(show.legend = FALSE, color=gtex.colors) + ylim(c(-1,1)) + coord_flip() + ggtitle(paste0(gene.names[ind],' V gene specific simple corshrink EZ')) + 
  geom_errorbar(aes(ymin=EffectSize-1.96*se, ymax=EffectSize+1.96*se), width=0.4, show.legend = FALSE, color=gtex.colors) + 
  theme_bw(base_size=12) + theme(axis.text.y = element_text(colour = gtex.colors, size = 6))


grid.arrange(p1, p2, p3, p4, p5, p6, nrow = 2)

Version	Author	Date
933fc27	zouyuxin	2019-02-08
5d17b16	zouyuxin	2019-01-29
02da57c	zouyuxin	2019-01-27

The pairwise sharing by magnitude

par(mfrow = c(1,2))
clrs=colorRampPalette(rev(c('darkred', 'red','orange','yellow','cadetblue1', 'cyan', 'dodgerblue4', 'blue','darkorchid1','lightgreen','green', 'forestgreen','darkolivegreen')))(200)

x           <- get_pairwise_sharing(m_simple_EZ)
colnames(x) <- colnames(get_lfsr(m_simple_EZ))
rownames(x) <- colnames(x)

corrplot::corrplot(x, method='color', type='upper', tl.col="black", tl.srt=45, tl.cex = 0.7, diag = FALSE, col=clrs, cl.lim = c(0,1), title = 'Simple EZ', mar=c(0,0,5,0))

x           <- get_pairwise_sharing(m_mle_EZ)
colnames(x) <- colnames(get_lfsr(m_mle_EZ))
rownames(x) <- colnames(x)
corrplot::corrplot(x, method='color', type='upper', tl.col="black", tl.srt=45, tl.cex = 0.7, diag = FALSE, col=clrs, cl.lim = c(0,1), title = 'MLE EZ', mar=c(0,0,5,0))

Version	Author	Date
58e443e	zouyuxin	2019-01-29
5d17b16	zouyuxin	2019-01-29

par(mfrow = c(1,2))
clrs=colorRampPalette(rev(c('darkred', 'red','orange','yellow','cadetblue1', 'cyan', 'dodgerblue4', 'blue','darkorchid1','lightgreen','green', 'forestgreen','darkolivegreen')))(200)

x           <- get_pairwise_sharing(m_Vgene_EZ_kushal)
colnames(x) <- colnames(get_lfsr(m_Vgene_EZ_kushal))
rownames(x) <- colnames(x)

corrplot::corrplot(x, method='color', type='upper', tl.col="black", tl.srt=45, tl.cex = 0.7, diag = FALSE, col=clrs, cl.lim = c(0,1), title = 'V gene specific Kushal EZ', mar=c(0,0,5,0))

x           <- get_pairwise_sharing(m_Vgene_EZ_simple)
colnames(x) <- colnames(get_lfsr(m_Vgene_EZ_simple))
rownames(x) <- colnames(x)
corrplot::corrplot(x, method='color', type='upper', tl.col="black", tl.srt=45, tl.cex = 0.7, diag = FALSE, col=clrs, cl.lim = c(0,1), title = 'V gene specific simple EZ', mar=c(0,0,5,0))

Version	Author	Date
933fc27	zouyuxin	2019-02-08
5d17b16	zouyuxin	2019-01-29

x           <- get_pairwise_sharing(m_Vgene_EZ_simple_corshrink)
colnames(x) <- colnames(get_lfsr(m_Vgene_EZ_simple_corshrink))
rownames(x) <- colnames(x)
corrplot::corrplot(x, method='color', type='upper', tl.col="black", tl.srt=45, tl.cex = 0.7, diag = FALSE, col=clrs, cl.lim = c(0,1), title = 'V gene specific simple corshrink EZ', mar=c(0,0,5,0))

Version	Author	Date
933fc27	zouyuxin	2019-02-08

sessionInfo()

R version 3.5.3 (2019-03-11)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Mojave 10.14.6

Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.5/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.5/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] gridExtra_2.3     ggplot2_3.2.0     kableExtra_1.1.0  knitr_1.23       
[5] mashr_0.2.21.0641 ashr_2.2-38      

loaded via a namespace (and not attached):
 [1] tidyselect_0.2.5  xfun_0.7          corrplot_0.84    
 [4] purrr_0.3.2       lattice_0.20-38   colorspace_1.4-1 
 [7] htmltools_0.3.6   viridisLite_0.3.0 yaml_2.2.0       
[10] rlang_0.4.0       mixsqp_0.1-119    pillar_1.4.2     
[13] withr_2.1.2       glue_1.3.1        foreach_1.4.4    
[16] plyr_1.8.4        stringr_1.4.0     munsell_0.5.0    
[19] gtable_0.3.0      workflowr_1.4.0   rvest_0.3.4      
[22] mvtnorm_1.0-10    codetools_0.2-16  evaluate_0.13    
[25] labeling_0.3      pscl_1.5.2        doParallel_1.0.14
[28] parallel_3.5.3    highr_0.8         Rcpp_1.0.1       
[31] readr_1.3.1       backports_1.1.4   scales_1.0.0     
[34] rmeta_3.0         webshot_0.5.1     truncnorm_1.0-8  
[37] abind_1.4-5       fs_1.3.1          hms_0.4.2        
[40] digest_0.6.19     stringi_1.4.3     dplyr_0.8.1      
[43] grid_3.5.3        rprojroot_1.3-2   tools_3.5.3      
[46] magrittr_1.5      lazyeval_0.2.2    tibble_2.1.3     
[49] crayon_1.3.4      whisker_0.3-2     pkgconfig_2.0.2  
[52] MASS_7.3-51.1     Matrix_1.2-15     SQUAREM_2017.10-1
[55] xml2_1.2.0        assertthat_0.2.1  rmarkdown_1.13   
[58] httr_1.4.0        rstudioapi_0.10   iterators_1.0.10 
[61] R6_2.4.0          git2r_0.25.2      compiler_3.5.3

GTEx V6 using pipeline

Yuxin Zou

2018-12-20

Flash model based on effects:

Flash model based on z scores:

Estimated null cor V

Results