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/20190429_RNASeqSTAR_quantifications.Rmd/
    Untracked:  docs/figure/20190502_Check_eQTLs.Rmd/

Unstaged changes:
    Deleted:    analysis/20190412_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 in all samples (~17,500 genes passed this filter). 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 5 genotype PCs (PCs 1-3 visually segregate admixture and large population structure while PCs4-5 take into account some closely related samples). \(x\) is coded as 0,1,2, \(U \sim MVN(0,\sigma^2 K)\) where \(K\) an identity matrix (therefore this is not really and lmm). I choose this model after coworkers suggested to not mix lmm (with kinship matrix) with genotype PCs because it is sort of odd to justify or unelegent. Dealing with both close relatedness and distant population structure in gwas models is an area of ongoing research. Some helpful links: https://www.biorxiv.org/content/biorxiv/early/2018/09/07/409953.full.pdf https://www.nature.com/articles/srep06874

I think in the future I will still refine this model, specifically to include both kinship matrix (lmm) as well as first 3 genotype PCs since the kinship matrix produced by KING or GEMMA clearly does not capture the population structure that is captured in the first 3 PCs that is also reflected in Admixture analysis.

Association testing was implemented in the R package ‘MatrixEQTL’. This resulted in ~5000 eSNP-gene pairs at FDR<0.1 (Benjamini Hodgeberg correction).

Here I want to check that the results of this analysis are reasonable, starting by checking boxplots of gene expression stratified by genotype for a handful of significant eSNP-gene pairs. Based on a previous analysis of phenotypes after regressing out the first 3 RNA-seq PCs (same as what were included in this model) results in sample 317 being an outlier. Here I want to check if this sample is now driving a lot of spurious associations.

# Read in genotypes for eQTLs
Genotypes <- read.table("../data/PastAnalysesDataToKeep/20140424_sig_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] 4699

length(unique(colnames(Genotypes)))

[1] 4699

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

# Read in phenotypes, from count table
CountTable <- read.table('../data/PastAnalysesDataToKeep/ForAssociationTesting.phenotypes.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))

# 5235 eQTLs at FDR<0.1
dim(eQTLs)

[1] 5235    6

# betas for these
hist(eQTLs$beta)

Past versions of checks-1.png

# 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)
}

MyBoxplot(MergedData, "FID", "ID.15.23454749.C.CA", "ENSPTRT00000076477.1")

Past versions of checks-2.png

MyBoxplot(MergedData, "FID", "ID.1.27687107.C.A", "ENSPTRT00000101617.1")

Past versions of checks-3.png

MyBoxplot(MergedData, "FID", "ID.19.47916617.CA.C", "ENSPTRT00000026036.6")

Past versions of checks-4.png

Interpretation: This is concerning; the top hits seem to be all false positives, driven by outliers. Let’s check some randomly sampled eqtls to see if all of the eqtls are like this.

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])))
  )
}

Past versions of checks-for-random-eqtls-1.png

Past versions of checks-for-random-eqtls-2.png

Past versions of checks-for-random-eqtls-3.png

Past versions of checks-for-random-eqtls-4.png

Past versions of checks-for-random-eqtls-5.png

Past versions of checks-for-random-eqtls-6.png

Past versions of checks-for-random-eqtls-7.png

Past versions of checks-for-random-eqtls-8.png

Past versions of checks-for-random-eqtls-9.png

Past versions of checks-for-random-eqtls-10.png

Past versions of checks-for-random-eqtls-11.png

Past versions of checks-for-random-eqtls-12.png

Error in data.frame(Genotype = DataFrame[[SNP.name]], Phenotype = DataFrame[[Gene.name]],  : 
  arguments imply differing number of rows: 0, 38
Error in data.frame(Genotype = DataFrame[[SNP.name]], Phenotype = DataFrame[[Gene.name]],  : 
  arguments imply differing number of rows: 0, 38

Past versions of checks-for-random-eqtls-13.png

Past versions of checks-for-random-eqtls-14.png

Past versions of checks-for-random-eqtls-15.png

Error in data.frame(Genotype = DataFrame[[SNP.name]], Phenotype = DataFrame[[Gene.name]],  : 
  arguments imply differing number of rows: 0, 38

Past versions of checks-for-random-eqtls-16.png

Past versions of checks-for-random-eqtls-17.png

Interpretation: better than previous analysis but still a lot of spurious looking associations are driven by outliers. At least it doesn’t seem like 317 is the outlier in most of these as I would have expected based on a of looking at the residuals after regressing out the same 3 RNA-seq PCs that were included in this model. Note that the outliers seem to be the really really low expressed values.

Next I may try the method of transformating the phenotype data by first standardizing across genes, then quantile normalizing across individuals as described in leafcutter paper.

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