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
    Ignored:    docs/assets/.DS_Store

Untracked files:
    Untracked:  analysis/20190502_Check_eQTLs.Rmd
    Untracked:  analysis_temp/20190412_Check_eQTLs.Rmd
    Untracked:  data/PastAnalysesDataToKeep/20190428_log10TPM.txt.gz

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

library(plyr)
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)

Past versions of checks-1.png

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

Past versions of checks-2.png

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

Past versions of checks-3.png

Past versions of checks-4.png

Past versions of checks-5.png

Past versions of checks-6.png

Past versions of checks-7.png

Past versions of checks-8.png

Past versions of checks-9.png

Past versions of checks-10.png

Past versions of checks-11.png

Past versions of checks-12.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-13.png

Past versions of checks-14.png

Past versions of checks-15.png

Past versions of checks-16.png

Past versions of checks-17.png

Past versions of checks-18.png

Past versions of checks-19.png

Past versions of checks-20.png

# List of chimp tested genes
ChimpTestedGenes <- rownames(read.table('../output/ExpressionMatrix.un-normalized.txt.gz', header=T, check.names=FALSE, row.names = 1))


ChimpToHumanGeneMap <- read.table("../data/Biomart_export.Hsap.Ptro.orthologs.txt.gz", header=T, sep='\t', stringsAsFactors = F)
kable(head(ChimpToHumanGeneMap))

# Of this ortholog list, how many genes are one2one
table(ChimpToHumanGeneMap$Chimpanzee.homology.type)

ortholog_many2many  ortholog_one2many   ortholog_one2one 
              2278              19917             140351 

OneToOneMap <- ChimpToHumanGeneMap %>%
  filter(Chimpanzee.homology.type=="ortholog_one2one")

# Read gtex heart egene list
GtexHeartEgenes <- read.table("../data/Heart_Left_Ventricle.v7.egenes.txt.gz", header=T, sep='\t', stringsAsFactors = F) %>%
  mutate(gene_id_stable = gsub(".\\d+$","",gene_id)) %>%
  filter(gene_id_stable %in% OneToOneMap$Gene.stable.ID) %>%
  mutate(chimp_id = plyr::mapvalues(gene_id_stable, OneToOneMap$Gene.stable.ID,  OneToOneMap$Chimpanzee.gene.stable.ID, warn_missing = F)) %>%
  filter(chimp_id %in% ChimpTestedGenes)


length(GtexHeartEgenes$gene_id_stable)

[1] 11586

length(OneToOneMap$Gene.stable.ID)

[1] 140351

length(intersect(GtexHeartEgenes$gene_id_stable, OneToOneMap$Gene.stable.ID))

[1] 11586

HumanSigGenes <- GtexHeartEgenes %>%
  filter(qval<0.05) %>%
  pull(chimp_id)

HumanNonSigGenes <- GtexHeartEgenes %>%
  filter(qval>0.05) %>%
  pull(chimp_id)

ChimpSigGenes <- GtexHeartEgenes %>%
  filter(chimp_id %in% eQTLs$gene) %>%
  pull(chimp_id)

ChimpNonSigGenes <- GtexHeartEgenes %>%
  filter(! chimp_id %in% eQTLs$gene) %>%
  pull(chimp_id)



matrix( c( length(intersect(ChimpSigGenes,HumanSigGenes)),length(intersect(HumanSigGenes,ChimpNonSigGenes)),
         length(intersect(ChimpSigGenes,HumanNonSigGenes)), length(intersect(ChimpNonSigGenes,HumanNonSigGenes)), nrow = 2))

     [,1]
[1,]   15
[2,] 4225
[3,]   18
[4,] 7328
[5,]    2

ContigencyTable <- matrix(c(15,4225,18,7328), nrow=2)

#Contigency table of one to one orthologs tested in both chimps and humans of whether significant in humans, or chimps, or both, or neither
ContigencyTable

     [,1] [,2]
[1,]   15   18
[2,] 4225 7328

fisher.test(ContigencyTable)

    Fisher's Exact Test for Count Data

data:  ContigencyTable
p-value = 0.2846
alternative hypothesis: true odds ratio is not equal to 1
95 percent confidence interval:
 0.677282 3.040443
sample estimates:
odds ratio 
  1.445263 

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   plyr_1.8.4       

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1       highr_0.8        cellranger_1.1.0 compiler_3.5.1  
 [5] pillar_1.3.1     git2r_0.24.0     workflowr_1.2.0  tools_3.5.1     
 [9] digest_0.6.18    lubridate_1.7.4  jsonlite_1.6     evaluate_0.13   
[13] nlme_3.1-137     gtable_0.3.0     lattice_0.20-38  pkgconfig_2.0.2 
[17] rlang_0.3.3      cli_1.1.0        rstudioapi_0.10  yaml_2.2.0      
[21] haven_2.1.0      xfun_0.6         withr_2.1.2      xml2_1.2.0      
[25] httr_1.4.0       hms_0.4.2        generics_0.0.2   fs_1.2.6        
[29] rprojroot_1.3-2  grid_3.5.1       tidyselect_0.2.5 glue_1.3.1      
[33] R6_2.4.0         readxl_1.1.0     rmarkdown_1.11   modelr_0.1.4    
[37] magrittr_1.5     whisker_0.3-2    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