Human genomic sequences and RefSeq gene coordinates (version hg19) were obtained from UCSC. All RNA-seq reads were mapped to human genome by Tophat2 (v2.0.12) (Trapnell et al., 2009). Relative expression level (RPKM value) was then calculated by Cufflinks (v1.2.1) using Refseq genes as reference annotation with “--GTF” option (Trapnell et al., 2010). Differentially-expressed genes in organoids to their parental hESCs were defined by 2-fold change (Figure 4B and 4C). GO analysis was then performed to the differentially-expressed genes by GOstats (v2.24.0) in Bioconductor. False discovery rate (FDR) was calculated to each GO term by Benjamini-Hochberg method with p.adjust function in R. GO terms with FDR<0.05 were used as statistical significance.