2.1 Getting started using dear_seq()

2 inputs are required to run dear_seq():

• The gene expression matrix
• The data.frame containing the variables of interested to be tested (e.g.  the group factor)

Both can be supplied in the formed of a SummarizedExperiment object for instance. A third optional input is the covariates that will be adjusted on (such as age for instance).

2.2 An example analysis

Tuberculosis (TB) is a disease caused by a bacterium called Mycobacterium tuberculosis. Bacteria typically infect the lungs, but they can also affect other parts of the body. Tuberculosis can remain in a quiescent state in the body. It is called latent tuberculosis infection (LTBI). It is an infection without clinical signs, bacteriological and radiological disease.

Berry et al. first identified a whole blood 393 transcript signature for active TB using microarray analysis. Then, Singhania et al. found a 373-genes signature of active tuberculosis using RNA-Seq, confirming microarray results, that discriminates active tuberculosis from latently infected and healthy individuals.

We sought to investigate how many of the 373 genes Singhania et al. found using edgeR might actually be false positives. As a quick example, we propose a partial re-analysis (see Gauthier et al. for the complete analysis) of the Singhania et al. dataset: we proceed to the Differential Expression Analysis (DEA) of the Active TB group against the Control one, omitting LTBI/Control and ActiveTB/LTBI comparisons. We focused on the Berry London dataset.

It results in 54 patients of whom 21 are active TB patients, 21 are latent TB patients and 12 are healthy control particiants.

For more details, check the article from Berry et al., 2010 here.

library(GEOquery)
GSE107991_metadata <- GEOquery::getGEO("GSE107991", GSEMatrix = FALSE)

get_info <- function(i) {
    name <- GEOquery::Meta(GSMList(GSE107991_metadata)[[i]])$source_name_ch1
    name <- gsub("Active_TB", "ActiveTB", name)
    name <- gsub("Test_set", "TestSet", name)
    c(unlist(strsplit(name, split = "_")), GEOquery::Meta(GSMList(GSE107991_metadata)[[i]])$title)
}

infos <- vapply(X = 1:length(GSMList(GSE107991_metadata)), FUN = get_info, FUN.VALUE = c("a", 
    "b", "c", "d", "e"))
infos_df <- cbind.data.frame(SampleID = names(GSMList(GSE107991_metadata)), t(infos), 
    stringsAsFactors = TRUE)
rownames(infos_df) <- names(GEOquery::GSMList(GSE107991_metadata))
colnames(infos_df)[-1] <- c("Cohort", "Location", "Set", "Status", "SampleTitle")

library(readxl)
GEOquery::getGEOSuppFiles("GSE107991", filter_regex = "edgeR_normalized_Berry_London")
London <- readxl::read_excel("GSE107991/GSE107991_edgeR_normalized_Berry_London.xlsx")
genes <- London$Genes
London <- as.matrix(London[, -c(1:3)])
rownames(London) <- genes

library(dearseq)
library(SummarizedExperiment)
col_data <- data.frame(Status = infos_df$Status)
col_data$Status <- stats::relevel(col_data$Status, ref = "Control")
rownames(col_data) <- infos_df$SampleTitle
se <- SummarizedExperiment(assays = London, colData = col_data)

res_dearseq <- dearseq::dear_seq(object = se[, se$Status != "LTBI"], variables2test = "Status", 
    which_test = "asymptotic", preprocessed = TRUE)

5.1 GEOquery package

In case the data you want to analyze is publicly available through Gene Expression Omnibus (GEO), you can access it with the GEOquery package, that can be installed with the following commands:

if (!requireNamespace("GEOquery", quietly = TRUE)) {
    if (!requireNamespace("BiocManager", quietly = TRUE)) {
        install.packages("BiocManager")
    }
    BiocManager::install("GEOquery")
}

5.2 readxl

The readxl package allows to easily imports data from Excel format and can be similarly installed with the following commands:

if (!requireNamespace("readxl", quietly = TRUE)) {
    install.packages("readxl")
}

More details can be found on Bioconductor and in Davis S, Meltzer P, (2007) GEOquery: a bridge between the Gene Expression Omnibus (GEO) and Bioconductor Bioinformatics 14:1846-1847.

5.3 Session Info

#> R version 4.1.0 RC (2021-05-10 r80283)
#> Platform: x86_64-w64-mingw32/x64 (64-bit)
#> Running under: Windows Server x64 (build 17763)
#> 
#> Matrix products: default
#> 
#> locale:
#> [1] LC_COLLATE=C                          
#> [2] LC_CTYPE=English_United States.1252   
#> [3] LC_MONETARY=English_United States.1252
#> [4] LC_NUMERIC=C                          
#> [5] LC_TIME=English_United States.1252    
#> 
#> attached base packages:
#> [1] stats4    parallel  stats     graphics  grDevices utils     datasets 
#> [8] methods   base     
#> 
#> other attached packages:
#>  [1] BiocSet_1.6.0               dplyr_1.0.6                
#>  [3] SummarizedExperiment_1.22.0 GenomicRanges_1.44.0       
#>  [5] GenomeInfoDb_1.28.0         IRanges_2.26.0             
#>  [7] S4Vectors_0.30.0            MatrixGenerics_1.4.0       
#>  [9] matrixStats_0.58.0          dearseq_1.4.0              
#> [11] readxl_1.3.1                GEOquery_2.60.0            
#> [13] Biobase_2.52.0              BiocGenerics_0.38.0        
#> 
#> loaded via a namespace (and not attached):
#>  [1] httr_1.4.2             sass_0.4.0             tidyr_1.1.3           
#>  [4] bit64_4.0.5            jsonlite_1.7.2         bslib_0.2.5.1         
#>  [7] assertthat_0.2.1       statmod_1.4.36         blob_1.2.1            
#> [10] GenomeInfoDbData_1.2.6 cellranger_1.1.0       yaml_2.2.1            
#> [13] pillar_1.6.1           RSQLite_2.2.7          lattice_0.20-44       
#> [16] glue_1.4.2             limma_3.48.0           digest_0.6.27         
#> [19] XVector_0.32.0         colorspace_2.0-1       plyr_1.8.6            
#> [22] htmltools_0.5.1.1      Matrix_1.3-3           pkgconfig_2.0.3       
#> [25] zlibbioc_1.38.0        purrr_0.3.4            scales_1.1.1          
#> [28] tibble_3.1.2           KEGGREST_1.32.0        generics_0.1.0        
#> [31] ggplot2_3.3.3          ellipsis_0.3.2         cachem_1.0.5          
#> [34] pbapply_1.4-3          magrittr_2.0.1         crayon_1.4.1          
#> [37] memoise_2.0.0          evaluate_0.14          fansi_0.4.2           
#> [40] xml2_1.3.2             tools_4.1.0            hms_1.1.0             
#> [43] BiocIO_1.2.0           formatR_1.9            lifecycle_1.0.0       
#> [46] stringr_1.4.0          munsell_0.5.0          DelayedArray_0.18.0   
#> [49] Biostrings_2.60.0      AnnotationDbi_1.54.0   compiler_4.1.0        
#> [52] jquerylib_0.1.4        rlang_0.4.11           grid_4.1.0            
#> [55] RCurl_1.98-1.3         CompQuadForm_1.4.3     bitops_1.0-7          
#> [58] rmarkdown_2.8          gtable_0.3.0           DBI_1.1.1             
#> [61] curl_4.3.1             R6_2.5.0               knitr_1.33            
#> [64] fastmap_1.1.0          bit_4.0.4              utf8_1.2.1            
#> [67] ontologyIndex_2.7      KernSmooth_2.23-20     readr_1.4.0           
#> [70] stringi_1.6.2          Rcpp_1.0.6             png_0.1-7             
#> [73] vctrs_0.3.8            tidyselect_1.1.1       xfun_0.23