{'fraction_non_alphanumeric': 0.03848066023665109, 'fraction_numerical': 0.017201948891319478, 'mean_word_length': 5.092065414900061, 'pattern_counts': {'":': 0, '<': 1, '<?xml version=': 0, '>': 2, 'https://': 0, 'lorem ipsum': 0, 'www.': 3, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 1, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Transcriptome Dynamics of Developing Photoreceptors in Three‐Dimensional Retina Cultures Recapitulates Temporal Sequence of Human Cone and Rod Differentiation Revealing Cell Surface Markers and Gene Networks

Abstract The derivation of three‐dimensional (3D) stratified neural retina from pluripotent stem cells has permitted investigations of human photoreceptors. We have generated a H9 human embryonic stem cell subclone that carries a green fluorescent protein (GFP) reporter under the control of the promoter of cone‐rod homeobox (CRX), an established marker of postmitotic photoreceptor precursors. The CRXp‐GFP reporter replicates endogenous CRX expression in vitro when the H9 subclone is induced to form self‐organizing 3D retina‐like tissue. At day 37, CRX+ photoreceptors appear in the basal or middle part of neural retina and migrate to apical side by day 67. Temporal and spatial patterns of retinal cell type markers recapitulate the predicted sequence of development. Cone gene expression is concomitant with CRX, whereas rod differentiation factor neural retina leucine zipper protein (NRL) is first observed at day 67. At day 90, robust expression of NRL and its target nuclear receptor NR2E3 is evident in many CRX+ cells, while minimal S‐opsin and no rhodopsin or L/M‐opsin is present. The transcriptome profile, by RNA‐seq, of developing human photoreceptors is remarkably concordant with mRNA and immunohistochemistry data available for human fetal retina although many targets of CRX, including phototransduction genes, exhibit a significant delay in expression. We report on temporal changes in gene signatures, including expression of cell surface markers and transcription factors; these expression changes should assist in isolation of photoreceptors at distinct stages of differentiation and in delineating coexpression networks. Our studies establish the first global expression database of developing human photoreceptors, providing a reference map for functional studies in retinal cultures. Stem Cells 2015;33:3504–3518

(D) GFP expression in hESC derived 3-D retina. The green signal at the basal side is autofluorescence, which is normally observed when the photograph is captured from floating aggregates. However, the real GFP+ cells can be easily distinguished on the apical side by cell shape (see arrowhead).
(E) GFP expression in NR section from 3-D d30 retina.
(F) Histogram showing a temporal increase in GFP expression. The bar graph indicates the number of GFP+ cells as differentiation proceeds from d37 to d90. Figure S2. Expression of photoreceptor proteins in 3-D retina cultures, related to figures 3 and 4. GFP expression shows CRX promoter activity, and therefore indicates developing photoreceptors.
(A-B) Immunostaining for OPN1SW (S-opsin) (A) and OPN1SW and ARR3 (B). Arrows indicate co-localization of GFP and ARR3 or OPN1SW. Arrowheads show cells with no co-localization.
(C) Immunostaining for OPN1M/L (M/L opsin) and ARR3. Arrows indicate colocalization of GFP and ARR3 or OPN1M/LW. Arrowheads show cells with no colocalization.
(D) Immunostaining for RHO (rhodopsin) and CRX. Arrows indicate co-localization of GFP and CRX. Arrowheads show cells that express RHO, but not CRX.
(E-G) Immunostaining of OTX2 (E), NR2E3 (F) and NRL (G). (F) includes a higher magnification image showing nuclear staining of NR2E3, which is distinct from the background (appearing as a line) on the apical side. (C) Top ten differentially expressed transcription factor genes, obtained by comparing gene profiles between d37 vs. d47, d37 vs. d67 and d37 vs. d90.   Gene expression pattern of flow-sorted CRXp-GFP cells from human 3-D retina is compared to that of flowed sorted photoreceptors from Nrlp-GFP/wild type (rods) and Nrlp-GFP/Nrl -/-(S-cone-like) mouse retina [1].

Supporting Information Movie
Movie S1. Migration of CRXp-GFP+ cells from basal to apical surface of NR
Human ESCs were passaged by incubating in 1 mg/ml collagenase IV (Invitrogen) at 37°C for 30 min and maintaining the cells in small clumps.

Alignment and Transcript Quantitation
Qualitative assessment of the FASTQ files was evaluated using FastQC (http://www.bioinformatics.babraham.ac.uk/projects/fastqc/). The reference assembly GRCh37.p13 with transcript annotation Gencode 19 was used for alignment and quantitation. Transcript quantitation was accomplished by streaming read alignment from Bowtie2 [2] to eXpress [3][4][5]. Bowtie2 alignment settings were used as previously described [4] and eXpress setting -r-stranded was used for directional transcript analysis. Bias corrected effective counts were used for downstream normalization and differential expression analysis. Read alignments used for viewing in Integrated Genomics Viewer (IGV) [6] were performed using a 2step alignment strategy employing TopHat2 [7] splice junction mapper keeping the alignment with the best alignment score. Reads were first aligned to the annotation GTF prior to aligning to the genome reference assembly using the default settings with the exception of the following options: -g 1 -library-type fr-firststrand -G <gtf file>. This study utilized the high-performance computational capabilities of the Biowulf Linux cluster at the National Institutes of Health, Bethesda, Md. (http://biowulf.nih.gov).

Sample Batch Correction and Data Normalization
Principal component analysis (PCA) was applied to the FPKM values calculated by eXpress [3] to explore if any batch effect exists among GFP (+) samples. The first principle component revealed that our data had a strong batch effect due to replicate production procedure (data not shown). Effective counts from eXpress were used to filter and to normalize the data via library size correction using the TMM method within the R (R Core Team, 2014, http://www.R-project.org/) package edgeR [8][9][10][11]. Prior to normalization transcripts showing at least 1 Count Per Million (CPM) in at least half of the samples within each time point were kept for further analysis. Design matrices were created for defining groups (time course) with and without batch definition, however, only the time course factor was used for TMM normalization. Limma [12] "RemoveBatchEffect" function was executed on log2 transformed CPM data. Batch corrected log2 CPM values were then converted to FPKM values prior to clustering. Affinity Propagation (AP) clustering [13,14] was used to assess the success of the procedure. Row-wise average of undifferentiated sample values was added as the first column to the processed data after differential expression analysis. The same procedure was also applied to samples from undifferentiated and GFP (-) samples.

Differential Expression Analysis
Differential expression was performed using the edgeR likelihood ratio test from on a generalized linear model (glm), with day 37 as the intercept. In the glm, the batch factor was used as a blocking factor. Significance was defined as having a false discovery rate (FDR) less than or equal to 0.05, and 4697 transcripts fit this criterion. Volcano plots of those comparisons were plotted to see the ideal fold change cutoff ( Figure 5E). The colors of the points those represent the mean expression of the gene between the time points compared. The key observation is that the mean expression of a transcript is getting closer to 0 (blue) as the fold change is getting closer to 0 although all the FDRs in the plot equal 0.05. We

Cluster Analysis
Two types of clustering were used throughout the analysis of the data: Affinity Propagation (AP) and hierarchical clustering. Hierarchical clustering was used to reorder the transcripts within a cluster generated by AP. We have used the Pearson's Correlation with simple Euclidean distance as the similarity measure as the input for AP via "CorSimMat()" function. The clustering function was executed for those 3157 significantly differentially expressed and 53 clusters were generated.

GO Analysis
Cluster results were analyzed with clusterProfiler [15]. Level 8 Biological Process (BP) terms, where GO term "visual perception" has been defined as the 8 th child term (by following the shortest path) of BP of GO Directed Acyclic Graph, were identified by getGOLevel function of GO.db R package. The results intersected with the ones produced by clusterProfiler, so that no GO terms left except the ones belonging to 8 th level in GO BP Directed Acyclic Graph. The terms left then was analyzed by their semantic similarity [16,17]. The resulting similarity matrix was analyzed with APcluster again to define the exemplar GO terms (where input parameter q was set as 0.85). The resulting exemplars where searched in the original data frame and the others were eliminated. The resulting data frame was plotted with a customized function written via ggplot2 components.

CRX Network Analysis
Several methods have been developed to infer high confidence networks from multiple high-throughput experimental data sets [18][19][20][21]. We implemented a rank combined approach based framework [18] that allows flexibility in selection of data sets and preferred method/algorithm for individual tasks, such as differential expression analysis for knockout experiments, computing pairwise co-expression scores, etc. This framework was developed in R programming language and implemented a method to compute transcription factor to gene association scores in a ChIP-seq experiment as described in [18]. Limma R package [12] was used to perform differential expression analysis of knockout versus wild-type experiments.
High confidence targets (targets with averaged ranked score ≥ 1.5) of CRX were inferred from RNA-Seq (P21 wild-type and P21 knockout [22]) and ChIP-seq (p56) [23] data sets. Functional and disease analysis of genes that were expressed (>= 1 FPKM in all samples in at least one time point) and non-expressed were performed individually through the use of QIAGEN's Ingenuity Pathway Analysis (IPA ® , QIAGEN Redwood City, www.qiagen.com/ingenuity) tool (data not shown).