20190424_Check_eQTLs

workflowr

Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/
    Ignored:    analysis/.DS_Store
    Ignored:    analysis_temp/.DS_Store
    Ignored:    data/.DS_Store
    Ignored:    data/PastAnalysesDataToKeep/.DS_Store
    Ignored:    docs/.DS_Store

Untracked files:
    Untracked:  analysis/20190502_Check_eQTLs.Rmd
    Untracked:  analysis_temp/20190412_Check_eQTLs.Rmd
    Untracked:  data/PastAnalysesDataToKeep/20190428_log10TPM.txt.gz
    Untracked:  docs/assets/
    Untracked:  docs/figure/20190426_CheckRNASeqPCs_2.Rmd/
    Untracked:  docs/figure/20190429_RNASeqSTAR_quantifications.Rmd/
    Untracked:  docs/figure/20190502_Check_eQTLs.Rmd/

Unstaged changes:
    Deleted:    analysis/20190412_Check_eQTLs.Rmd
    Modified:   analysis/20190428_Check_eQTLs.Rmd

library(tidyverse)
library(knitr)
library(data.table)

My second pass at eqtl mapping was as follows: Genotypes filtered for MAF>0.5 & GenotypingRate>0.9, HWE-pvalue<10e-7.5. ~9.5M variants passed these filter. Gene expression as log(TPM), filtering for genes with >0 TPM (transcript level TPM summed to get gene TPM) in all samples (~14000 genes passed this filter). TPM was then standardized across individuals and quantile normalized to a normal distribution across genes. Sample MD_And was dropped from analysis because it was previously shown to be an outlier that caused spurious associations. Association testing used the following linear mixed model for each cis-variant-gene-pair (cis definied as <1Mb from gene):

\[ Y =Wα+xβ+u+ε \]

where \(Y\) is gene expression as \(log(TPM)\), \(W\) covariates include first three genotype principal components (to account for population structure) as well as 3 RNA-seq PCs, Sex, an interceptm and 3 genotype PCs. \(x\) is coded as 0,1,2, \(U \sim MVN(0,\sigma^2 K)\) where \(K\) an centered kinship matrix made from gemma software.

Association testing was implemented in the R package ‘MatrixEQTL’. This resulted in 0 eSNP-gene pairs at FDR<0.1 (Benjamini Hodgeberg correction). QQ-plots didn’t show any P-value inflation relative to the null. From previous (unshown) tests, this lack of significant tests is likely due to the quantile normalizatino step diminishing signal. Even though none of the results are significant after multiple hypothesis correction, I still want to check some of the top nominally significant results to check that they aren’t false positives driven by extreme outliers as I had seen in previous analyses.

# Read in genotypes for eQTLs
Genotypes <- read.table("../data/PastAnalysesDataToKeep/20190428_28_Top5000eqtls.genotypes.raw", header=T, check.names = F, stringsAsFactors = F)
colnames(Genotypes) <- sub("_.*", "", colnames(Genotypes))
# Genotypes[!duplicated(as.list(Genotypes))]
kable(Genotypes[1:10,1:10])

#Make sure there aren't duplicate columns
length(colnames(Genotypes))

[1] 4890

length(unique(colnames(Genotypes)))

[1] 4890

# Read in eQTLs from MatrixEQTL output (already filtered for FDR<0.1)
eQTLs <- read.table("../data/PastAnalysesDataToKeep/20190428_Top5000eqtls.txt.gz", header=T)
kable(head(eQTLs))

# This count table is the log10(TPM); no standardization or quantile normalization
# Read in phenotypes, from count table
CountTable <- read.table('../output/log10TPM.StandardizedAndNormalized.txt', header=T, check.names=FALSE, row.names = 1) %>% 
  t() %>%
  as.data.frame() %>%
  rownames_to_column(var = "FID")
kable(CountTable[1:10, 1:10])

MergedData <- left_join(Genotypes, CountTable, by="FID")

#eqtls, ordered from most significant at top
kable(head(eQTLs))

# betas for these
hist(eQTLs$beta, breaks=20)

Not sure what to make of that multi-modal distribution of betas (Two peaks around 0 might be expected, but I wasn’t expecting another peak at ~5). Maybe most of these eqtls are the snps in LD that are associated with the same gene…

#how many snps
length(unique(eQTLs$SNP))

[1] 4941

#how many genes
length(unique(eQTLs$gene))

[1] 1103

Average of ~5 snps per gene, so I’m still perplexed by the multiple modes in the distribution of effect sizes.

In any case, now I want to check some boxplots

# Expression box plot, stratified by genotype. For a fe of top of the list snp-gene pairs (most significant)
MyBoxplot <- function(DataFrame, Labels.name, SNP.name, Gene.name){
  data.frame(Genotype = DataFrame[[SNP.name]],
    Phenotype = DataFrame[[Gene.name]],
    FID=DataFrame[[Labels.name]]) %>%
    ggplot(aes(x=factor(Genotype), y=Phenotype, label=FID)) +
    geom_boxplot(outlier.shape = NA) +
    geom_text(position=position_jitter(width=0.25), alpha=1, size=2) +
    scale_y_continuous(name=paste("log TPM", Gene.name)) +
    xlab(SNP.name)
}

#Box plot of top3 eqtls
MyBoxplot(MergedData, "FID", "ID.3.10249149.C.CT", "ENSPTRG00000043259")

MyBoxplot(MergedData, "FID", "ID.19.55922879.C.CT", "ENSPTRG00000011430")

MyBoxplot(MergedData, "FID", "ID.22.4961539.C.A", "ENSPTRG00000014048")

Interpretation: These look more believable than previous analyses without quantile normalization in the sense that these don’t look driven by single outliers. Let’s check some randomly sampled eqtls.

set.seed(1)
RandomSampleOfEqtls <- eQTLs %>% sample_n(20) %>% select(SNP, gene, beta)
kable(RandomSampleOfEqtls)

for(i in 1:nrow(RandomSampleOfEqtls)) {
  try(
    print(MyBoxplot(MergedData, "FID", as.character(RandomSampleOfEqtls$SNP[i]), as.character(RandomSampleOfEqtls$gene[i])))
  )
}

Session information

sessionInfo()

R version 3.5.1 (2018-07-02)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS  10.14

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] data.table_1.12.0 knitr_1.22        forcats_0.4.0    
 [4] stringr_1.4.0     dplyr_0.8.0.1     purrr_0.3.2      
 [7] readr_1.3.1       tidyr_0.8.2       tibble_2.1.1     
[10] ggplot2_3.1.0     tidyverse_1.2.1  

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1       highr_0.8        cellranger_1.1.0 plyr_1.8.4      
 [5] pillar_1.3.1     compiler_3.5.1   git2r_0.24.0     workflowr_1.2.0 
 [9] tools_3.5.1      digest_0.6.18    lubridate_1.7.4  jsonlite_1.6    
[13] evaluate_0.13    nlme_3.1-137     gtable_0.3.0     lattice_0.20-38 
[17] pkgconfig_2.0.2  rlang_0.3.3      cli_1.1.0        rstudioapi_0.10 
[21] yaml_2.2.0       haven_2.1.0      xfun_0.6         withr_2.1.2     
[25] xml2_1.2.0       httr_1.4.0       hms_0.4.2        generics_0.0.2  
[29] fs_1.2.6         rprojroot_1.3-2  grid_3.5.1       tidyselect_0.2.5
[33] glue_1.3.1       R6_2.4.0         readxl_1.1.0     rmarkdown_1.11  
[37] modelr_0.1.4     magrittr_1.5     whisker_0.3-2    backports_1.1.3 
[41] scales_1.0.0     htmltools_0.3.6  rvest_0.3.2      assertthat_0.2.1
[45] colorspace_1.4-1 labeling_0.3     stringi_1.4.3    lazyeval_0.2.2  
[49] munsell_0.5.0    broom_0.5.1      crayon_1.3.4