20190421_RegressOutRNASeqPC

workflowr

Ignored files:
    Ignored:    .DS_Store
    Ignored:    .Rhistory
    Ignored:    .Rproj.user/

Untracked files:
    Untracked:  analysis/20190327_MakeCovariateFiles.Rmd
    Untracked:  analysis/20190421_RegressOutRNASeqPCs.Rmd
    Untracked:  analysis_temp/
    Untracked:  docs/figure/20190412_Check_eQTLs.Rmd/
    Untracked:  docs/figure/20190424_Check_eQTLs.Rmd/

Unstaged changes:
    Deleted:    analysis/20190412_Check-Kinship-Matrices.Rmd
    Modified:   analysis/20190412_Check_eQTLs.Rmd
    Modified:   analysis/index.Rmd

library(corrplot)
library(ggfortify)
library(readxl)
library(tidyverse)
library(psych)
library(ggrepel)
library(knitr)
library(reshape2)
library(gplots)
library(matrixStats)

RNA-seq data for each individual was pseudo-mapped/quantified by kallisto and merged into a matrix of TPM values. Genes were filtered by some criteria (see link) to generate a table of phenotypes for association testing. Here I will assess if there is still an obvious batch affect (compared to my previous PCA analysis of the TPM count table) after I regress out the first (or first few PCs), since I will be included some of these PCs as covariates in the gene expression model.

# Read in count table, filtering out rows that are all zeros. Here the count-table is log transformed TPM values, with genes filtered just as they would be used for association testing
CountTable <- read.table(gzfile('../output/ForAssociationTesting.phenotypes.txt'), header=T, check.names=FALSE, row.names = 1)

kable(head(CountTable))

# Read Covariates file
Covariates <- read.table('../output/Covariates/3RS_10GT.covariates.txt', header=T, check.names = F, stringsAsFactors = F, row.names=1)

# Read admixture coefficients (K=4), and first 3 principle components, since some form of population substructure will likely be included in the expression modeling as a covariate.
AdmixtureCoeff <- read.table("../output/PopulationStructure/Admixture/MergedForAdmixture.4.Q.labelled") %>% 
  dplyr::rename(Individual.ID=V2) %>%
  select(-V1, -V3, -V4, -V5, -V6) %>%
  dplyr::rename(Admix.Western=V9, Admix.Eastern=V10, Admix.Central=V8, Admix.NigeriaCameroon=V7) #Renaming the admixture clusters after looking at plots with known subspecies
kable(head(AdmixtureCoeff))

GenotypePCs <- read.table("../output/PopulationStructure/pca.eigenvec", header=T) %>%
  select(IID, PC1, PC2, PC3) %>%
  dplyr::rename(Individual.ID=IID, GenotypePC1=PC1, GenotypePC2=PC2, GenotypePC3=PC3)
kable(head(GenotypePCs))

# Read in other metadata
OtherMetadata <- as.data.frame(read_excel("../data/Metadata.xlsx"))
kable(head(OtherMetadata))

#Merge all metadata tables
Metadata <- OtherMetadata %>%
  left_join(GenotypePCs, by=c("Individual.ID")) %>%
  left_join(AdmixtureCoeff, by=c("Individual.ID"))
kable(head(Metadata))

Check some of the gene expression phenotype distributions

#function to plot datapoints as violin plot with individual points labelled. I find this slightly more useful than histograms for plotting phenotype because I can see the distribution and also label points
MyPlot <- function(DataIn){
  ggplot(mapping=aes(x=1, y=DataIn, label=names(DataIn))) +
  geom_violin() +
  geom_text(position=position_jitter(width=0.25), alpha=1, size=2)
}

#Plot expression phenotypes for a few genes
MyPlot(unlist(CountTable[1,]))

MyPlot(unlist(CountTable[2,]))

MyPlot(unlist(CountTable[3,]))

MyPlot(unlist(CountTable[4,]))

MyPlot(unlist(CountTable[5,]))

MyPlot(unlist(CountTable[6,]))

MyPlot(unlist(CountTable[7,]))

MyPlot(unlist(CountTable[8,]))

MyPlot(unlist(CountTable[9,]))

MyPlot(unlist(CountTable[10,]))

#What is distribution of expression levels for all genes in the table
# Histogram of rowMeans. (Arithmetic mean of the log-transformed data)
hist(rowMeans(CountTable))

Plot correlation matrix… Should closely resemble the correlation matrix in this previous analysis, except that the previous correlation matrix included the top5000 expressed genes, while this one includes the top5000 genes of the genes that are included in the phenotype table for association testing.

CorMatrix <- CountTable %>%
  mutate(sumVar = rowSums(.)) %>%
  arrange(desc(sumVar)) %>%
  head(5000) %>%
  select(-sumVar) %>%
  scale() %>%
  cor(method = c("spearman"))

RNAExtractionDate <- as.character(unclass(factor(plyr::mapvalues(row.names(CorMatrix), from=Metadata$Individual.ID, to=Metadata$RNA.Extract_date))))

The following `from` values were not present in `x`: MD_And

RNA.Library.prep.batch <- as.character(unclass(factor(plyr::mapvalues(row.names(CorMatrix), from=Metadata$Individual.ID, to=Metadata$RNA.Library.prep.batch))))

The following `from` values were not present in `x`: MD_And

# Heatmap of correlation. Row colors for RNA extraction batch, column colors for RNA library prep batch
heatmap.2(CorMatrix, trace="none", ColSideColors=RNAExtractionDate, RowSideColors = RNA.Library.prep.batch)

# What is mean correlation
mean(CorMatrix)

[1] 0.8059961

Now regress out the first RNA-seq PC and recheck correlation matrix

# Regress out PC1

GetResiduals <- function(Y, X){
  df <- data.frame(c(Y),X)
  mylm <- lm(Y ~ ., data=df)
  return(mylm$residuals)
}

ResidualCountTable <- t(apply(CountTable,1,GetResiduals, t(Covariates[c("PC1.RS"),])))

CorMatrix <- ResidualCountTable %>% as.data.frame() %>%
  mutate(sumVar = rowSums(.)) %>%
  arrange(desc(sumVar)) %>%
  head(5000) %>%
  select(-sumVar) %>%
  scale() %>%
  cor(method = c("spearman"))

RNAExtractionDate <- as.character(unclass(factor(plyr::mapvalues(row.names(CorMatrix), from=Metadata$Individual.ID, to=Metadata$RNA.Extract_date))))

The following `from` values were not present in `x`: MD_And

RNA.Library.prep.batch <- as.character(unclass(factor(plyr::mapvalues(row.names(CorMatrix), from=Metadata$Individual.ID, to=Metadata$RNA.Library.prep.batch))))

The following `from` values were not present in `x`: MD_And

# Heatmap of correlation. Row colors for RNA extraction batch, column colors for RNA library prep batch
heatmap.2(CorMatrix, trace="none", ColSideColors=RNAExtractionDate, RowSideColors = RNA.Library.prep.batch)

# What is mean correlation
mean(CorMatrix)

[1] 0.2167298

Interpretation: The regressing first PC1 seems to take away the obvious batch effect that was present in the lower left of the original correlation matrix (The samples that were in that batch no longer cluster together in the lower left), however, regressing out first PC alone still leaves 2 clear clusters of samples.

In my first pass at eQTL mapping I including 3 PCs, let’s see how the correlation matrix clusters after regressing out all of these PCs. (Not iteratively regressing out, but all at once; take residuals from this model Y~PC1+PC2+PC3)

# Regress out PC1, PC2, and PC3

ResidualCountTable <- t(apply(CountTable,1,GetResiduals, t(Covariates[c("PC1.RS", "PC2.RS", "PC3.RS"),])))

CorMatrix <- ResidualCountTable %>% as.data.frame() %>%
  mutate(sumVar = rowSums(.)) %>%
  arrange(desc(sumVar)) %>%
  head(5000) %>%
  select(-sumVar) %>%
  scale() %>%
  cor(method = c("spearman"))

RNAExtractionDate <- as.character(unclass(factor(plyr::mapvalues(row.names(CorMatrix), from=Metadata$Individual.ID, to=Metadata$RNA.Extract_date))))

The following `from` values were not present in `x`: MD_And

RNA.Library.prep.batch <- as.character(unclass(factor(plyr::mapvalues(row.names(CorMatrix), from=Metadata$Individual.ID, to=Metadata$RNA.Library.prep.batch))))

The following `from` values were not present in `x`: MD_And

# Heatmap of correlation. Row colors for RNA extraction batch, column colors for RNA library prep batch
heatmap.2(CorMatrix, trace="none", ColSideColors=RNAExtractionDate, RowSideColors = RNA.Library.prep.batch)

# What is mean correlation
mean(CorMatrix)

[1] 0.1172888

Clearly there is still residual structure in the data even after regressing out the top 3 PCs. Now I want to figure out if this residual structure is caused by any technical factors that I could justifiably adjust for (eg RIN score) or if it is genetic. I will now perform PCA on the residuals and see if they correlate to any of the observed technical factors or genetic population substructure that I previously looked at before regressing out PCs.

#First get a list of the same genes that were used in the original PCA analysis (top 2500 expressed genes)

GenesForPCA <- CountTable %>%
  rownames_to_column('gene') %>%
  mutate(sumVar = rowSums(select(.,-gene))) %>%
  arrange(desc(sumVar)) %>%
  head(2500) %>%
  pull(gene)

pca_results <- ResidualCountTable %>% as.data.frame() %>%
  rownames_to_column('gene') %>%
  filter(gene %in% GenesForPCA) %>%
  column_to_rownames('gene') %>%
  mutate(Variance = rowVars(as.matrix(.))) %>%
  filter(Variance>0) %>%
  select(-Variance) %>%
  t() %>%
  prcomp(center=T, scale. = T)

summary(pca_results)

Importance of components:
                           PC1     PC2     PC3     PC4     PC5     PC6
Standard deviation     42.6985 8.70349 7.77771 7.20705 6.58620 6.20179
Proportion of Variance  0.7322 0.03042 0.02429 0.02086 0.01742 0.01545
Cumulative Proportion   0.7322 0.76262 0.78691 0.80777 0.82519 0.84064
                           PC7     PC8     PC9    PC10   PC11    PC12
Standard deviation     5.38525 5.23972 5.02747 4.90940 4.7344 4.44945
Proportion of Variance 0.01165 0.01103 0.01015 0.00968 0.0090 0.00795
Cumulative Proportion  0.85228 0.86331 0.87346 0.88314 0.8921 0.90009
                          PC13    PC14   PC15    PC16    PC17    PC18
Standard deviation     4.25279 4.06900 4.0245 3.97194 3.82780 3.73073
Proportion of Variance 0.00726 0.00665 0.0065 0.00634 0.00588 0.00559
Cumulative Proportion  0.90736 0.91401 0.9205 0.92685 0.93273 0.93832
                          PC19    PC20    PC21   PC22    PC23    PC24
Standard deviation     3.70848 3.59657 3.57082 3.4581 3.40354 3.29980
Proportion of Variance 0.00552 0.00519 0.00512 0.0048 0.00465 0.00437
Cumulative Proportion  0.94384 0.94904 0.95416 0.9590 0.96362 0.96799
                          PC25    PC26    PC27    PC28    PC29    PC30
Standard deviation     3.23734 3.12198 3.02463 2.94894 2.85489 2.75365
Proportion of Variance 0.00421 0.00391 0.00367 0.00349 0.00327 0.00305
Cumulative Proportion  0.97220 0.97611 0.97979 0.98328 0.98655 0.98960
                          PC31    PC32    PC33    PC34     PC35      PC36
Standard deviation     2.67728 2.57653 2.50205 2.41617 2.01e-14 7.464e-15
Proportion of Variance 0.00288 0.00267 0.00251 0.00234 0.00e+00 0.000e+00
Cumulative Proportion  0.99248 0.99514 0.99766 1.00000 1.00e+00 1.000e+00
                            PC37      PC38
Standard deviation     5.681e-15 5.184e-15
Proportion of Variance 0.000e+00 0.000e+00
Cumulative Proportion  1.000e+00 1.000e+00

#Scree plot
screeplot(pca_results)

# Merge with metadata
Merged <- merge(pca_results$x, Metadata, by.x = "row.names", by.y = "Individual.ID", all=TRUE)
kable(head(Merged))

PCs_to_test <- Merged[,2:11]
kable(head(PCs_to_test))

# Grab potential continuous confounders that make sense to test
Continuous_confounders_to_test <- Merged[, c("RIN", "Age", "GenotypePC1", "GenotypePC2", "GenotypePC3", "Admix.Western", "Admix.Eastern", "Admix.Central", "Admix.NigeriaCameroon")]
kable(head(Continuous_confounders_to_test))

# Test
Spearman_test_results <- corr.test(Continuous_confounders_to_test, PCs_to_test, adjust="none", method="spearman")

MinP_floor <- floor(log10(min(Spearman_test_results$p)))
# Plot
Spearman_test_results$p %>%
  melt() %>%
  dplyr::rename(Pvalue = value, Principle.Component=Var2, Potential.Confounder=Var1) %>%
  ggplot(aes(x=Potential.Confounder, y=Principle.Component, fill=Pvalue)) + 
    geom_tile() +
    geom_text(aes(label = signif(Pvalue, 2))) +
    scale_fill_gradient(limits=c(10**MinP_floor, 1), breaks=10**seq(MinP_floor,0,1), trans = 'log', high="white", low="red" ) +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

# Grab potential categorical confounders that make sense to test
Categorical_confounders_to_test <- Merged[,c("Viral.status", "SX","RNA.Extract_date", "RNA.Library.prep.batch", "Source")]
kable(head(Categorical_confounders_to_test))

# Viral status will need to be reformatted to make factors that make sense for testing (example: HBV+, HBV- HCV+, HCV- are factors that make sense). Let's assume that NA means negative status.
Categorical_confounders_to_test$HBV_status <- grepl("HBV+", Categorical_confounders_to_test$Viral.status)
Categorical_confounders_to_test$HCV_status <- grepl("HCV+", Categorical_confounders_to_test$Viral.status)
Categorical_confounders_to_test <- Categorical_confounders_to_test[, -1 ]
kable(head(Categorical_confounders_to_test))

# Do one-way anova test as a loop.
# First initialize results matrix
Pvalues <- matrix(ncol = dim(PCs_to_test)[2], nrow = dim(Categorical_confounders_to_test)[2])
colnames(Pvalues) <- colnames(PCs_to_test)
rownames(Pvalues) <- colnames(Categorical_confounders_to_test)
for (confounder in seq_along(Categorical_confounders_to_test)) {
    for (PC in seq_along(PCs_to_test)) {
      res.aov <- aov(PCs_to_test[[PC]] ~ Categorical_confounders_to_test[[confounder]])
      pval <- summary(res.aov)[[1]][["Pr(>F)"]][1]
      Pvalues[confounder, PC] <- pval
    }
}
# Plot
MinP_floor <- floor(log10(min(Pvalues)))
Pvalues %>%
  melt() %>%
  dplyr::rename(Pvalue = value, Principle.Component=Var2, Potential.Confounder=Var1) %>%
  ggplot(aes(x=Potential.Confounder, y=Principle.Component, fill=Pvalue)) + 
    geom_tile() +
    geom_text(aes(label = signif(Pvalue, 2))) +
    scale_fill_gradient(limits=c(10**MinP_floor, 1), breaks=10**seq(MinP_floor,0,1), trans = 'log', high="white", low="red" ) +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

Interpretation: The only observed technical factors (or measured population structure statistics) that correlate with any of the PCs on the residuals are Admix.Central and HBV status. Given that HBV status is correlated with such late PCs that explain so little variance (see scree plot above), and there are only 3 HBV+ observations, I think including HBV in the model is not justified. Also note that the correlation matrix of the residuals after regressing out 3PCs generally shows two clusters, neither of which segregate with the 3HBV+ observations. The other significant correlation is between PC9 (which explains little variance) and Admix.Central coefficient. This trait is already included implicitly in the eQTL calling model as a genotype PC covariate.

Lastly, I want to look at the residuals after regressing out those first three PCs.

#Plot residuals for a few genes
MyPlot(ResidualCountTable[1,])

MyPlot(ResidualCountTable[2,])

MyPlot(ResidualCountTable[3,])

MyPlot(ResidualCountTable[4,])

MyPlot(ResidualCountTable[5,])

MyPlot(ResidualCountTable[6,])

MyPlot(ResidualCountTable[7,])

MyPlot(ResidualCountTable[8,])

MyPlot(ResidualCountTable[9,])

MyPlot(ResidualCountTable[10,])

#Plot of first two PCs of residuals
pca_results$x %>%
  as.data.frame() %>%
  rownames_to_column('sample') %>%
  ggplot(aes(x=PC1, y=PC2, label=sample)) +
  geom_text()

Interpretation: After regressing out the PCs and looking at residuals, sample 317 seems to be consistently an outlier. This is also reflected in PC space of the residuals. It was not an outlier in the original phenotypes. Was it an outlier in the original PCs that were regressed out?

# The PCs were read in from the covariates file

t(Covariates) %>%
  as.data.frame() %>%
  rownames_to_column('sample') %>%
  ggplot(aes(x=PC1.RS, y=PC2.RS, label=sample)) +
  geom_text()

t(Covariates) %>%
  as.data.frame() %>%
  rownames_to_column('sample') %>%
  ggplot(aes(x=PC2.RS, y=PC3.RS, label=sample)) +
  geom_text()

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] matrixStats_0.54.0 gplots_3.0.1       reshape2_1.4.3    
 [4] knitr_1.22         ggrepel_0.8.0      psych_1.8.10      
 [7] forcats_0.4.0      stringr_1.4.0      dplyr_0.8.0.1     
[10] purrr_0.3.2        readr_1.3.1        tidyr_0.8.2       
[13] tibble_2.1.1       tidyverse_1.2.1    readxl_1.1.0      
[16] ggfortify_0.4.5    ggplot2_3.1.0      corrplot_0.84     

loaded via a namespace (and not attached):
 [1] Rcpp_1.0.1         lubridate_1.7.4    lattice_0.20-38   
 [4] gtools_3.8.1       assertthat_0.2.1   rprojroot_1.3-2   
 [7] digest_0.6.18      R6_2.4.0           cellranger_1.1.0  
[10] plyr_1.8.4         backports_1.1.3    evaluate_0.13     
[13] highr_0.8          httr_1.4.0         pillar_1.3.1      
[16] rlang_0.3.3        lazyeval_0.2.2     rstudioapi_0.10   
[19] gdata_2.18.0       rmarkdown_1.11     labeling_0.3      
[22] foreign_0.8-71     munsell_0.5.0      broom_0.5.1       
[25] compiler_3.5.1     modelr_0.1.4       xfun_0.6          
[28] pkgconfig_2.0.2    mnormt_1.5-5       htmltools_0.3.6   
[31] tidyselect_0.2.5   gridExtra_2.3      workflowr_1.2.0   
[34] crayon_1.3.4       withr_2.1.2        bitops_1.0-6      
[37] grid_3.5.1         nlme_3.1-137       jsonlite_1.6      
[40] gtable_0.3.0       git2r_0.24.0       magrittr_1.5      
[43] scales_1.0.0       KernSmooth_2.23-15 cli_1.1.0         
[46] stringi_1.4.3      fs_1.2.6           xml2_1.2.0        
[49] generics_0.0.2     tools_3.5.1        glue_1.3.1        
[52] hms_0.4.2          parallel_3.5.1     yaml_2.2.0        
[55] colorspace_1.4-1   caTools_1.17.1.1   rvest_0.3.2       
[58] haven_2.1.0