Check RNA-seq PCs

workflowr

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

Untracked files:
    Untracked:  analysis/20190321_Check-Kinship-And-PopulationStructure.Rmd
    Untracked:  analysis/20190325_MergingRNASeqLanes.Rmd
    Untracked:  analysis/20190326_Admixture.Rmd
    Untracked:  analysis/20190326_PCA.Rmd
    Untracked:  analysis/20190327_MakeFamAndCovariateFiles.Rmd
    Untracked:  analysis/20190327_MakeFamPhenotypeFile.Rmd
    Untracked:  docs/figure/20190321_Check-Kinship-And-PopulationStructure.Rmd/
    Untracked:  docs/figure/20190325_MergingRNASeqLanes.Rmd/
    Untracked:  docs/figure/20190326_PCA.Rmd/

Unstaged changes:
    Deleted:    ._workflowr.yml.swp
    Modified:   analysis/20190320_Check-RNAseq-PCs.Rmd
    Modified:   analysis/index.Rmd

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

RNA-seq data for each individual was pseudo-mapped/quantified by kallisto and merged into a matrix of TPM values.

# Read in count table, filtering out rows that are all zeros
CountTable <- read.table(gzfile('../output/CountTable.tpm.txt.gz'), header=T, check.names=FALSE, row.names = 1) %>%
  rownames_to_column('gene') %>% #hack to keep rownames despite dplyr filter step
  filter_if(is.numeric, any_vars(. > 0)) %>%
  column_to_rownames('gene')

kable(head(CountTable))

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

First, look at TPM distribution of a couple genes across samples, looking for normality. Plot normal-QQ plot of TPM. With and without a couple potential transformations.

RowNumber<-5000

# raw tpm
Phenotype <- as.numeric(CountTable[RowNumber,])
hist(Phenotype)

qqnorm(Phenotype)
qqline(Phenotype)

# log(tpm)
Phenotype <- log(as.numeric(CountTable[RowNumber,]))
hist(Phenotype)

qqnorm(Phenotype)
qqline(Phenotype)

Now I will more systematically test for normality for all genes after a few transformations…

# histogram of P-values of Shapiro test for normality of raw TPM for each gene.
df.shapiro <- CountTable %>%
  apply(1, shapiro.test)
Pvals<-unlist(lapply(df.shapiro, function(x) x$p.value))
hist(Pvals)

# ...After log-transform+pseudocount.
df.shapiro <- CountTable %>%
  +0.1 %>% log() %>%
  apply(1, shapiro.test)
Pvals<-unlist(lapply(df.shapiro, function(x) x$p.value))
hist(Pvals)

# Maybe the inflation of non-normal expression phenotypes is due to lots of rows with too many (0 + pseudocount) values.
# Retry after filtering all genes with any 0 counts, then log-transform.
df.shapiro <- CountTable %>%
  filter_all(all_vars(.>0)) %>% log() %>%
  apply(1, shapiro.test)
Pvals<-unlist(lapply(df.shapiro, function(x) x$p.value))
hist(Pvals)

# ...After filtering all genes with any 0 counts, then log-transform, and filter for only top 500 expressed genes
df.shapiro <- CountTable %>%
  filter_all(all_vars(.>0)) %>% log() %>%
  mutate(sumVar = rowSums(.)) %>%
  arrange(desc(sumVar)) %>%
  head(500) %>%
  select(-sumVar) %>%
  apply(1, shapiro.test)
Pvals<-unlist(lapply(df.shapiro, function(x) x$p.value))
hist(Pvals)

Data is generally not normal, even after considering only highly expressed genes after log-transformation. Probably too much structure/covariates that must be accounted for. For now I will keep exploring the data with log-transformed data, and consider if/how to more carefully transform the data later (eg quantile normalization) .

Plot correlation matrix… For this purpose I will consider log(TPM + 0.1 pseudocount) for top 5000 expressed genes.

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

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

# 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.7765719

Now perform PCA, plot a few visualizations…

# pca with log-transformed count table (+ 0.1 pseudocount)
pca_results <- CountTable %>%
  +0.1 %>%
  mutate(sumVar = rowSums(.)) %>%
  arrange(desc(sumVar)) %>%
  head(2500) %>%
  select(-sumVar) %>%
  log() %>%
  t() %>%
  prcomp(center=T, scale. = T)

summary(pca_results)

Importance of components:
                           PC1      PC2      PC3     PC4     PC5     PC6
Standard deviation     39.2240 14.91153 10.40179 8.84746 8.19430 7.48033
Proportion of Variance  0.6154  0.08894  0.04328 0.03131 0.02686 0.02238
Cumulative Proportion   0.6154  0.70435  0.74763 0.77894 0.80580 0.82818
                           PC7     PC8     PC9    PC10    PC11    PC12
Standard deviation     6.53488 6.34844 6.20140 5.37654 4.99539 4.78658
Proportion of Variance 0.01708 0.01612 0.01538 0.01156 0.00998 0.00916
Cumulative Proportion  0.84526 0.86138 0.87677 0.88833 0.89831 0.90747
                          PC13    PC14    PC15    PC16    PC17    PC18
Standard deviation     4.70654 4.46198 4.16461 3.87579 3.72778 3.58367
Proportion of Variance 0.00886 0.00796 0.00694 0.00601 0.00556 0.00514
Cumulative Proportion  0.91634 0.92430 0.93124 0.93725 0.94280 0.94794
                          PC19    PC20    PC21   PC22    PC23   PC24
Standard deviation     3.44880 3.37702 3.21974 3.0002 2.95093 2.7392
Proportion of Variance 0.00476 0.00456 0.00415 0.0036 0.00348 0.0030
Cumulative Proportion  0.95270 0.95726 0.96141 0.9650 0.96849 0.9715
                          PC25    PC26    PC27    PC28    PC29    PC30
Standard deviation     2.64912 2.57087 2.51016 2.46120 2.41857 2.38343
Proportion of Variance 0.00281 0.00264 0.00252 0.00242 0.00234 0.00227
Cumulative Proportion  0.97430 0.97694 0.97946 0.98189 0.98423 0.98650
                          PC31    PC32    PC33    PC34    PC35    PC36
Standard deviation     2.30791 2.20390 2.15539 2.09649 2.07821 1.92604
Proportion of Variance 0.00213 0.00194 0.00186 0.00176 0.00173 0.00148
Cumulative Proportion  0.98863 0.99057 0.99243 0.99419 0.99592 0.99740
                          PC37    PC38      PC39
Standard deviation     1.83670 1.76876 2.137e-14
Proportion of Variance 0.00135 0.00125 0.000e+00
Cumulative Proportion  0.99875 1.00000 1.000e+00

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

Past versions of pca-1.png

ggplot(Merged, aes(x=PC1, y=PC2, color=factor(RNA.Library.prep.batch), label=Row.names)) + 
  geom_point() +
  geom_text_repel(size=2.5)

Past versions of pca-2.png

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

Past versions of pca-3.png

Can already see maybe something with the RNA extraction batch in first PC… Now I am going to look more systematically for significant correlations between potential observed confounders in the Metadata and the first 10 principle components. Will use Pearsons’s correlation to test continuous continuous confounders, will use anova for categorical confounders.

First the continuous 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", "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))

And the categorical confounders…

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

Plot some of the significant PC-metadata associations for a better look…

ggplot(Merged, aes(x=factor(RNA.Extract_date), y=PC6, label=Label)) + 
  geom_boxplot() +
  geom_jitter(position=position_jitter(width=.2, height=0))

ggplot(Merged, aes(x=factor(RNA.Library.prep.batch), y=PC6)) + 
  geom_boxplot() +
  geom_jitter(position=position_jitter(width=.1, height=0))

ggplot(Merged, aes(factor(RNA.Extract_date))) + 
  geom_bar(aes(fill = factor(RNA.Library.prep.batch)))

Strongest effect seems to be related to batch… RNA library prep batch and RNA extraction batch (date) both covary with a top PC. Have to check with Claudia that the last batch was a different extraction method (Trizol vs RNEasy).

ggplot(Merged, aes(x=factor(Source), y=PC3)) + 
  geom_boxplot() +
  geom_jitter(position=position_jitter(width=.1, height=0))

ggplot(Merged, aes(x=RIN, y=PC3, label=Label, color=factor(Source))) +
  geom_point() +
  geom_text_repel(size=2.5)

Also, Source covaries with PC3, though we don’t have many observations for a lot of Source categories so it may not make be a good idea to directly model source.

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] gplots_3.0.1    reshape2_1.4.3  knitr_1.22      ggrepel_0.8.0  
 [5] psych_1.8.10    forcats_0.4.0   stringr_1.4.0   dplyr_0.8.0.1  
 [9] purrr_0.3.2     readr_1.3.1     tidyr_0.8.2     tibble_2.1.1   
[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] 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       whisker_0.3-2      rmarkdown_1.11    
[22] labeling_0.3       foreign_0.8-71     munsell_0.5.0     
[25] broom_0.5.1        compiler_3.5.1     modelr_0.1.4      
[28] xfun_0.6           pkgconfig_2.0.2    mnormt_1.5-5      
[31] htmltools_0.3.6    tidyselect_0.2.5   gridExtra_2.3     
[34] workflowr_1.2.0    crayon_1.3.4       withr_2.1.2       
[37] bitops_1.0-6       grid_3.5.1         nlme_3.1-137      
[40] jsonlite_1.6       gtable_0.3.0       git2r_0.24.0      
[43] magrittr_1.5       scales_1.0.0       KernSmooth_2.23-15
[46] cli_1.1.0          stringi_1.4.3      fs_1.2.6          
[49] xml2_1.2.0         generics_0.0.2     tools_3.5.1       
[52] glue_1.3.1         hms_0.4.2          parallel_3.5.1    
[55] yaml_2.2.0         colorspace_1.4-1   caTools_1.17.1.1  
[58] rvest_0.3.2        haven_2.1.0