Introduction to zFPKM Transformation

Ron Ammar

2021-10-26

## Registered S3 method overwritten by 'printr':
##   method                from     
##   knit_print.data.frame rmarkdown

Summary

Perform the zFPKM transform on RNA-seq FPKM data. This algorithm is based on the publication by Hart et al., 2013 (Pubmed ID 24215113). The reference recommends using zFPKM > -3 to select expressed genes. Validated with ENCODE open/closed promoter chromatin structure epigenetic data on six of the ENCODE cell lines. It works well for gene level data using FPKM or TPM, but does not appear to calibrate well for transcript level data.

Identifying active genes for subsequent analysis

We calculate zFPKM for existing FPKM from gse94802.

library(dplyr)
library(GEOquery)
library(stringr)
library(SummarizedExperiment)
library(tidyr)

getSpecificGEOSupp <- function(url) {
  temp <- tempfile()
  download.file(url, temp)
  out <- read.csv(gzfile(temp), row.names=1, check.names=FALSE)
  out <- select(out, -MGI_Symbol)
  return(out)
}

gse94802_fpkm <- "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE94nnn/GSE94802/suppl/GSE94802_Minkina_etal_normalized_FPKM.csv.gz"
gse94802_counts <- "ftp://ftp.ncbi.nlm.nih.gov/geo/series/GSE94nnn/GSE94802/suppl/GSE94802_Minkina_etal_raw_counts.csv.gz"

if (file.exists("gse94802.rds")) {
  esetlist <- readRDS("gse94802.rds")
} else {
  esetlist <- getGEO("gse94802")
}

doe <- pData(esetlist[[1]])

colData <- DataFrame(
  condition=ifelse(str_detect(doe$title, regex("control", ignore_case=TRUE)), "control", "mutant"),
  sample_id=str_match(doe$title, "rep\\d_(.+)")[, 2],
  row.names=str_match(doe$title, "rep\\d_(.+)")[, 2])

se <- SummarizedExperiment(assays=SimpleList(fpkm=getSpecificGEOSupp(gse94802_fpkm),
                                             counts=getSpecificGEOSupp(gse94802_counts)),
                           colData=colData)

# clear namespace
rm(esetlist, gse94802_fpkm, gse94802_counts, doe, colData, getSpecificGEOSupp)

We compute zFPKM.

library(zFPKM)
assay(se, "zfpkm") <- zFPKM(se)

We can also plot the Guassian fit to the FPKM data for which the z-scores are based.

zFPKMPlot(se)

To determine which genes are active, we compute the median expression within each group.

activeGenes <- assay(se, "zfpkm") %>%
  mutate(gene=rownames(assay(se, "zfpkm"))) %>%
  gather(sample_id, zfpkm, -gene) %>%
  left_join(select(as.data.frame(colData(se)), sample_id, condition), by="sample_id") %>%
  group_by(gene, condition) %>%
  summarize(median_zfpkm=median(zfpkm)) %>%
  ungroup() %>%
  mutate(active=(median_zfpkm > -3)) %>%
  filter(active) %>%
  select(gene) %>%
  distinct()

seActive <- SummarizedExperiment(
  assays=SimpleList(counts=as.matrix(assay(se, "counts")[activeGenes$gene, ])),
  colData=colData(se))

In the following DE analysis, we only use genes that were active in either group.

library(limma)
library(edgeR)

# Generate normalized log2CPM from counts AFTER we filter for protein-coding
# genes that are detectably expressed.
dge <- DGEList(counts=assay(seActive, "counts"))
dge <- calcNormFactors(dge)
design <- model.matrix(~ 0 + condition, data=colData(seActive))
vq <- voomWithQualityWeights(dge, design, plot=TRUE)

fit <- lmFit(vq, design)
contrastMatrix <- makeContrasts(conditioncontrol - conditionmutant, levels=design)
fit <- contrasts.fit(fit, contrastMatrix)
fit <- eBayes(fit, robust=TRUE)
deGenes <- topTable(fit, number=Inf)

References

Hart T, Komori HK, LaMere S, Podshivalova K, Salomon DR. Finding the active genes in deep RNA-seq gene expression studies. BMC Genomics. 2013 Nov 11;14:778. doi: 10.1186/1471-2164-14-778.

## R version 4.1.1 (2021-08-10)
## Platform: x86_64-pc-linux-gnu (64-bit)
## Running under: Ubuntu 20.04.3 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.14-bioc/R/lib/libRblas.so
## LAPACK: /home/biocbuild/bbs-3.14-bioc/R/lib/libRlapack.so
## 
## locale:
##  [1] LC_CTYPE=en_US.UTF-8       LC_NUMERIC=C              
##  [3] LC_TIME=en_GB              LC_COLLATE=C              
##  [5] LC_MONETARY=en_US.UTF-8    LC_MESSAGES=en_US.UTF-8   
##  [7] LC_PAPER=en_US.UTF-8       LC_NAME=C                 
##  [9] LC_ADDRESS=C               LC_TELEPHONE=C            
## [11] LC_MEASUREMENT=en_US.UTF-8 LC_IDENTIFICATION=C       
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] edgeR_3.36.0                limma_3.50.0               
##  [3] zFPKM_1.16.0                tidyr_1.1.4                
##  [5] SummarizedExperiment_1.24.0 GenomicRanges_1.46.0       
##  [7] GenomeInfoDb_1.30.0         IRanges_2.28.0             
##  [9] S4Vectors_0.32.0            MatrixGenerics_1.6.0       
## [11] matrixStats_0.61.0          stringr_1.4.0              
## [13] GEOquery_2.62.0             Biobase_2.54.0             
## [15] BiocGenerics_0.40.0         dplyr_1.0.7                
## [17] printr_0.2                 
## 
## loaded via a namespace (and not attached):
##  [1] Rcpp_1.0.7             locfit_1.5-9.4         lattice_0.20-45       
##  [4] assertthat_0.2.1       digest_0.6.28          utf8_1.2.2            
##  [7] R6_2.5.1               backports_1.2.1        evaluate_0.14         
## [10] ggplot2_3.3.5          highr_0.9              pillar_1.6.4          
## [13] zlibbioc_1.40.0        rlang_0.4.12           data.table_1.14.2     
## [16] jquerylib_0.1.4        Matrix_1.3-4           checkmate_2.0.0       
## [19] rmarkdown_2.11         labeling_0.4.2         statmod_1.4.36        
## [22] readr_2.0.2            RCurl_1.98-1.5         munsell_0.5.0         
## [25] DelayedArray_0.20.0    compiler_4.1.1         xfun_0.27             
## [28] pkgconfig_2.0.3        htmltools_0.5.2        tidyselect_1.1.1      
## [31] tibble_3.1.5           GenomeInfoDbData_1.2.7 fansi_0.5.0           
## [34] crayon_1.4.1           tzdb_0.1.2             bitops_1.0-7          
## [37] grid_4.1.1             jsonlite_1.7.2         gtable_0.3.0          
## [40] lifecycle_1.0.1        DBI_1.1.1              magrittr_2.0.1        
## [43] scales_1.1.1           stringi_1.7.5          farver_2.1.0          
## [46] XVector_0.34.0         xml2_1.3.2             bslib_0.3.1           
## [49] ellipsis_0.3.2         generics_0.1.1         vctrs_0.3.8           
## [52] tools_4.1.1            glue_1.4.2             purrr_0.3.4           
## [55] hms_1.1.1              fastmap_1.1.0          yaml_2.2.1            
## [58] colorspace_2.0-2       knitr_1.36             sass_0.4.0