Dataset selection and quality control

In order to avoid over-fitting of signatures derived from human expression data, we decided to use a normal tissue expression dataset from a non-human mammal species, because it has been shown that tissue-preferential expression patterns tend to be conserved between mammal species. We identified a dataset of Canis lupus familiaris (dog), which is publicly available in Gene Expression Omnibus (GDS4164).

In this study, the authors examined 39 samples from 10 pathologically normal tissues (liver, kidney, heart, lung, brain, lymph node, spleen, jejunum, pancreas, and skeletal muscle) of four dogs (with one pancreas sample missing). We downloaded the data, and performed minimal pre-processing: for multiple probesets that map to same genes, we kept the one with the highest average expression level and removed the rest. The resulting dataset contained expression of 16797 genes. BioQC was applied to the dataset to test whether there are major contamination issues. The results, including tissues reported by the authors, and the BioQC tissue signatures with the highest and second-highest scores, are reported in the following table:

By comparing the tissue labels provided by the authors and the predictions of BioQC, we notice that in most cases the two match well (despite of ontological differences). In three cases (sample ID GSM502573, GSM502594, and GSM502596) though, there seem to be intriguing differences, which might be explained by different sampling procedures or immune cell infiltration. We will however in this vignette not further explore them. These three samples are removed from the simulation procedures.

An example of weighted mixing: heart and jejunum

As an example, we take average expression of heart and jejunum samples, and mix them by different compositions. This allows us comparing enrichment scores and their ranks when the expression profiles of heart and jejunum are mixed in silico:

Figure 2: Results of a mixing case study. Left panel: BioQC enrichment scores of small intestine and cardiac muscle varying upon different proportions of jejunum; Right panel: ranks of enrichment scores varying upon different proportions of jejunum.

We observe that with as little as 5% contamination of heart tissue in jejunum samples (rightmost in the right panel), the rank of heart signature jumps from 34 to 9; 10% and 20% contamination will further enhance the rank to 4 and 3 respectively. If we start from the other end, namely assuming jejunum contamination in heart samples, the BioQC algorithms ranks jejunum the 7th only when there are more than 25% contamination. If we set enrichment score equal or over 3 as the threshold of calling a suspected contamination event (\(p<0.001\) in the one-sided Wilcoxon-Mann-Whitney test), it takes about 10% heart in jejunum tissue or about 30% jejunum tissue in heart to make a call. It means the sensitivity of contamination detection is not symmetric between tissues: contamination by tissues with distinct expression patterns (such as heart) are easier to be detected than contamination by tissues with less distinct expression patterns (such as small intestine).

While it is difficult to quantify the absolute sensitivity of contamination detection, it is apparent that if the enrichment score of a unforeseen tissue is very high (or ranked high), one may suspect potential contamination. Also, if there are replicates of samples from the same tissue, a higher value in one sample compared with the other samples suggests a contamination or infiltration incident.

Pairwise Mixing

Following the heart-jejunum example, we performed all 45 pairwise mixing experiments, producing weighted linear combinations of gene expression profiles of each pair of tissues (excluding self-mixing). The results are summaried in a heatmap:

Figure 3: Results of the pairwise mixing experiment. Each cell represents the minimal percentage of tissue of the column as contamination in the tissue of the row that can be detected by BioQC. No values are available for cells on the diagonal because self-mixing was excluded. Heart and skeletal muscle are very close to each other and therefore their detection limit is not considered.

The heatmap visualization summarizes the detection limit of contamination of each pair of tissues. Take the cell in row 1 column 2 from top left: its value (0.15) means that if there are 15% or more contamination by heart in the brain sample, BioQC will be able to detect it (with the threshold enrichment score \(\geqslant3\) or the rank \(\leqslant10\)), because the enrichment score is equal to or larger than 3, or the heart tissue signature ranks in the top 10 of all tissue signatures.

Take another cell in row 2 column 1 from top left: its value (0.5) means that if there are 50% or more contanmination by brain in a heart sample, BioQC will be able to detect it. Here we observe the asymmetry again that we observed before with the heart/jejenum example: while it is realtively easy to identify heart contamination of a brain sample, it is more difficult to identify brain contamination of a heart sample in this dataset.

The average detection limits of tissues as contamination sources are listed in the following table. The values are derived from median values of each column in the heatmap, except for diagonal and missing elements.

Interestingly, brain samples are a special case: if they contaminate other tissues, it is more difficult to identify (but not other way around). It can be partially explained by the experiment design: Briggs et al. profiled the whole cerebrum, whereas in BioQC there are 22 distinct gene sets assigned to distinct brain regions. Though the prefrontal cortex signature scored highest in the canine brain samples, its score is relative low (7.45), and the genes in the signature are not too far away from the background genes:

Figure 4: Average expression of tissue-preferential genes in respective tissues. For each tissue (e.g. brain), we calculate the median ratio of gene expression level of specific genes over the median expression level of background genes. The value reflects the specificity of tissue-specific genes in respective tissues. Likely due to the sampling of different brain regions, the score of brain ranks the lowest.

Therefore only a strong contamination by brain in this dataset will be detected by the given threshold. We expect that if prefrontal cortex instead of cerebrum sample was profiled, the mixing profile of brain will be similar to other organs. This needs to be tested in other datasets.

Apart from that, most in silico contamination events are detectable in this dataset, with median detection limit around 0.2. This suggests that BioQC is sensitive towards tissue heterogeneity in physiological settings.

Comparing BioQC with Principal Component Analysis (PCA)

In the context of the dog transcriptome dataset, we can compare the results of principal component analysis (PCA) with that of BioQC:

Sample separation revealed by principal component analysis (PCA) of the dog transcriptome dataset.

PCA sugggests that samples from each tissue tend to cluster together, in line with the BioQC results. In addition, PCA reveals that tissues with cells of similar origins cluster together, such as skeletal muscle and heart. As expected, one brain sample and two lung samples seem to be different from other samples of the same cluster, which are consistent with the BioQC findings; however, unlike BioQC, PCA does not provide information on what are potential contamination/infiltration casues.

Therefore, we believe BioQC complements existing unsupervised methods to inspect quality of gene expression data.

R Session Info

sessionInfo()

## R version 4.4.1 (2024-06-14)
## Platform: x86_64-pc-linux-gnu
## Running under: Ubuntu 24.04.1 LTS
## 
## Matrix products: default
## BLAS:   /home/biocbuild/bbs-3.20-bioc/R/lib/libRblas.so 
## LAPACK: /usr/lib/x86_64-linux-gnu/lapack/liblapack.so.3.12.0
## 
## 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       
## 
## time zone: America/New_York
## tzcode source: system (glibc)
## 
## attached base packages:
## [1] stats4    stats     graphics  grDevices utils     datasets  methods  
## [8] base     
## 
## other attached packages:
##  [1] hgu133plus2.db_3.13.0 rbenchmark_1.0.0      gplots_3.2.0         
##  [4] gridExtra_2.3         latticeExtra_0.6-30   lattice_0.22-6       
##  [7] org.Hs.eg.db_3.20.0   AnnotationDbi_1.68.0  IRanges_2.40.0       
## [10] S4Vectors_0.44.0      testthat_3.2.1.1      limma_3.62.0         
## [13] RColorBrewer_1.1-3    BioQC_1.34.0          Biobase_2.66.0       
## [16] BiocGenerics_0.52.0   knitr_1.48           
## 
## loaded via a namespace (and not attached):
##  [1] KEGGREST_1.46.0         gtable_0.3.6            xfun_0.48              
##  [4] bslib_0.8.0             caTools_1.18.3          vctrs_0.6.5            
##  [7] tools_4.4.1             bitops_1.0-9            RSQLite_2.3.7          
## [10] highr_0.11              blob_1.2.4              pkgconfig_2.0.3        
## [13] KernSmooth_2.23-24      desc_1.4.3              lifecycle_1.0.4        
## [16] GenomeInfoDbData_1.2.13 compiler_4.4.1          deldir_2.0-4           
## [19] Biostrings_2.74.0       brio_1.1.5              statmod_1.5.0          
## [22] GenomeInfoDb_1.42.0     htmltools_0.5.8.1       sass_0.4.9             
## [25] yaml_2.3.10             crayon_1.5.3            jquerylib_0.1.4        
## [28] cachem_1.1.0            gtools_3.9.5            locfit_1.5-9.10        
## [31] digest_0.6.37           rprojroot_2.0.4         fastmap_1.2.0          
## [34] grid_4.4.1              cli_3.6.3               magrittr_2.0.3         
## [37] edgeR_4.4.0             UCSC.utils_1.2.0        bit64_4.5.2            
## [40] rmarkdown_2.28          XVector_0.46.0          httr_1.4.7             
## [43] jpeg_0.1-10             bit_4.5.0               interp_1.1-6           
## [46] png_0.1-8               memoise_2.0.1           evaluate_1.0.1         
## [49] rlang_1.1.4             Rcpp_1.0.13             glue_1.8.0             
## [52] DBI_1.2.3               pkgload_1.4.0           jsonlite_1.8.9         
## [55] R6_2.5.1                zlibbioc_1.52.0