Eye in a Disk: eyeIntegration Human Pan-Eye and Body Transcriptome Database Version 1.0

Abstract

Purpose: We develop an accessible and reliable RNA sequencing (RNA-seq) transcriptome database of healthy human eye tissues and a matching reactive web application to query gene expression in eye and body tissues.

Methods: We downloaded the raw sequence data for 1375 RNA-seq samples across 54 tissues in the Genotype-Tissue Expression (GTEx) project as a noneye reference set. We then queried several public repositories to find all healthy, nonperturbed, human eye-related tissue RNA-seq samples. The 916 eye and 1375 GTEx samples were sent into a Snakemake-based reproducible pipeline we wrote to quantify all known transcripts and genes, removes samples with poor sequence quality and mislabels, normalizes expression values across each tissue, perform 882 differential expression tests, calculate GO term enrichment, and output all as a single SQLite database file: the Eye in a Disk (EiaD) dataset. Furthermore, we rewrote the web application eyeIntegration (available in the public domain at https://eyeIntegration.nei.nih.gov) to display EiaD.

Results: The new eyeIntegration portal provides quick visualization of human eye-related transcriptomes published to date by database version, gene/transcript, 19 eye tissues, and 54 body tissues. As a test of the value of this unified pan-eye dataset, we showed that fetal and organoid retina are highly similar at a pan-transcriptome level, but display distinct differences in certain pathways and gene families, such as protocadherin and HOXB family members.

Conclusions: The eyeIntegration v1.0 web app serves the pan-human eye and body transcriptome dataset, EiaD. This offers the eye community a powerful and quick means to test hypotheses on human gene and transcript expression across 54 body and 19 eye tissues.

From anterior to posterior along the light trajectory, the human eye is composed of the cornea, lens, retina, RPE, and choroid. The differentiation, maturation, and function of these tissues is mediated through spatial- and temporal-specific transcript and gene expression patterns, also known as the transcriptome. Today, RNA-sequencing (RNA-seq) is the predominant technology for quantifying the transcriptome. Analysis of the transcripts' expression across tissue, time, and perturbation allows researchers to decipher the genetic controls of eye development and function. To this end, a wide variety of human tissue sources have been used to assess gene function, including primary tissue (fetal and postmortem), differentiated stem cells, immortalized cell lines, and most recently, organoids. These tissue types have been deeply sequenced across the cornea,^1–7 lens,⁸ retina,^9–17 and RPE (choroid).^14,17–34 

The Genotype-Tissue Expression (GTEx) Gene Expression Web App Lacks Eye-Specific Tissues

The GTEx Project has generated RNA-seq data across dozens of postmortem human tissues from hundreds of unique donors, and presents the gene and transcript level data in a comprehensive and user-friendly web app (available in the public domain at https://gtexportal.org/); however eye tissues have not been included.^35,36 Recently Ratnapriya et al.³⁷ reported on a huge set of postmortem retina, normal and with varying degrees of AMD and the GTEx project is providing the data as a download link. These data, as of June 2019, are not available in the interactive GTEx visualizations.³⁷ The Sequence Read Archive (SRA) and European Nucleotide Archive (ENA) are the primary repositories for all raw sequence data and two groups have quantified large portions of the RNA-seq data, including some human eye tissues, from the SRA: recount2 and ARCHS4.^38,39 To date, no curation of the sample level metadata has been done; therefore, it is challenging to parse out which eye tissues are present and even more difficult to determine whether any samples were chemically or genetically perturbed. More targeted web resources that allow researchers to quickly assess gene expression in eye tissues include iSYTE, EXPRESS, and retina.Tigem.it.^16,40,41 However iSYTE only includes lens samples, EXPRESS is limited to a subset of mouse lens and retina samples, and retina.Tigem.it is retina only. Thus, we aimed our efforts at developing an easily accessible and reliable RNA-seq based transcriptome database of healthy human eye tissues and a matching reactive web application to query gene expression in eye and body tissues. 

The eyeIntegration App Interactively Serves Huge GTEx and Human Eye Tissue Datasets (Eye in a Disk [EiaD])

The eyeIntegration web resource (available in the public domain at https://eyeIntegration.nei.nih.gov), originally released in 2017 at version 0.6, provides the largest set of transcriptomes from hand-curated human eye tissues along with hundreds of GTEx tissue samples.⁴² This interactive web app allows for quick transcript and gene comparisons across many eye tissues and dozens of other body tissues. The dataset that the original eyeIntegration web app served was created with a series of scripts, several of which were run interactively to manually assess quality control for the samples. The interactive nature of some of the steps precluded efficient and regular data updates for the data. 

To better meet the needs of the eye research community we have rewritten the bioinformatic pipeline that creates the eye and body RNA-seq dataset to allow for regular, versioned updates for eyeIntegration. We call this reproducible and versioned transcriptome dataset “Eye in a Disk” (EiaD). The pipeline automates the EiaD creation, ensures full reproducibility of the results, allow for external data comparison, provides consistent sample quality control, and improves efficiency for future sample updates. The 2019 EiaD dataset contains several new tissue types, full gene product quantification, along with hundreds of new samples and improved sample labeling. The eyeIntegration web app also has been rewritten to provide many new features, including versioned EiaD datasets, custom URL shortcut creation, new visualizations, improved data table searching, easy download of core datasets, and local install of the entire interactive resource with three commands. Additionally, we are prototyping new tools to display single cell RNA-seq (scRNA-seq) data to provide researchers access to cell type–specific information about gene expression across murine retinal development. 

The EiaD Dataset Can Be Used to Identify Potential Avenues to Improve Retina Organoid Maturation

Retina organoids are an increasingly popular means to model human retina development. We used our pan-study EiaD dataset to show that, at a pan-transcriptome level, organoids are highly similar to early fetal retina tissue. We also showed that important temporal gene expression patterns in the fetal retina tissue are recapitulated in the organoids. As the organoid differentiation methods do not yet produce fully mature retina, we focused on identifying differentially expressed processes between organoid retina and embryonic retina, and detected, for example, identifying protocadherin and HOXB family gene expression differences that suggest targetable pathways to improve and benchmark organoid differentiation methods. 

Methods

Identification of Potential Eye Samples

We exhaustively searched the SRA with the SRAdb R package for eye related tissues using the query ‘cornea|retina|RPE|macula|fovea|choroid|sclera|iris|lens|eye' across all columns and rows in the ‘SRA' table.^43,44 As the SRAdb is being deprecated, we also ran searches on the SRA and Gene Expression Omnibus (GEO) web pages with as follows: ((“Homo sapiens”[orgn:__txid9606]) AND (transcriptomic[Source]) AND (“2019/01/01”[Publication Date] : “3000”[Publication Date]) AND (retina[Text Word] OR RPE[Text Word] OR macula[Text Word] OR fovea[Text Word] OR choroid[Text Word] OR sclera[Text Word] OR iris[Text Word] OR lens[Text Word] OR cornea[Text Word] OR ‘trabecular meshwork'[Text Word] OR ‘canals of schlemm'[Text Word] OR ‘cillary body'[Text Word] OR ‘optic nerve'[Text Word] OR ‘laminar cribosa'[Text Word] OR retina[Title] OR RPE[Title] OR macula[Title] OR fovea[Title] OR choroid[Title] OR sclera[Title] OR iris[Title] OR lens[Title] OR cornea[Title] OR ‘trabecular meshwork'[Title] OR ‘canals of schlemm'[Title] OR ‘cillary body'[Title] OR ‘optic nerve'[Title] OR ‘laminar cribosa'[Title]))”. We hand selected relevant studies and selected healthy, control or unmodified samples spanning primary adult tissue, primary fetal tissue, induced pluripotent stem cell (iPSC)–derived tissue, stem cell–derived organoids, and immortalized cell lines. To compare gene expression in the eye against expression in other body tissues, we obtained samples from 54 different body tissues from the GTEx project. Using SRA metadata from each study, we extracted sample and run accessions, library type, tissue of origin, and subtissue of origin. Any of the preceding information missing from the SRA metadata was added by hand, when available. Stem cell-derived tissues and cell lines are marked as subtissues of the tissue they model. 

Raw Data Download and Quantification

We downloaded the relevant SRA files for each sample directly from the NCBI ftp server using the file transfer software Aspera. SRA files were converted to FASTQ format using the tool fastq-dump from the SRAtoolkit software package.⁴³ Samples only available in the BAM format were converted to FASTQ format using SAMTools.⁴⁵ Sample transcriptomes were quantified using the alignment free quantification software Salmon, using transcriptomic index built from gencode v28 protein coding transcript sequences using the transcriptomic aligner Salmon.^46,47 Using the resulting expression quantification, we identified lowly or unused transcripts within the gencode annotation, and removed transcripts that accounted for 5% or less of the total expression for its parent gene as per Soneson et al.⁴⁸ Samples were requantified against a transcriptomic index built on the filtered transcript sequences. The Salmon count values were quantified as (transcript) length scaled Transcripts Per Million (TPM) to the transcript and gene level using tximport.⁴⁹ 

Quality Control

We first removed samples with a Salmon calculated mapping rate less than 40%. This value was selected as being the far left tail of the distribution of mapping rates across samples (Supplementary Fig. S2). We removed lowly expressed genes by calculating the median expression across all samples for each gene and kept genes that had a median count >200 across all samples. To reduce the noise from experimental variability between each study, we normalized samples by sequence library size using the calcNormFactors function from the edgeR R package, and then quantile smoothed expression data using the R package qsmooth at the tissue level.^50,51 In a change from our previous eyeIntegration work,⁴² we now correct our counts for mapping rate and tissue type with the limma batchEffects function.⁵² The transformed values are used for the box plot and t-SNE visualizations. 

To identify outliers we followed an approach similar to a method reported by Wright et al.⁵³ Briefly, we first selected the 3000 genes with the highest variance across all samples and then for each subtissue type Display Formula\(\def\upalpha{\unicode[Times]{x3B1}}\)\(\def\upbeta{\unicode[Times]{x3B2}}\)\(\def\upgamma{\unicode[Times]{x3B3}}\)\(\def\updelta{\unicode[Times]{x3B4}}\)\(\def\upvarepsilon{\unicode[Times]{x3B5}}\)\(\def\upzeta{\unicode[Times]{x3B6}}\)\(\def\upeta{\unicode[Times]{x3B7}}\)\(\def\uptheta{\unicode[Times]{x3B8}}\)\(\def\upiota{\unicode[Times]{x3B9}}\)\(\def\upkappa{\unicode[Times]{x3BA}}\)\(\def\uplambda{\unicode[Times]{x3BB}}\)\(\def\upmu{\unicode[Times]{x3BC}}\)\(\def\upnu{\unicode[Times]{x3BD}}\)\(\def\upxi{\unicode[Times]{x3BE}}\)\(\def\upomicron{\unicode[Times]{x3BF}}\)\(\def\uppi{\unicode[Times]{x3C0}}\)\(\def\uprho{\unicode[Times]{x3C1}}\)\(\def\upsigma{\unicode[Times]{x3C3}}\)\(\def\uptau{\unicode[Times]{x3C4}}\)\(\def\upupsilon{\unicode[Times]{x3C5}}\)\(\def\upphi{\unicode[Times]{x3C6}}\)\(\def\upchi{\unicode[Times]{x3C7}}\)\(\def\uppsy{\unicode[Times]{x3C8}}\)\(\def\upomega{\unicode[Times]{x3C9}}\)\(\def\bialpha{\boldsymbol{\alpha}}\)\(\def\bibeta{\boldsymbol{\beta}}\)\(\def\bigamma{\boldsymbol{\gamma}}\)\(\def\bidelta{\boldsymbol{\delta}}\)\(\def\bivarepsilon{\boldsymbol{\varepsilon}}\)\(\def\bizeta{\boldsymbol{\zeta}}\)\(\def\bieta{\boldsymbol{\eta}}\)\(\def\bitheta{\boldsymbol{\theta}}\)\(\def\biiota{\boldsymbol{\iota}}\)\(\def\bikappa{\boldsymbol{\kappa}}\)\(\def\bilambda{\boldsymbol{\lambda}}\)\(\def\bimu{\boldsymbol{\mu}}\)\(\def\binu{\boldsymbol{\nu}}\)\(\def\bixi{\boldsymbol{\xi}}\)\(\def\biomicron{\boldsymbol{\micron}}\)\(\def\bipi{\boldsymbol{\pi}}\)\(\def\birho{\boldsymbol{\rho}}\)\(\def\bisigma{\boldsymbol{\sigma}}\)\(\def\bitau{\boldsymbol{\tau}}\)\(\def\biupsilon{\boldsymbol{\upsilon}}\)\(\def\biphi{\boldsymbol{\phi}}\)\(\def\bichi{\boldsymbol{\chi}}\)\(\def\bipsy{\boldsymbol{\psy}}\)\(\def\biomega{\boldsymbol{\omega}}\)\(\def\bupalpha{\unicode[Times]{x1D6C2}}\)\(\def\bupbeta{\unicode[Times]{x1D6C3}}\)\(\def\bupgamma{\unicode[Times]{x1D6C4}}\)\(\def\bupdelta{\unicode[Times]{x1D6C5}}\)\(\def\bupepsilon{\unicode[Times]{x1D6C6}}\)\(\def\bupvarepsilon{\unicode[Times]{x1D6DC}}\)\(\def\bupzeta{\unicode[Times]{x1D6C7}}\)\(\def\bupeta{\unicode[Times]{x1D6C8}}\)\(\def\buptheta{\unicode[Times]{x1D6C9}}\)\(\def\bupiota{\unicode[Times]{x1D6CA}}\)\(\def\bupkappa{\unicode[Times]{x1D6CB}}\)\(\def\buplambda{\unicode[Times]{x1D6CC}}\)\(\def\bupmu{\unicode[Times]{x1D6CD}}\)\(\def\bupnu{\unicode[Times]{x1D6CE}}\)\(\def\bupxi{\unicode[Times]{x1D6CF}}\)\(\def\bupomicron{\unicode[Times]{x1D6D0}}\)\(\def\buppi{\unicode[Times]{x1D6D1}}\)\(\def\buprho{\unicode[Times]{x1D6D2}}\)\(\def\bupsigma{\unicode[Times]{x1D6D4}}\)\(\def\buptau{\unicode[Times]{x1D6D5}}\)\(\def\bupupsilon{\unicode[Times]{x1D6D6}}\)\(\def\bupphi{\unicode[Times]{x1D6D7}}\)\(\def\bupchi{\unicode[Times]{x1D6D8}}\)\(\def\buppsy{\unicode[Times]{x1D6D9}}\)\(\def\bupomega{\unicode[Times]{x1D6DA}}\)\(\def\bupvartheta{\unicode[Times]{x1D6DD}}\)\(\def\bGamma{\bf{\Gamma}}\)\(\def\bDelta{\bf{\Delta}}\)\(\def\bTheta{\bf{\Theta}}\)\(\def\bLambda{\bf{\Lambda}}\)\(\def\bXi{\bf{\Xi}}\)\(\def\bPi{\bf{\Pi}}\)\(\def\bSigma{\bf{\Sigma}}\)\(\def\bUpsilon{\bf{\Upsilon}}\)\(\def\bPhi{\bf{\Phi}}\)\(\def\bPsi{\bf{\Psi}}\)\(\def\bOmega{\bf{\Omega}}\)\(\def\iGamma{\unicode[Times]{x1D6E4}}\)\(\def\iDelta{\unicode[Times]{x1D6E5}}\)\(\def\iTheta{\unicode[Times]{x1D6E9}}\)\(\def\iLambda{\unicode[Times]{x1D6EC}}\)\(\def\iXi{\unicode[Times]{x1D6EF}}\)\(\def\iPi{\unicode[Times]{x1D6F1}}\)\(\def\iSigma{\unicode[Times]{x1D6F4}}\)\(\def\iUpsilon{\unicode[Times]{x1D6F6}}\)\(\def\iPhi{\unicode[Times]{x1D6F7}}\)\(\def\iPsi{\unicode[Times]{x1D6F9}}\)\(\def\iOmega{\unicode[Times]{x1D6FA}}\)\(\def\biGamma{\unicode[Times]{x1D71E}}\)\(\def\biDelta{\unicode[Times]{x1D71F}}\)\(\def\biTheta{\unicode[Times]{x1D723}}\)\(\def\biLambda{\unicode[Times]{x1D726}}\)\(\def\biXi{\unicode[Times]{x1D729}}\)\(\def\biPi{\unicode[Times]{x1D72B}}\)\(\def\biSigma{\unicode[Times]{x1D72E}}\)\(\def\biUpsilon{\unicode[Times]{x1D730}}\)\(\def\biPhi{\unicode[Times]{x1D731}}\)\(\def\biPsi{\unicode[Times]{x1D733}}\)\(\def\biOmega{\unicode[Times]{x1D734}}\)\(T\) and each sample Display Formula\(i\) in Display Formula\(T\), we first calculated Display Formula\({r_i}\), the average correlation between Display Formula\(i\) and all other samples in Display Formula\(T\). Next, we calculated Display Formula\({D_i}\), where Display Formula\({D_i} = {{\left( {{r_i} - r} \right)} \over {median\left( {{r_i} - r } \right)}}\) and Display Formula\(r \) is the grand mean of all Display Formula\({r_i}\) for Display Formula\(i\) in Display Formula\(T\). We removed samples with Display Formula\({D_i} \lt \) −17.5; we determined this threshold by generating a tSNE plot of our samples, and visually identifying outliers in adult retina tissue. The (Display Formula\(max\left( {D_i} \right)\)) among these outliers was −17.58 and from this we chose −17.5 as our outlier threshold. 

To calculate Pearson correlation (Display Formula\({R^2}\)) between GTEx-calculated TPM gene values and our GTEx TPM gene values, we downloaded “GTEx_Analysis_2016-01-15_v7_RNASeQCv1.1.8_gene_tpm.gct.gz” and matched against our GTEx TPM values, running the Pearson correlation with Display Formula\(log2\left( {TPM + 0.01} \right)\) values as per Zhang et al. with the cor function in R.⁵⁴ 

Differential Gene Expression Analysis and GO Term Enrichment

We used the nontransformed length scaled TPM values to determine differential gene expression between different subtissue types. First, we generated a synthetic body set to serve as single representative subtissue type for pan-body gene expression by randomly sampling GTEx tissues. We used the voom function from the limma R package to convert gene expression to precision weights, and then performed pairwise differential expression tests for all combinations of eye subtissues (using mapping rate as a covariate), the synthetic body tissue, and human body tissues using an empirical Bayes test.^52,55 We extracted significant genes (FDR P < 0.01) for all 882 comparisons and used these to calculate GO enrichment. The significant gene list for each eye subtissue was split into upregulated and down regulated sets and each set was tested for enrichment using the enrichGO function from the clusterProfiler R package (q-value < 0.01).⁵⁶ 

eyeIntegration Web App and R Package

The data generated in the above steps is consolidated into a SQLite database, with the original dataset for eyeIntegration and the new 2019 EiaD dataset each getting a separate database file. The code that creates the eyeIntegration web app is written in Shiny and R and has been wrapped into an R package (available in the public domain at https://github.com/davemcg/eyeIntegration_app/) that can be deployed on a local computer or a web server (available in the public domain at https://eyeIntegration.nei.nih.gov). The app can be deployed on a local computer with 50 GB of free disk space by running three commands in R: “devtools::install_github(‘davemcg/eyeIntegration_app')”, “eyeIntegrationApp::get_eyeIntegration_datasets()”, and “eyeIntegrationApp::run_eyeIntegration()”. 

Snakemake Reproducible Pipeline

While the sample search and metadata parsing in a semicurated process, the processing from the raw data to the creation of the SQLite EiaD database underlying eyeIntegration is wrapped in a Snakemake pipeline, which ensures full reproducibility of the results.⁵⁷ We make the code for the pipeline available at https://github.com/davemcg/EiaD_build. 

scRNA-seq Processing

The eyeIntegration site, as of June 2019, hosts two large scRNA-seq datasets from Macosko et al.⁵⁸ and Clark et al.⁵⁹ We use the processed gene count data directly from each group, as well as their cluster assignments, which specify what cell type each individual cell is. The count data are mean averaged to the cell type, age, and gene level for the single cell expression section of eyeIntegration. We also displayed t-SNE and UMAP-based two-dimensional visualizations of the Macosko and Clark data, respectively, in the web app. For detail so the t-SNE processing we did on the Macosko dataset, see the methods of Bryan et al.⁴² 

Power Calculation

We used the ssizeRNA R package to calculate power (p) across samples (n) at an FDR of 0.05.⁶⁰ Important parameters for ssize RNA include the variability (dispersions for the samples and genes), which were calculated directly from our EiaD length scaled TPM values by the edgeR packages estimateCommonDisp and estimateTagwiseDisp. The code to calculate the power is given as ‘power_calc.R.' 

Manuscript as Code and Reproducibility

The figures, tables, and most numbers, are all created and laid out in an R markdown document that interweaves code and text. The knitr and pandoc program is used to lay out the figures and tables and output a docx file. The code that generates this study can be found at https://github.com/davemcg/eyeIntegration_v1_app_manuscript. 

The relevant code-bases (https://github.com/davemcg/eyeIntegration_v1_app_manuscript, https://github.com/davemcg/EiaD_build) and the EiaD dataset itself have been deposited into Zenodo with accession 10.5281/zenodo.3238677 to ensure the data can be accessed in the future, even should eyeIntegration and GitHub become inaccessible in the future. 

Results

EiaD 2019 Contains 24 New Human Eye RNA-seq Studies, 448 New Retina AMD Samples, 207 New Eye Samples, and 16 Total Eye Subtissue Types

Our query on May 8, 2019 to the SRA found 107 potentially relevant studies. We removed nonpertinent studies and selected healthy or unmodified tissue from each relevant study for a total, including 46 studies, 30 of which are new to the 2019 EiaD dataset. The 2019 EiaD dataset contains 835 human eye tissue samples and also includes 1314 GTEx samples across 54 tissues for easy comparison (Table; Supplementary Table S1). The 2019 EiaD contains six undifferentiated iPSC, 56 cornea, four lens, 648 retina, and 121 RPE (choroid) samples; in total we have added 655 new samples to the 2019 EiaD (Fig. 1). We refer to native-tissue extracted RPE as RPE (choroid) because it is not possible to remove the choroid from the RPE without culturing. 

Table

 

EiaD Contains a Large Set of Diverse Eye Tissues, Including Embryonic Stem Cells (ESC)

 

Figure 1

 

Substantial increase in eye tissue count and type from 2017 (180, lighter color) to 2019 (835, darker color) EiaD. We also improved the metadata labeling, the cornea samples (green) now delineates endothelial and epithelial tissues and the retina samples (orange) distinguish retina organoid and RGC from stem cells. Counts for each bar plot given in the boxes. The y-axis is a log2 transformed count of samples passing our QC filters.

Stem cell-derived cornea, stem cell-derived lens, and fetal retina are three new types of subtissues that are now available in EiaD. We also have substantially improved the granularity of the cornea tissue metadata, now delineating whether the tissue is from the endothelium or epithelium (Fig. 1); previously these had been grouped together as adult tissue. Another substantial addition to the 2019 EiaD are nonprotein coding genes; while protein-coding is the most common gene and transcript typse, there are dozens of different noncoding classes. The 2017 version of eyeIntegration only quantified protein coding genes and transcripts. We now quantify expression across 41 gene and 45 transcript types, including protein coding, retained intron, lincRNA, antisense, and pseuodogenes (Supplementary Table S2). 

We also have added the large retina AMD postmortem Ratnapriya et al.³⁷ cohort to EiaD 2019. This cohort contains hundreds of samples ranging from non-AMD (Minnesota Grading System [MGS] 1) to severe AMD (MGS 4). While eyeIntegration is intended to be a source for normal tissues, we have made an exception for this study, as this is a large cohort and AMD is a common disease. We found our corrections methods did not group the non-AMD Ratnapriya et al.³⁷ samples with our other collected retina samples (see Retina MGS in Supplementary Fig. S3). This may be related to the lower mapping rate of the Ratnapriya et al.³⁷ data (see Retina MGS in Supplementary Fig. S2). 

467 More GTEx Samples and Nine New GTEx to 2019 EiaD

Our previous dataset for eyeIntegration version 0.6 held approximately 20 samples per GTEx tissue type. We ran power calculations to assess our ability to detected ≥1 log2(Fold Change) in gene expression between two conditions to determine whether this is a sufficient number of samples (Supplementary Fig. S4). Our calculations suggested, for example, that we had 83% power to detect a 1 log2(Fold Change) difference in gene expression with two groups of 20 samples. To increase our power to make significant eye-to-body comparisons, we added approximately 10 more samples per GTEx tissue types (which at 30 samples, would give approximately 90% power). We also took this opportunity to add bladder, bone marrow, cervix uteri, fallopian tube, ovary, prostate, testis, uterus, and vagina GTEx tissue samples (Supplementary Table S1). 

Rigorous Quality Control and Reproducible Workflow System Ensures High Quality Transcriptomes That Consistently Cluster Together by Tissue Type

We built an automated pipeline for processing and analyzing all data for the web app using the program Snakemake, a python-based workflow management system that allows for efficient parallel execution of the analysis, facilitates reproduction by others, and simplifies long-term maintenance of the EiaD data (Fig. 2; Supplementary Fig. S5).⁵⁷ To create a high quality final dataset across the 2291 initial samples (Supplementary Table S3) and 67,315,523,736 reads, we developed a rigorous quality control procedure as part of our analysis, considering a sample's read mapping rate and median count level as well as behavior relative to samples of the same subtissue type (see Methods). To identify tube mislabeling or sample extraction issues, we used sample-level gene correlation metrics (see Methods) to identify variability within samples of the same subtissue and ensure overall consistency in data processing (Fig. 3). After these steps 81 eye samples and 61 GTEx samples were removed. 

Figure 2

 

Raw RNA-seq data from the SRA is run through our pipeline to create the EiaD, which is used by eyeIntegration app to serve interactive gene expression visualizations across 73 tissues.

Figure 3

 

t-SNE two-dimensional transcriptome profiles by sample demonstrate effective quality control and transcriptome processing. Colors match different tissue types and shapes of points define the origin of the tissues.

To ensure there are no substantial differences in quantification of gene TPM values, we calculated the Display Formula\({R^2}\) between GTEX and EiaD generated TPM values for our shared GTEx samples (see Methods); we computed an Display Formula\({R^2}\) of 0.89. Zhang et al.⁵⁴ reported that RNA-seq quantifications done between alignment-free methods (used in EiaD) and alignment-based methods (used by GTEx) get a Display Formula\({R^2}\) ranging from 0.89 to 0.93. As Zhang et al.⁵⁴ compared quantification methods with identical gene references (we use Gencode GRCh38 gene models and GTEx uses hg19) and did not scale TPM score differently, our result falls in line with expectations. 

After our quality control and processing workflow, we found that samples of the same tissue type and origin cluster well together (Fig. 3; Supplementary Table S4). For example, in the retina group, primary adult tissue clusters tightly and distinctly from other cell types, and retinal organoids and fetal retina samples cluster together. Our ability to uniformly cluster data by known biological source independent of study origin demonstrates that our workflow can effectively account for technical variation between studies. 

While t-SNE is a powerful algorithm for grouping samples, it is not consistent for determining the relationships between clusters.⁶¹ PCA is more useful in this regard. We ran a PCA dimensionality reduction (Supplementary Fig. S6) on all samples, finding that the eye tissues still generally group together and apart from all other human body tissues. Adult retina is most similar to the brain tissue. RPE and cornea are most similar to blood, bone marrow, and skin. 

The eyeIntegration Web App Provides Interactive Visual Portal to All Data

The EiaD 2019 dataset is used directly by the eyeIntegration web app (available in the public domain at https://eyeIntegration.nei.nih.gov). The web app was designed to provide a simple interface that has the same general concept – select specific genes and tissue and view relevant information. The web-app is divided into four general categories: expression, two-dimension sample relationships, gene networks, and data tables. 

Custom Gene and Tissue Expression Boxplots

The ‘Expression' tab of the webpage provides a wealth of information about gene- and transcript-level expression for eye and body tissues, giving the user the ability to compare the expression of different genes within a single tissue, as well as the expression of genes across multiple tissues (Fig. 4A). The user first selects either the 2017 or 2019 gene or transcript EiaD dataset, then Hugo Gene Nomeclature Committee (HGNC) genes names (or ENSEMBL transcripts), then tissues. A boxplot then is generated after hitting the “Re(Draw) Plot” button with overlaid individual data points. On mouse-over, the metadata for the individual sample is displayed. A tabular report is generated based on selected genes and tissues: a table with links to Ensembl, GeneCard, and OMIM for each gene for quick referencing, and a table containing expression levels for each selected gene in each selected tissue. The tables can be arranged or sorted to the user's preference and can be easily downloaded for local use. 

Figure 4

 

Screenshots from eyeIntegration web app. (A) Pan-tissue gene expression box plots with accompanying data tables. The data tables display the rank (lower is more highly expressed) of each gene in each sub tissue, decile of the rank (10 is the highest decile of expression), and gene's mean log2(TPM + 1) score for each sub tissue. (B) Heatmap visualization.

Heatmap built by the R package ComplexHeatmaps based on expression can be drawn for selected genes and tissues and gene expression can be compared across many genes and tissues (Fig. 4B).⁶² Finally, a session can be saved or shared by building a custom link for the session with the “Build URL Shortcut” button. 

Differential Expression and Gene Ontology Enrichment Tests Allow Quick Comparison of Gene Differences Between Groups

We performed multiple differential comparisons at the subtissue level within all eye tissues and against a pan-body synthetic set comprised of a stratified sample of all tissues present in our subset of the GTEx dataset, allowing quick identification of eye-specific genes across 882 different comparisons. We expanded the differential tests in the 2019 EiaD by adding the GTEx tissues as direct comparisons to our eye subtissues. The user can view the results selecting ‘Differential' under the ‘Expression' tab (Supplementary Fig. S7D). As with ‘Expression,' the user can select which version of the web app to draw data from as well as select for gene- or transcript-level comparisons. The user additionally has the option to select different gene classes to examine, for example, protein coding, lincRNA. 

The results of differential expression are presented in a tabular format showing log₂ fold change, average expression, and P values. Depending on the comparison, there are 1 to 33,380 differentially expressed genes (SupplementaryFiles). The table can be easily searched for any given gene, viewed and ordered to the user's preference, and downloaded in CSV format. Differential expression can be visualized through fold change bar graphs with the ‘Pan-tissue plots' selection under ‘Expression.' Additionally, we performed GO enrichment for all differential comparisons. Enriched GO terms are presented first as a word cloud, for quick comparison of GO enrichment. We provide tables, with similar viewing options as the differential expression table, for enriched GO terms in each class of a given differential comparison. 

Murine scRNA-seq Enables Testing of Retina Cell Type Specific Expression

We incorporated scRNA-seq data from murine retina across two studies.^58,59 This allows researchers to quickly examine gene expression across individual cell types in the retina. Single cell gene expression data are visualized through a heatmap showing the expression of a gene across multiple retinal cell types and different developmental time points, from embryonic day (E)11 to postnatal day (P)14 (when available), and a table of expression values is generated containing the expression data used to draw the heatmap (Fig. 4C). We also provide t-SNE/UMAP–based clustering using cell type–specific labeling created by the publishing authors (Fig. 4D, see Methods). The plots show all cell types present at a given developmental stage, and highlights cells expressing a gene above a user-selected given level. 

EiaD 2019 Suggests That iPSC-Derived Organoids and Fetal Retina Have Closely-Related Transcriptomes

There are currently two major approaches to studying developing human retina: postmortem fetal tissue and stem-cell derived organoids. We looked at how well these approaches to studying developing retina compare at a transcriptomic level, for tissue–organoid relationships and how well they correlate across early development. 

To evaluate how the tissues and organoids compare at a transcriptome level, we looked at the same t-SNE plot from Figure 3 and focused in on the three types of retina tissue (adult, fetal, and organoid; Fig. 5A). Here, we saw three distinct groupings: adult retina (1), developing fetal retina and stem cell-derived organoid (2), and undifferentiated and early differentiating stem cells (3). We identified several organoid samples in cluster 3, but these share one important difference from the rest of the organoid samples in cluster 3: they have been differentiating for less than 30 days (shape ‘X'). All of the organoid retina samples in cluster 2 were older than 50 days. 

Figure 5

 

Organoid retina, stem cell retina, and fetal retina tissue have highly similar transcriptomes. The zoom inset (A) shows the retina samples. The “Subtissue Cluster” shading shows the cluster membership of the three major groups. The shapes of the points show the different origin types – notable types include the square for adult, the ‘X' for organoid under 30 days of differentiation, and the diamond for organoid over 30 days of differentiation. (B) Major markers of retina progenitor, photoreceptors (cone and rod), and RGCs have similar gene expression patterns across development in retina fetal tissue and organoids.

To assess how similarly the fetal and organoid retina develop through time, we plotted expression of retinal progenitors, photoreceptors, and retinal ganglion markers by time in days (Fig. 5B). Each row is a gene marker of either retinal progenitor, photoreceptor, or RGC. The rows are hierarchically clustered to put more similar expression patterns closer together, as denoted by the height of the dendrogram. We split the organoid tissues into three groups: Kaewkhaw et al.¹² GFP+ and GFP- samples, and Eldred et al.⁶³ samples. The Kaewkhaw samples are flow sorted for a GFP marker (GFP+) under the control of the CRX promoter, an important regulator of photoreceptor development. GFP+ cells would be enriched in photoreceptor populations. We saw that the retinal progenitor, photoreceptor, and RGC groups are largely clustered together, with patterns of expression consistent across the fetal retinal and organoid groups. 

Differential Gene Expression of Organoid Retina Versus Fetal Tissue Identifies Sets of Genes Relating to Patterning (HOXB Family), Cell Adhesion (Protocadherin Family), and RGC Identity (BRN3/POU4F, NEFL, GAP43, SNCG)

To identify specific changes between retinal organoid and fetal retina tissue, we performed differential gene expression and GO term enrichment analyses. The GO term enrichment identified cell adhesion (protocadherins) and patterning (HOXB family) as enriched gene sets in retinal organoids. As there is some evidence suggesting that protocadherins influence RGC viability and we noticed that several RGC markers appeared to have lower expression in the organoids compared to the fetal tissue Figure 5B, we looked more closely into RGC marker expression.⁶⁴ 

We plotted HOXB family, protocadherin family, and RGC genes in a heatmap visualization, with columns as age in days of fetal or organoid retina. Rows are genes, split by the three different groups of genes and are internally clustered by how similar the expression patterns are. We observed that there are strong, consistent gene expression differences in these three groups of genes between fetal retina and the organoid samples (Supplementary Fig. S8). We also plotted the differential expression values between all organoids and all fetal retina samples; all genes across all three groups are significantly differentially expressed with an FDR corrected P < 0.01. 

Limitations of the RNA-seq Quantification in EyeIntegration

Salmon quantification, while highly performant and accurate, has a higher variance for lower read depth samples and shorter transcripts.⁵⁴ Extra care should be taken with comparisons with lower counts of samples (cornea, RGC) as smaller sample numbers decrease the confidence in differential expression. We do not recommend you directly compare our TPM values with your counts data as there are many important variables that will differ. Instead run our Snakemake pipeline (available in the public domain at https://www.github.com/davemcg/EiaD_build), adding your samples. Finally, we would like to remind any users that RNA-seq methods measure mRNA levels, but the functional unit is the protein; Westerns are still the gold standard with which to evaluate expression and localization. 

Data Accessibility

Individual data files for gene expression and sample metadata can be downloaded from the ‘Data' tab on the web app. All data and code used to generate the web app can be installed from the R command line by running devtools::install_github(‘davidmcg/eyeIntegration_app'). The code for the EiaD data processing pipeline can be found at https://github.com/davemcg/EiaD_build. 

Discussion

EiaD 2019 contains a large set of carefully curated, reproducibly processed human eye RNA-seq datasets alongside a human body tissue comparison set from the GTEx project. It is available for local install as an R package at https://www.github.com/davemcg/eyeIntegration_app and it is served via a web app, eyeIntegration at https://eyeIntegration.nei.nih.gov. The web app offers a wide range of user-driven visualizations to compare expression of genes across dozens of human body and eye tissues. Furthermore, murine scRNA-seq datasets have been incorporated, allowing for examination of retina cell type–specific gene expression. Several human and nonhuman primate studies have been posted in the past year on the preprint server bioRxiv and as the raw data becomes publicly available, we will be updating this section of eyeIntegration.^65–67 

If you wish to have your data added to EiaD in the future, we suggest you: (1) deposit data into GEO/SRA; (2) use clear, descriptive, consistent, and detailed metadata for each sample; and (3) (optional) contact the corresponding author. Contacting the corresponding author is only necessary if you feel your data should be included in EiaD and were deposited into the SRA before May 8, 2019. 

As human fetal tissue is difficult to obtain and, thus, not very amenable for chemical or genetic modification, it is crucial for organoid-based models to be developed. Our merging of these datasets and analysis at the transcriptome level (compared to cross-analyzing using a limited number of known marker genes) indicated that these two approaches successfully recapitulate fetal retina tissue, to a first approximation, at the whole transcriptome level. However as organoids do not develop to full function, it is important to look at how gene expression differs between retinal organoid and fetal tissue so as to suggest areas for improvement. 

We used our large dataset to narrow in on three core processes that differ significantly and substantially between retinal organoids and fetal retina. First we showed that the HOXB family is overexpressed in the organoids. The homeobox family is well known to initiate polarity of the embryo during early development.⁶⁸ Retinoic acid is applied at approximately day 20 in culture to help differentiate stem cell to organoids and also is known to activate gene members of the HOXB family. The lack of HOXB expression at any age in fetal retina and the broad chromatin and gene expression changes HOXB family members can mediate suggests that HOXB activity may be unwanted for organoid maturation. 

Next, we detected several protocadherins more highly expressed in the fetal tissue, relative to the organoids. Protocadherins mediate cell-to-cell connections and, in the developing mice, are shown to be important for spinal internneurons and RGC survival.^64,69 We would predict that decreased protocadherin expression reduces the number and maturation of RGC. Indeed, we observed that many canonical RGC markers, while present in detectable levels in the organoids, are significantly underexpressed relative to fetal tissue. This result suggested that modifying culture conditions to promote protocadherin expression may result in higher RGC yield and survival. 

We built the EiaD dataset and the accompanying web app, eyeIntegration in the hopes that easily accessible gene expression across tissue space and time will be a useful tool for hypothesis generation and refinement in eye research. Wrapping all of the data processing steps in a Snakemake pipeline has several important advantages for the community: our code is publicly available for review, our analyses are reproducible, future sample updates can be streamlined in with less effort, and because all the processing is in modular pieces it is easier to add new analysis steps. In the future, we plan on regularly adding new samples to EiaD, offering de novo eye tissue transcriptomes, expanding the single cell RNA-seq expression tooling, adding nonhuman eye samples, and epigenetic datasets. 

Acknowledgements

The authors thank the dozens of groups who provided the raw data required to create this project. We keep a running list of the projects and associated citations at https://github.com/davemcg/eyeIntegration_app/blob/master/inst/citations.md and strongly encourage anyone who uses EiaD and eyeIntegration to cite relevant projects. We also thank Kapil Bharti, Robert Hufnagel, and Brian Brooks for their continuous set of critiques and suggestions in the development of eyeIntegration app over the past 2 years, and Tiziana Cogliati for editing this manuscript. We also thank the two anonymous reviewers for their careful reading and constructive criticisms. Finally, this work used the computational resources of the National Institutes of Health (NIH; Bethesda, MD) HPC Biowulf cluster (available in the public domain at http://hpc.nih.gov). 

Supported by the Intramural Research Program of the National Eye Institute. 

Disclosure: V. Swamy, None; D. McGaughey, None 

Amended August 6, 2019: In the third paragraph of the article, the URL https://eyeIntegration.nei.nih was corrected to be https://eyeIntegration.nei.nih.gov.