20190502_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/20190428_Check_eQTLs.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
    Modified:   analysis/index.Rmd

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

My third 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(CPM), filtering for genes with 80% of samples > 10 reads on the gene. CPM 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 4 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’. I can see clear enrichment of small P-values compared to a single-pass label-permutated null. (Insert images here).

For this analysis I somewhat arbitraily choose a P-value threshold to further examine hits of 1e-6 since this is where I see deviation from the permutated null. At this threshold, I estimate FDR~20% (num of hits in real data versus permuted.)

Permuted null.

Real data.

# Read in genotypes for eQTLs
Genotypes <- read.table("../data/PastAnalysesDataToKeep/20190502_SigQTLs.genotypes.txt.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] 289

length(unique(colnames(Genotypes)))

[1] 289

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

SampleList <- read.table("../output/ForAssociationTesting.temp.fam")$V2

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

MergedData <- left_join(Genotypes, CountTable, by=c("IID" = "FID"))  %>% 
  as.data.frame()

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

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

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

[1] 303

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

[1] 43

# 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("log10(CPM)", Gene.name)) +
    xlab(SNP.name)
}

MyBoxplot(MergedData, "IID", as.character("ID.1.126465756.T.C"), as.character("ENSPTRG00000001061"))

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

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

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

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     backports_1.1.3  scales_1.0.0    
[41] htmltools_0.3.6  rvest_0.3.2      assertthat_0.2.1 colorspace_1.4-1
[45] labeling_0.3     stringi_1.4.3    lazyeval_0.2.2   munsell_0.5.0   
[49] broom_0.5.1      crayon_1.3.4