20190429_RNASeqSTAR_quantifications

workflowr

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

Untracked files:
    Untracked:  analysis/20190428_Check_eQTLs.Rmd
    Untracked:  analysis/20190429_RNASeqSTAR_quantifications.Rmd
    Untracked:  data/PastAnalysesDataToKeep/20190428_log10TPM.txt.gz
    Untracked:  docs/figure/20190429_RNASeqSTAR_quantifications.Rmd/

Unstaged changes:
    Modified:   analysis/20190424_Check_eQTLs.Rmd
    Modified:   analysis/20190426_CheckRNASeqPCs_2.Rmd

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

All my previous analyses (for example, here) were based on RNA-seq quantifications using the kallisto mapper/quantifier. I have noticed at low expression levels especially, there is a tendency for large variation due to a single outlier sample. I was wondering if this is a result of the kallisto softare quantification. So I remapped all the RNA-seq data with STAR aligner and here I will compare, looking specifically for outliers.

# Read in count table from kallisto pseudoalignment/quantification (not raw counts but transformed to log10 TPM)
CountTableKallisto <- read.table(gzfile('../data/PastAnalysesDataToKeep/20190428_log10TPM.txt.gz'), header=T, check.names=FALSE, row.names = 1)
kable(CountTableKallisto[1:10,1:10])

# Read in count table from STAR alignments
CountTableSTAR <- read.table(gzfile('../data/PastAnalysesDataToKeep/20190429_STAR.CountTable.txt.gz'), header=T, check.names=FALSE, row.names = 1)
kable(CountTableSTAR[1:10,1:10])

# Some initial checks of STAR  count table before doing any transformations
dim(CountTableSTAR)

[1] 31373    39

# Millions of mapped reads per sample
kable(colSums(CountTableSTAR)/1000000)

# XY scatter of first two samples
qplot(log10(CountTableSTAR$`338`), log10(CountTableSTAR$`295`))

#Histogram of average gene expressions measured in counts.
qplot(log10(rowMeans(CountTableSTAR)))

`stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

Warning: Removed 4416 rows containing non-finite values (stat_bin).

# Cumulative plot of mean reads per gene. Added psuedo count
plot(ecdf(log10(rowMeans(CountTableSTAR) + 0.001)), xlab="log10Counts", ylab="Fraction of genes with at less than x counts on average")

# Now test some filtering methods, mostly based on minimum read counts:

# 80% of samples must have at least 10 reads
Q <- rowQuantiles(as.matrix(CountTableSTAR), probs=0.2)
# Now some filtering:
GeneSetFilter1 <- names(which(Q>10))

# 100% of samples must have at least 1 read
GeneSetFilter2<-CountTableSTAR %>%
  rownames_to_column('gene') %>%
  filter_if(is.numeric, all_vars(.>0)) %>%
  pull(gene)

# 100% of samples must have >0 TPM (quantified by kallisto)
GeneSetFilter3<-rownames(CountTableKallisto)

length(GeneSetFilter1)

[1] 14059

length(GeneSetFilter2)

[1] 14755

length(GeneSetFilter3)

[1] 14824

length(intersect(GeneSetFilter1, GeneSetFilter2))

[1] 13967

length(intersect(GeneSetFilter2, GeneSetFilter3))

[1] 13574

length(intersect(GeneSetFilter1, GeneSetFilter3))

[1] 13119

length(intersect(intersect(GeneSetFilter1, GeneSetFilter2), GeneSetFilter3))

[1] 13094

All of these filter methods result in highly overlapping sets of ~14000 genes. Conceptually I prefer method 1, where 80% of samples must contain at least 10 reads, thus allowing some samples to have much less or even 0. This method will require a pseudocount to deal with log-transform.

Now do some basic transformations. I will just convert to log10(CPM).

CountTableFiltered <- CountTableSTAR %>% as.data.frame() %>%
  rownames_to_column('gene') %>%
  filter(gene %in% GeneSetFilter1) %>%
  column_to_rownames('gene')

# table of log10COM
log10CPM_table <- log10((CountTableFiltered + 1) / colSums(CountTableFiltered))
kable(log10CPM_table[1:10,1:10])

Now make correlation matrix from STAR quantifications

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

CorMatrix <- log10CPM_table %>%
  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)

# A couple of pairwise comparisons
qplot(log10CPM_table$`4X0095`, log10CPM_table$`4X0212`)

qplot(log10CPM_table$`295`, log10CPM_table$Little_R)

Now make the same plots from the kallisto quantifications

CorMatrix <- CountTableKallisto %>%
  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)

qplot(CountTableKallisto$`4X0095`, CountTableKallisto$`4X0212`)

qplot(CountTableKallisto$`295`, CountTableKallisto$Little_R)

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