Check RNA-seq PCs

workflowr

Ignored files:
    Ignored:    .Rproj.user/
    Ignored:    code/snakemake_workflow/bahblah

Untracked files:
    Untracked:  ._workflowr.yml.swp
    Untracked:  analysis/20190320_Check-RNAseq-PCs.Rmd
    Untracked:  docs/figure/

Unstaged changes:
    Modified:   _workflowr.yml
    Modified:   analysis/index.Rmd

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

# Read in count table, filtering for TPM>1
CountTable <- read.table('../output/CountTable.tpm.txt', header=T, check.names=FALSE, row.names = 1) %>%
  filter(rowSums(.) > 1)

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

Now perform PCA, plot a few visualizations

# pca with log-transformed count table (+ 0.1 pseudocount)
pca_results <- prcomp(t(log(CountTable+0.1)), center=T, scale. = T)

# Plot PCs
autoplot(pca_results)

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

# Plot a couple PCs with a couple potential covariates
ggplot(Merged, aes(x=PC2, y=PC1, color=factor(RNA.Extract_date), label=Row.names)) + 
  geom_point() +
  geom_text_repel(size=2.5)

ggplot(Merged, aes(x=PC2, y=PC1, color=RIN, label=Row.names)) + 
  geom_point() +
  geom_text_repel(size=2.5)

Now I am going to look more systematically for significant correlations between potential observed confounders in the Metadata and the first 10 PCs. Will use Spearman’s correlation to test continuous continuous confounders, will use anova for categorical confounders.

# Grab first 10 PCs
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")]
kable(head(Continuous_confounders_to_test))

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

# Plot
Spearman_test_results$p %>%
  melt() %>%
  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, 1))) +
    scale_fill_gradient(limits=c(0.001, 1), breaks=c(0.001,0.01,0.1,1), trans = 'log', high="white", low="red" ) +
    theme_classic()

# 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
Pvalues %>%
  melt() %>%
  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, 1))) +
    scale_fill_gradient( trans = 'log', high="white", low="red" ) +
    theme_classic() +
    theme(axis.text.x = element_text(angle = 45, hjust = 1))

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] bindrcpp_0.2.2  reshape2_1.4.3  knitr_1.21      ggrepel_0.8.0  
 [5] psych_1.8.10    forcats_0.4.0   stringr_1.3.1   dplyr_0.7.8    
 [9] purrr_0.2.5     readr_1.3.1     tidyr_0.8.2     tibble_1.4.2   
[13] tidyverse_1.2.1 readxl_1.1.0    ggfortify_0.4.5 ggplot2_3.1.0  
[17] corrplot_0.84  

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