High-resolution transcriptome analysis with long-read RNA sequencing

Hyunghoon Cho1,‡, Joe Davis2, Xin Li3, Kevin S. Smith3, Alexis Battle1,2,∗,
Stephen B. Montgomery1,2,3,∗

1 Department of Computer Science, Stanford University, Stanford, CA, USA
2 Department of Genetics, Stanford University, Stanford, CA, USA
3 Department of Pathology, Stanford University, Stanford, CA, USA
$\ast$ E-mail<EMAIL_ADDRESS>and<EMAIL_ADDRESS>
$\ddagger$ Present address: Computer Science and Artificial Intelligence
Laboratory, Massachusetts Institute of Technology, Cambridge, MA, USA

## Abstract

RNA sequencing (RNA-seq) enables characterization and quantification of
individual transcriptomes as well as detection of patterns of allelic
expression and alternative splicing. Current RNA-seq protocols depend on high-
throughput short-read sequencing of cDNA. However, as ongoing advances are
rapidly yielding increasing read lengths, a technical hurdle remains in
identifying the degree to which differences in read length influence various
transcriptome analyses. In this study, we generated two paired-end RNA-seq
datasets of differing read lengths (2x75 bp and 2x262 bp) for lymphoblastoid
cell line GM12878 and compared the effect of read length on transcriptome
analyses, including read-mapping performance, gene and transcript
quantification, and detection of allele-specific expression (ASE) and allele-
specific alternative splicing (ASAS) patterns. Our results indicate that,
while the current long-read protocol is considerably more expensive than
short-read sequencing, there are important benefits that can only be achieved
with longer read length, including lower mapping bias and reduced ambiguity in
assigning reads to genomic elements, such as mRNA transcript. We show that
these benefits ultimately lead to improved detection of cis-acting regulatory
and splicing variation effects within individuals.

## Author Summary

Ongoing improvements to next generation sequencing technologies are leading to
longer sequencing read lengths, but a thorough understanding of the impact of
longer reads on RNA sequencing analyses is lacking. To provide a better
understanding of this issue, we generated and compared two RNA sequencing
datasets of differing read lengths (2x75 bp and 2x262 bp) and investigated the
impact of read length on various aspects of analysis, including the
performance of currently available read-mapping tools, gene and transcript
quantification, and detection of allele-specific expression patterns. Our
results indicate that, while the scalability of read-mapping tools and the
cost-effectiveness of long read protocol is an issue that requires further
attention, longer reads enable more accurate quantification of diverse aspects
of gene expression, including individual-specific patterns of allele-specific
expression and alternative splicing.

## Introduction

The application of next-generation sequencing (NGS) to RNA has provided more
complete means to annotate and quantify transcriptomes. Specifically, it has
improved the characterization of many aspects of RNA biology including the
detection of transcription start sites [1, 2, 3], allele-specific expression
[4], alternative splicing events [5], fusion transcripts [6], RNA-editing [7],
and antisense transcription [8]. However, the analytical properties of RNA-seq
data with respect to rapidly changing technological parameters, such as
continually increasing read lengths, remain challenging to ascertain. While it
would be ideal to sequence the full length of each mRNA molecule, current NGS
technologies are limited to analyzing short fragments of cDNA, and only a
limited number of bases can be read from each fragment with reasonable
accuracy. Paired-end sequencing methods were developed to partially address
this limitation, and, notably, it has been shown that the paired-end reads
enable more accurate estimations of transcript abundances [9]. Moreover, it is
established that longer reads will improve performance in several of the noted
applications of RNA-seq, such as transcript assembly, transcript
quantification, and gene fusion detection [10, 11]. However, as public RNA-seq
data of variable read lengths are increasingly available to the research
community, ongoing assessment of the effect of read length on various aspects
of transcriptome analysis will be necessary to understand the practical
implications of integrating and comparing such data. In this study, we
generated paired-end RNA-sequencing datasets for the lymphoblastoid cell line
GM12878 with 262 bp reads (L262) and 75 bp reads (L75) and carried out a
comparative analysis.

One aspect of RNA-seq analysis we focused on was its ability to quantify gene
expression. Although RNA-seq is considered to be a more robust quantification
method than hybridization-based microarray assays due to its wide dynamic
range and ability to detect novel sequences, various sources of bias inherent
to RNA-seq have also been identified, including GC-content [12], transcript
length [13], overdispersion [14] and mappability [15]. While biases arising
from GC-content and transcript length can be countered with appropriate
normalization methods [12, 13], the inability of short reads to effectively
disambiguate fragments from non-unique regions of the genome can only be
mitigated by using longer reads. To evaluate the benefit of longer reads, we
compared the number of reads that uniquely mapped to each gene in both L262
and L75 datasets and investigated the differences that are distinguishable
from random biological variation. We found that the proportion of reads in
L262 that unambiguously aligned to a pseudogene is roughly three times higher
than that of L75, with a similar observation in coding genes with low
mappability.

Quantification of transcript-level expression from RNA-seq data is a
considerably harder problem than that of gene-level expression, since
alternative splicing patterns introduce ambiguity in mapping each read to a
particular transcript. To estimate transcript abundances from short-read data,
several methods based on statistical inference have been developed such as
Cufflinks [16], TIGAR [17], RSEM [18] and FluxCapacitor
(http://ux.sammeth.net/capacitor.html). While these methods aim to make an
informed guess using available information, one could improve the outcome or
even make these tools obsolete by using longer reads to reduce or remove the
overall ambiguity in mapping reads to transcripts. In this study, we looked at
the number of mRNA isoforms each read maps to and found that long reads lead
to better specificity. In particular, 32.98% of exonic reads in L262 were
unambiguously assigned to a single isoform, compared to 25.04% for L75. To
further evaluate this class of methods, we evaluated the degree to which two
different methods, Cufflinks and FluxCapacitor, predicted transcript
abundances from the L75 data that were directly measured in the L262 data.

RNA-sequencing also supports the quantification of allele-specific expression
(ASE) and allele-specific alternative splicing (ASAS), providing direct
evidence of cis-regulatory effects within an individual. Analysis of allele-
specific patterns has the potential to reveal additional regulatory variants
and mechanisms, but analysis in RNA-seq is limited to evaluation of reads
containing at least one heterozygous site, which is needed to differentiate
between the two copies of a gene. This significantly reduces the effective
library size and consequentially many tests are statistically underpowered.
Therefore, long reads are useful due to a higher chance of containing a
heterozygous site. Here, we report the degree of improved detection of ASE and
ASAS from longer read data. We further identified compelling cases of ASAS,
which give us a glimpse of the complex regulatory mechanisms within
individuals.

One important practical consideration regarding read length is cost-
efficiency. Specifically, for any given analytical goal, it is not a priori
clear whether it is more cost-effective to use a long read protocol, or a
short read protocol with higher depth. We provide some insight into this issue
by carrying out a cost comparison between L262 and L75, generated on a MiSeq
and HiSeq, respectively. Based on current sequencing prices, it is more cost-
effective to select a short-read protocol for many analyses. Specifically,
obtaining a target number of reads with isoform-specificity can be achieved at
lower total cost with short reads, even though a larger number of reads must
be sequenced; a similar result is seen for obtaining haplotype-specific reads
for ASE analysis. On the other hand, the reduction in mapping bias, which
ultimately leads to more accurate quantification, is only accessible via
longer read lengths and cannot be achieved by simply increasing the number of
short reads. More broadly, some elements (genes, transcripts, or loci) in low
mappability regions would remain under-represented in short read data,
regardless of depth or cost. As the price gap between long-read and short-read
RNA sequencing platforms becomes smaller, we expect the benefits of longer
read lengths to become more apparent.

## Materials and Methods

### Sample preparation and sequencing

GM12878 is a CEPH/UTAH family EBV transformed peripheral blood B lymphocyte
cell line (Cat# XC01463) purchased from The Coriell Institute (Camden, NJ) and
grown in RPMI 1640 supplemented with 10% FCS and penicillin/streptomycin in
humidified 5% CO2 at a concentration of approximately $1x10^{6}$ cells/ml.
Total RNA was isolated from $1x10^{7}$ GM12878 cells using Trizol. RNA quality
was assessed with the Agilent Bioanalyzer 2100 and RIN scores above 9 were
used for cDNA production. One microgram total RNA was used to isolate poly A
purified mRNA which was used for library construction using the Illumina
TruSeq RNA Preparation Kit. RNA fragmentation time was 8 minutes for samples
to be used for 75 bp paired-end reads and 1 minute for 262 bp paired-end
reads. Strand specificity was performed using dUTP during second strand
synthesis. GM12878 cDNA was indexed with Illumina adapters for use in
multiplex HiSeq runs. GM12878 cDNA sequenced by an Illumina HiSeq yielded
25,933,924 75 bp paired-end reads (L75). For MiSeq sequencing, 6 pM cDNA was
denatured and loaded into a 500 cycle MiSeq sequencing cartridge that yielded
9,524,186 262 bp paired-end reads (L262). All data is freely available from
GEO under the accession identifier GSE57862.

### Read alignment

Based on the observation that many read pairs in L262 come from cDNA fragments
shorter than the read length (262 bp) (Supplementary Figure 1), we identified
and merged such pairs by searching for a 13 bp adapter sequence, which
immediately follows the fragment sequence when the fragment is shorter than
the read length. During merge, if the two reads disagreed on a base, we chose
the base with a higher quality score. This resulted in a bipartition of L262
into a single-end library of 4,470,824 reads (L262S) and a paired-end library
of 5,051,896 read pairs (L262P). 1,466 read pairs were discarded for being
shorter than 20 bp after the merge or having more than 10 mismatches. L262P
and L262S were aligned to the reference separately and then aggregated for
comparison with L75. Most of the read pairs in L75 came from fragments longer
than the read length (Supplementary Figure 1), and thus we did not merge reads
in L75. We aligned each library using the following tools with default
settings: TopHat2 (version 2.0.5) [19], GSNAP (version 2012-07-20) [20],
Bowtie2 (version 2.0.0-beta7) [21], and STAR (version 2.3.0) [22]. For tools
that required a list of known splice sites, we used the reference gene
annotations from GENCODE 15 [23].

### Concordance analysis

The VCF file containing positions genotyped for GM12878 was obtained from the
1000 Genomes Project (http://www.1000genomes.org). After discarding
insertions, deletions (indels), and sites with incomplete genotypes, roughly
5.4M sites remained for analysis. For each dataset of aligned reads, genotyped
sites covered by at least 10 reads (20-40K sites) were analyzed for
concordance between the DNA-seq based genotype from the VCF file and the
observed genotype from the RNA-seq reads. For each genotyped site $i$, let
$n_{i}$ be the total number of reads covering the site, and $x_{i}$ be the
number of mismatched reads, which is defined as reads that contain an allele
not supported by the known genotype. Under the null hypothesis that there is
no mapping error, the number of mismatched reads follows
$\text{Binomial}(n_{i},\epsilon)$, where $\epsilon$ is the per-base sequencing
error rate. The discordance $p$-value is defined as $P(X\geq x_{i})$, where
$X\sim\text{Binomial}(n_{i},\epsilon)$. Intuitively, the more significant this
$p$-value is, the more confident we are that mapping error affected RNA-seq
reads at the corresponding site. We estimated $\epsilon$ separately for each
dataset by dividing the total number of mismatches across all sites by the
total number of aligned bases (Supplementary Figure 3). Then, for each mapper,
we compared the distribution of discordance $p$-values between L262 and L75
using a one-sided Wilcoxon rank-sum test. We controlled for differences in
statistical power and technical variance by stratifying the sites into four
groups based on read depth and performing a separate comparison for each bin:
1) 10-99, 2) 100-499, 3) 500-999, and 4) $\geq$1000 reads. The read depth
distribution was highly similar across read lengths for each bin
(Supplementary Figure 4).

### Mapping reads to genomic features

To assign each read to a genomic feature — gene or mRNA isoform — for the
purposes of quantifying gene abundance or analyzing isoform disambiguation, we
constructed a feature set based on the annotations from GENCODE 15 [23] for
each locus in the reference genome representing the features that contain the
corresponding nucleotide. Then, for each aligned read, we took the
intersection of all nonempty feature sets associated with the bases covered by
the read. Only considering the nonempty sets allows the reads to have bases
that are not explained by the reference annotations, as long as the reads
contain other informative bases. If the resulting set of features (which we
refer to as _consistent_ features) has a single element, it implies that the
read pair is unambiguously mapped to a feature. For mRNA isoforms, we also
incorporated information about splice junctions to obtain a more precise
feature mapping; when a read alignment included one or more splicing events,
we only considered the set of isoforms whose splicing patterns agree with the
read.

### Detection of allele-specific expression patterns

Phased genotypes for NA12878 were inferred by family relationships using Ped-
IBD [24], while discarding sites with Mendelian inconsistencies. For each
aligned read, we examined all heterozygous sites that appear in the read and
assigned the read to either the paternal or maternal haplotype only if all
heterozygous sites have the corresponding alleles. Gene-level ASE was assessed
with a binomial test based on the number of maternal and paternal reads
assigned to a given gene, assuming a balanced distribution as the null. We
detected ASAS by testing for differential usage of individual exon blocks
using a procedure developed in [9]. First, the exonic regions of each gene
were divided into “exon blocks” such that each block does not contain any of
the exon boundaries of known mRNA isoforms. Then we tested whether each exon
block exhibits a differential usage pattern using a chi-square test on a
2-by-2 contingency table where the rows indicate whether a read is maternal or
paternal and the columns indicate whether a read overlaps with the exon block
of interest (as illustrated in Supplementary Figure 5). We excluded human
leukocyte antigen (HLA) complex genes from the analysis to avoid cases with
possibly incorrect genotypes, and also disregarded genes with low total read
counts (less than 20 for ASE and 40 for ASAS). Significant cases were called
at 10% false discovery rate using the Benjamini and Hochberg procedure [25].

### Calculation of sequencing costs

We estimated the sequencing costs for L262 and L75 as follows. For L262, a 500
cycle MiSeq sequencing cartridge cost $880, which yielded a total of 9,524,186
read pairs ($\sim$4.99 billion bases). For L75, Centrillion Biosciences
charged us $2400 for a single 2x75 lane, which we used to sequence six
samples. Thus, a single L75 sample alone costs roughly $400, which yielded
25,933,924 read pairs ($\sim$3.89 billion bases). The resulting cost estimate
per million bases is $$0.176$ for L262 and $$0.103$ for L75. The cost of
library preparation (approximately $60 for each library) is not included in
the estimate.

## Results

### Performance of read-mapping tools

To analyze the impact of read length on the performance of currently available
read-mapping tools, both L262 and L75 were aligned to the reference genome
(GRCh37/hg19) using Bowtie2, TopHat2, STAR, and GSNAP (Methods). We emphasize
that our goal is not to compare the performance of different methods on a
particular dataset; rather, we aim to analyze how different methods scale to
longer reads. Figure 1a shows the percentage of reads uniquely mapped to the
reference genome. For tools that adequately handle reads crossing over one or
more splice junctions (TopHat2, STAR, and GSNAP), the alignment rate of L262
was consistently higher than that of L75, suggesting that longer reads can
improve overall alignment rate. The difference was most dramatic for GSNAP,
for which the rate was higher by 10.39 percentage points (pp) for L262 than
for L75. In contrast, Bowtie2 suffered greatly for L262, achieving a rate
10.07 pp lower than that of L75. This implies that the appropriate handling of
splice junctions is a crucial requirement for a mapping tool to be effectively
applied to libraries with longer read lengths, likely because longer reads
have a higher chance of crossing over one or more splice junctions
(Supplementary Figure 6). To illustrate, in the GSNAP output, 50.18% of
uniquely mapped reads in L262 contain at least one splice junction, while the
same is true for only 38.24% of uniquely mapped reads in L75. Higher alignment
rate of longer reads, however, comes at a computational cost: the execution
time of splice junction-aware mapping tools was much longer for L262 than for
L75 (Figure 1b). Most notably, TopHat2’s average runtime per read was
approximately 58 times longer for L262, suggesting that TopHat2’s underlying
alignment algorithm is more sensitive to read length than the other tools we
tested. STAR was the fastest tool for both libraries, but it used a
considerably larger amount of memory than others (data not shown). These
results demonstrate varying degrees of scalability of currently available
read-mapping tools.

### RNA and DNA genotype concordance

The higher proportion of mapped reads observed in L262 could in theory reflect
an undesirable increase in uncertain, _incorrectly_ mapped reads. To address
this possibility, we evaluate evidence of read mismapping as indicated by
genotype calls from RNA-seq reads. We sought to compare the effect of read
length on RNA/DNA genotype concordance by analyzing separately the concordance
of SNPs and indels. For SNPs, using L75 and L262 mappings with STAR, GSNAP,
and TopHat2, we calculated a _discordance $p$-value_ at each genotyped site,
which represents the level of deviation from the expected genotype due to
inaccurate mapping, while controlling for dataset-specific sequencing error
rate (Methods). The distribution of discordance $p$-values was significantly
higher (more concordant) in L262 across all read depth bins for GSNAP and in
all but one bin for STAR (Figure 2). This suggests that, in addition to an
increase in the _proportion_ of uniquely aligned reads, the longer read
protocol resulted in a lower rate of mapping errors. For TopHat2, the
improvement in concordance for long reads was not as clear, but we suspect
this is in part due to the lack of statistical power caused by low alignment
rate. Next, we sought to assess whether read length influences genotype
concordance estimates for indels of varying size (Supplementary Figure 8). In
general, we did not observe major effects of read length, and STAR and GSNAP
perform reasonably well in mapping indels as large as 20 bp in length. STAR
appears to perform equally well in mapping insertions and deletions,
performing best with the long reads whereas GSNAP appears to have little gain
and an overall bias in mapping deletions. TopHat2, on the other hand, is very
limited in its ability to call indels, maxing out in its ability to call an
indel at 2 bp in length. We hypothesize that the lack of improvement of indel
genotype concordance for longer reads across all methods indicates that
current limitations are predominantly algorithmic. Based on the mapping rate,
the GSNAP output was used for the remainder of the study.

### Quantification of gene and isoform abundance

To analyze the effect of read length on the quantification of gene abundance,
we calculated the number of reads that unambiguously mapped to each gene
(Methods) and compared the relative abundance between L262 and L75. Overall,
the normalized read counts were highly correlated between the two libraries (r
= 0.94), but we found an excess of genes that are more represented in L262
(Figure 3a). For instance, among 13,969 genes observed in both libraries (with
a fraction of reads $>10^{-6}$), 139 genes displayed more than tenfold
increase in L262, whereas only one gene (a pseudogene named RP5-857K21.11)
showed more than tenfold decrease in L262. Among the 139 genes that were more
represented in L262, 130 were pseudogenes, which suggested that longer reads
are more effective at capturing regions of high sequence similarity with other
regions in the genome. In fact, 3.64% of uniquely mapped reads in L262 are
from pseudogenes, which is considerably larger than 1.25% for L75
(Supplementary Figure 7). Among the few protein coding genes that displayed
more than tenfold increase in L262 were GLUD2 and NACA2, which are known to be
highly homologous to GLUD1 and NACA, respectively. To further test whether
long reads are more effective at capturing non-unique regions, we used the
mappability score obtained from the UCSC Genome Browser [26] as a measure of
sequence uniqueness. Sorting the genes according to average mappability score
revealed that genes with low mappability score tend to be more represented in
L262 than in L75 (Figure 3c). The same trend could be observed when we
restricted our analysis to protein coding genes. Leaving out pseudogenes and
genes with an average mappability score of less than one left 25,597 genes
(out of 56,680), and the correlation between normalized read counts of the two
libraries for this subset was considerably higher (r = 0.97, Figure 3b).

We next asked how much more effective long reads are at disambiguating mRNA
isoforms. For each read in the library, we calculated the number of isoforms
that can be assigned to any read (Methods) and compared its distribution
between the two libraries. Overall, the reads in L262 were assigned to fewer
different isoforms than those in L75 (Figure 4a). In particular, 32.98% of
reads in L262 that were assigned to at least one isoform were unambiguously
assigned to a single isoform, which is 7.94% higher than 25.04% observed for
L75. This agrees with intuition that longer reads have increased likelihood of
crossing a splice junction or containing informative bases that allow
disambiguation among candidate isoforms. Furthermore, the same pattern could
be observed within a library; by stratifying the reads in L262 according to
their _effective alignment length_ (the number of bases in the reference
covered by an aligned read pair; distribution shown in Supplementary Figure
2), we observed that reads with longer alignments tend to have fewer matching
isoforms indicating better specificity of the originating transcript (Figure
4b). In fact, if we only take the subset of 104,910 read pairs in L262 with
effective alignment length of 524 bp (the longest possible), the unambiguous
mapping rate increases to 40.96%.

As longer read data provides the opportunity to test the efficacy of
transcript abundance methods, we analyzed the relative ability of both
Cufflinks and FluxCapacitor to infer transcript abundances from L75 data.
Here, we used the L75 data mapped with GSNAP and ran both methods. We then
compared the correlation of estimated transcript abundances to those directly
measured in the L262 data. Transcript abundances from L262 were obtained by
identifying reads which unambiguously mapped to single transcripts,
summarizing counts for each transcript, and further dividing each by the
number of bases that are unique to that transcript. The latter provides a
means to account for the fact that counts are biased by the effective
proportion of the transcript length that is unique. Using this straightforward
approach to evaluate this class of tools, we identified better correlation
with Cufflinks ($R^{2}=0.49$) than FluxCapacitor ($R^{2}=0.34$; Figure 5).

### Detection of allele-specific expression patterns

Studies of allele-specific expression from RNA-seq data require
differentiating between alleles, which is possible for reads over one or more
heterozygous loci. Given full haplotype information for GM12878, we evaluated
the fractions for which the parental haplotype in unambiguous (Methods). We
were able to assign 7.01% of the reads in L262 to either maternal or paternal
haplotype as opposed to 4.26% for L75 (Figure 6a). This increased rate of
haplotype differentiation for longer reads could also be seen within L262 when
the reads were stratified by effective alignment length (Figure 6b). Noting
that the number of heterozygous sites in a given interval approximately
follows a Poisson distribution, we fit a Poisson model to calculate the
expected proportion of allele-specific reads for libraries with a read length
greater than 262 (Figure 6c). Our fitted model predicts that approximately a
quarter of the reads that cover 1,000 bases in the reference will contain at
least one heterozygous site.

Next, to evaluate the extent to which reference allele bias in read-mapping
affects the detection of allele-specific patterns in our data, we ran the
following simulation. Given the reads we classified as either paternal or
maternal, we flipped the bases at heterozygous sites. We then remapped the
modified reads and calculated, for each gene, the proportion of reads that
mapped to the same location as before, which we refer to as mapping retention
rate. When we ran this procedure for both L262 and L75, we found that L262
leads to significantly higher mapping retention rate overall (Supplementary
Figure 9). The fraction of genes with retention rate lower than 0.9 was 19.27%
for L75, while it was only 7.09% for L262. This suggests that allelic
imbalance estimates measured by long reads are more reliable than those from
short reads, and specifically less subject to false positives arising from
mapping artifacts. For the remaining parts of this study, we disregarded genes
with mapping retention rate lower than 0.9 for each library to avoid calling
allelic effects as significant that may only reflect mapping errors.

Because longer read length leads to more reads with heterozygous sites and
fewer genes marked as untestable by the mapping retention rate filter, L262 is
expected to be statistically more powerful in discovering allele-specific
expression patterns. Indeed, when we checked the number of significant cases
of gene-level ASE (Methods) at FDR = 10%, we got a considerably larger number
with L262 than L75, where the latter was randomly subsampled in two ways for
comparison by matching either the number of bases or the number of reads
(Figure 7a). We also observed fewer significant cases (180 at FDR = 10%) in a
simulated short-read library constructed by truncating each read in L262 to 75
bp. Interestingly, 206 out of 381 significant cases found with L262 were not
called significant by any of the subsampled L75 libraries, and among those, 24
did not pass the mapping retention rate filter in L75. One potential
explanation for cases that are only detected in L262 is related to
mappability: the mean mappability scores of genes with significant ASE in L262
that were not detected in L75 were significantly lower than that of genes
called significant in both (Wilcoxon rank-sum $p$-value: $3.70\times
10^{-4}$). To further validate significant ASE called with L262, we cross-
referenced our significant results with a list of known eQTLs from a large
population based study [4]. 46 out of 3,328 tested genes had a heterozygous
site at a known eQTL SNP for the gene, and 12 of those cases appeared in the
final set of 381 significant ASEs, showing a considerable enrichment
($\chi^{2}$ $p$-value = $0.018$). These results are further confirmed using a
standard, per-locus ASE test (as in [4]). In this case, we also observed a
higher rate of ASE detection in L262 (Supplementary Figure 10).

We next looked for a more complex form of allelic effects, allele-specific
alternative splicing, by assessing haplotype-specific differences in isoform
abundance with the exon inclusion-exclusion test developed in [9] (Methods),
in which the ratio of reads that support the inclusion or exclusion of an exon
block is compared between the two haplotypes. At FDR = 10%, we obtained 69
significant cases of differential exon usage with L262, which supported ASAS
of 38 distinct genes. Subsampled libraries of L75 detected significantly fewer
cases in general (Figure 7), resulting in a median of 41 differentially used
exon blocks (24 genes) for the base-subsampled L75 and 6.5 exon blocks (4.5
genes) for the read-subsampled L75. A simulated short-read library constructed
by truncating each read in L262 to 75 bp also resulted in fewer discoveries
(37 exon blocks, representing 25 genes). The coverage plots of the top two
cases detected from L262 – SNHG5 and IFI44L – with raw $p$-value $<10^{-12}$
show clear haplotype-specific difference in distribution of reads across exons
of the genes (Figure 8). Among other significant findings, ZNF83 (raw
$p$-value = $1.44\times 10^{-4}$) and CASP3 (raw $p$-value = $3.22\times
10^{-4}$) are supported by prior knowledge; they contain a differentially used
exon block that is known to be affected by an sQTL [27] where this individual
has a heterozygous genotype (Figure 8). CASP3 was not called significant in
any of the subsampled L75 libraries. These examples of ASAS represent likely
cases of a cis-acting genetic variant altering splicing in this individual.

### Sequencing cost comparison

While we have shown the degree to which longer read length provides a number
of key benefits for transcriptome analysis, it is important to note that it is
generally more expensive. Based on actual costs, we estimated the sequencing
cost per million bases of L262 to be $$0.176$ in contrast to $$0.103$ of L75
(Methods), which implies that one 262 bp paired-end read costs approximately
the same as six 75 bp paired-end reads. The sheer number of additional short
reads one can obtain in lieu of increasing read length renders some of the
benefits of longer reads cost-ineffective. In particular, the number of reads
unambiguously mapped to a genomic location, an mRNA isoform, or a parental
haplotype would all be significantly higher if we were to use a short read
dataset that is six times larger than the cost-matched long read dataset. In
terms of the ability to detect allele-specific patterns, L262 was
approximately $364 more expensive than the base-subsampled L75, which resulted
in 24% fewer cases in ASE detection and 44% fewer in ASAS detection on
average. When we matched the costs between the two libraries by using all
reads in L75 and 45.45% of the reads in L262, we found that L262 actually
detects 60% fewer ASE cases and 88% fewer ASAS cases on average than L75,
which detects 418 ASE and 106 ASAS cases (Supplementary Figure 11).
Nevertheless, among the cases detected by the cost-matched L262, 9% of ASE
cases and 11% of ASAS cases were not detected by L75 because they did not pass
the mapping retention rate filter. This implies that long reads, even with
lower total power than short reads at the same cost, more accurately capture
cases that are susceptible to mapping errors. Clearly, these results depend on
the current sequencing costs; as the price of long reads becomes more
comparable to that of short reads, we expect the benefits of longer read
lengths to become more apparent. In addition, using longer reads fundamentally
reduces the mapping bias and leads to more accurate quantification of genomic
elements, which cannot be achieved by simply obtaining a larger short read
dataset.

## Discussion

Longer read lengths in RNA-sequencing offer a number of advantages ultimately
supporting more accurate quantification of diverse transcriptional effects, as
demonstrated here through a comparative analysis of two libraries from the
same cell line with different read lengths. First, the unique alignment rate
of long reads is observed to be higher, with GSNAP uniquely mapping 93.44% of
long reads in contrast to 83.05% of short reads. Intuitively, longer sequences
are less likely to occur more than once in the genome, even with the
prevalence of regions sharing sequence similarity, such as gene families and
pseudogenes. In addition to the higher mapping rate, our analysis of genotype
concordance suggests that longer reads are significantly less affected by
mapping biases.

Next, there are differences evident in the estimates of gene expression
quantified from the two libraries, which in principle could be influenced by
differences in read mapping. In our study, genes with evidence of differential
expression by read length were primarily genes with low mappability, and the
proportion of long reads that unambiguously mapped to a pseudogene was roughly
three times higher than that of short reads, with a similar observation in
protein coding genes with low uniqueness. This suggests that the improved
mapping of longer reads also provides more accurate estimates of gene
expression. Furthermore, it reinforces that mappability is an important
parameter to consider when comparing data between datasets with varying read
lengths (such as for differential expression analyses).

Long reads can also reduce the ambiguity between different transcripts arising
from the same gene. Because such transcripts share large fractions of their
sequence in common, transcript-level quantification is a considerably harder
task than gene-level quantification. Only reads which span a splice junction
or which include bases simply not found in a particular transcript (such as an
exon excluded from one isoform) can directly disambiguate between related
transcripts. Here, we demonstrated that the longer read protocol indeed
allowed improved disambiguation; the proportion of exonic reads unambiguously
assigned to a single isoform was 32.98% for L262 and 25.04% for L75. We
further provide data and a straightforward means for comparing transcript
abundance methods by comparing their inferences from L75 to directly measured
abundances in L262; using this approach we have identified a better
correlation for Cufflinks compared to FluxCapacitor.

The use of RNA-sequencing to detect allele-specific transcriptional effects
offers the potential to pinpoint cis-regulatory effects [4, 28], but relies on
identifying which allele (or haplotype) each read arises from. Long reads not
only increase the proportion of reads unambiguously assigned to each haplotype
(7% for L262 compared to 4% for L75), but also diminish the susceptibility of
ASE analyses to mapping biases as shown by our mapping retention rate
analysis. Thus, longer reads provide a more robust quantification of allele-
specific expression. We observed that the longer read dataset actually detects
more significant cases than the short read dataset with the same total number
of bases, and identified two particularly compelling cases of ASAS supported
by a known sQTL, one of which could not be detected with short reads.
Detecting such allelic effects could dramatically improve our understanding of
the cis-regulation of alternative splicing within an individual. However,
initial efforts at quantifying allele-specific alternative splicing have
largely relied on indirect signals such as apparent differences in insert size
[4], or have been restricted to a small number loci with high coverage near
splice junctions [29]. Longer reads, which are both more likely to contain
heterozygous loci and more likely to disambiguate among related transcripts
and isoforms, thus offer a particularly large advantage in quantifying ASAS
and identifying significant cis-regulated alternative splicing effects.

Overall, the increasing availability of long read RNA-seq libraries offers
significant advantages in analysis of the transcriptome. While the current
long read sequencing technology is considerably more expensive than short read
sequencing, more accurate quantification of various aspects of transcriptome
due to lower mapping bias can only be achieved with longer read lengths. Where
identification of alternative splicing, fusion events, and novel transcripts
have traditionally required sophisticated computational methods and
statistical inference, longer reads will bring us increasingly close to direct
measurement of these events. Further distinguishing the benefits, long reads
allow more accurate SNP genotype calls from RNA-seq data and, critically,
increasingly detailed analysis of allele-specific expression. As public data
increasingly includes RNA-seq data of different read lengths, understanding
these advantages and caveats will lead to better tools and comparisons of
heterogeneous data.

## Acknowledgments

This work was supported by the Edward Mallinckrodt Jr. Foundation and the
National Institutes of Health (R01 MH10814).

## References

  * 1. Plessy C, Bertin N, Takahashi H, Simone R, Salimullah M, et al. (2010) Linking promoters to functional transcripts in small samples with nanoCAGE and CAGEscan.  Nature Methods 7: 528–534.
  * 2. Ni T, Corcoran DL, Rach EA, Song S, Spana EP, et al. (2010) A paired-end sequencing strategy to map the complex landscape of transcription initiation.  Nature Methods 7: 521–527.
  * 3. Valen E, Pascarella G, Chalk A, Maeda N, Kojima M, et al. (2009) Genome-wide detection and analysis of hippocampus core promoters using DeepCAGE.  Genome Research 19: 255–265.
  * 4. Montgomery SB, Sammeth M, Gutierrez-Arcelus M, Lach RP, Ingle C, et al. (2010) Transcriptome genetics using second generation sequencing in a Caucasian population.  Nature 464: 773–777.
  * 5. Sultan M, Schulz MH, Richard H, Magen A, Klingenhoff A, et al. (2008) A global view of gene activity and alternative splicing by deep sequencing of the human transcriptome.  Science 321: 956–60.
  * 6. Levin JZ, Berger MF, Adiconis X, Rogov P, Melnikov A, et al. (2009) Targeted next-generation sequencing of a cancer transcriptome enhances detection of sequence variants and novel fusion transcripts.  Genome Biology 10: R115.
  * 7. Peng Z, Cheng Y, Tan BCM, Kang L, Tian Z, et al. (2012) Comprehensive analysis of RNA-Seq data reveals extensive RNA editing in a human transcriptome.  Nature Biotechnology 30: 1–10.
  * 8. Cloonan N, Forrest ARR, Kolle G, Gardiner BBA, Faulkner GJ, et al. (2008) Stem cell transcriptome profiling via massive-scale mRNA sequencing.  Nature Methods 5: 613–619.
  * 9. Salzman J, Jiang H, Wong WH (2011) Statistical Modeling of RNA-Seq Data.  Statistical Science 26: 62–83.
  * 10. Shendure J, Ji H (2008) Next-generation DNA sequencing.  Nature Biotechnology 26: 1135–1145.
  * 11. Ozsolak F, Milos PM (2011) RNA sequencing: advances, challenges and opportunities.  Nature Reviews Genetics 12: 87–98.
  * 12. Benjamini Y, Speed TP (2012) Summarizing and correcting the GC content bias in high-throughput sequencing.  Nucleic Acids Research 40: 1–14.
  * 13. Oshlack A, Wakefield MJ (2009) Transcript length bias in RNA-seq data confounds systems biology.  Biology Direct 4: 14.
  * 14. Marioni JC, Mason CE, Mane SM, Stephens M, Gilad Y (2008) Rna-seq: An assessment of technical reproducibility and comparison with gene expression arrays.  Genome Research 18: 1509-1517.
  * 15. Schwartz S, Oren R, Ast G (2011) Detection and Removal of Biases in the Analysis of Next-Generation Sequencing Reads.  PLoS ONE 6: 12.
  * 16. Trapnell C, Williams BA, Pertea G, Mortazavi A, Kwan G, et al. (2010) Transcript assembly and quantification by RNA-Seq reveals unannotated transcripts and isoform switching during cell differentiation.  Nature Biotechnology 28: 511–515.
  * 17. Nariai N, Hirose O, Kojima K, Nagasaki M (2013) TIGAR: transcript isoform abundance estimation method with gapped alignment of RNA-Seq data by variational Bayesian inference.  Bioinformatics Oxford England .
  * 18. Li B, Dewey CN (2011) Rsem: accurate transcript quantification from rna-seq data with or without a reference genome.  BMC bioinformatics 12: 323.
  * 19. Kim D, Pertea G, Trapnell C, Pimentel H, Kelley R, et al. (2013) TopHat2: accurate alignment of transcriptomes in the presence of insertions, deletions and gene fusions.  Genome Biology 14: R36.
  * 20. Wu TD, Nacu S (2010) Fast and SNP-tolerant detection of complex variants and splicing in short reads.  Bioinformatics Oxford England 26: 873–81.
  * 21. Langmead B, Salzberg SL (2012) Fast gapped-read alignment with Bowtie 2.  Nature Methods 9: 357–360.
  * 22. Dobin A, Davis CA, Schlesinger F, Drenkow J, Zaleski C, et al. (2012) STAR: ultrafast universal RNA-seq aligner.  Bioinformatics : bts635–.
  * 23. Harrow J, Frankish A, Gonzalez JM, Tapanari E, Diekhans M, et al. (2012) GENCODE: The reference human genome annotation for The ENCODE Project.  Genome Research 22: 1760–1774.
  * 24. Li X, Yin X, Li J (2010) Efficient identification of identical-by-descent status in pedigrees with many untyped individuals.  Bioinformatics 26: i191–i198.
  * 25. Benjamini Y, Hochberg Y (1995) Controlling the false discovery rate: a practical and powerful approach to multiple testing.  Journal of the Royal Statistical Society Series B (Methodological) : 289–300.
  * 26. Meyer LR, Zweig AS, Hinrichs AS, Karolchik D, Kuhn RM, et al. (2013) The UCSC Genome Browser database: extensions and updates 2013.  Nucleic acids research 41: D64–9.
  * 27. Battle A, Mostafavi S, Zhu X, Potash JB, Weissman MM, et al. (2013) Characterizing the genetic basis of transcriptome diversity through rna-sequencing of 922 individuals.  Genome research .
  * 28. Pickrell JK, Marioni JC, Pai AA, Degner JF, Engelhardt BE, et al. (2010) Understanding mechanisms underlying human gene expression variation with RNA sequencing.  Nature 464: 768–72.
  * 29. Lappalainen T, Sammeth M, Friedländer MR, Hoen PA, Monlong J, et al. (2013) Transcriptome and genome sequencing uncovers functional variation in humans.  Nature 501: 506–511.
  * 30. Thorvaldsdóttir H, Robinson JT, Mesirov JP (2013) Integrative Genomics Viewer (IGV): high-performance genomics data visualization and exploration.  Briefings in bioinformatics 14: 178–92.

## Figure Legends

Cho, H. et al

Figure 1:  The effect of read length on read-mapping performance. We compared
the percentage of reads uniquely mapped (a) and the average runtime per a
million reads (b) of Bowtie2, TopHat2, STAR, and GSNAP on L262 and L75. Only
the splice-mappers – TopHat2, STAR, and GSNAP – achieved higher alignment
rates for L262 compared to L75. The increase in runtime going from L75 to L262
varied greatly across the mappers, TopHat2 being the most sensitive among the
four mappers tested.  Figure 2:  The effect of read length on RNA genotype
concordance. Results of comparisons of discordance p-values (-log10
transformed) between L262 and base-subsampled L75 stratified by mapper and
read depth. Sites were stratified into four groups based on read depth in each
dataset: 1) 10-99, 2) 100-499, 3) 500-999, and 4) $\geq$1000 reads. For each
mapper, at each of the four read depths, a one-sided Wilcoxon rank-sum test
was conducted to determine whether the $p$-values for the L262 dataset were
significantly higher compared to L75. The red horizontal line represents a
significance level of $\alpha=0.05$. Each symbol represents an independently
subsampled data for L75.  Figure 3:  Most of the difference in gene
quantification arise from poorly mappable genes. (a, b) Log scatter plots of
fraction of reads mapped to each gene between L262 and L75. Only plotting
genes that are not pseudogenes and have perfect mappability scores leads to a
scatter plot with near perfect correlation. (c) Genes that are observed in
both libraries and are at least 500 bp long were divided into five groups
according to the mean mappability score, which is obtained by summarizing the
75 bp mappability score track on UCSC Genome Browser. For each group, we
computed the fraction of genes that had a higher proportion of reads in L262
than in L75. The groups corresponding to lower mappability scores showed
higher proportion of genes more represented in L262. The same trend was
observed even when we restricted our analysis to protein coding genes.  Figure
4:  Longer reads are consistent with fewer number of mRNA isoforms. If all the
bases of a read that maps to a gene are part of the same isoform, this common
isoform (which does not have to be unique) is said to be consistent with the
read. We compared the cumulative distribution of the number of consistent
isoforms of reads in L75 versus reads in L262 (a) and across read groups
within L262 (b) divided according to the number of reference bases covered by
the alignment.  Figure 5:  Comparison of Cufflinks and FluxCapacitor on
transcript quantification. We used transcript counts from L262 to assess
performance of (a) Cufflinks and (b) FluxCapacitor on L75. Cufflinks and
FluxCapacitor were respectively applied on base-subsampled L75 to generate
FPKM/RPKM counts. On L262, we only counted reads that unambiguously mapped to
each transcript, and divided each by the number of bases that are unique to
that transcript. We restricted our analysis to 41,715 protein coding
transcripts with at least 10 unique bases. Using this approach, we identified
that Cufflinks has significantly better correlation to the long-read data than
FluxCapacitor ($p$-value $<10^{-163}$).  Figure 6:  Longer reads have a higher
chance of containing a heterozygous site that allows us to differentiate the
haplotype of origin.  We compared the proportion of reads unambiguously
assigned to maternal or paternal haplotypes between L75 and L262 (a) and
within L262 (b), where reads in L262 were grouped by the number of reference
bases covered by the alignment of (possibly overlapping) paired-end reads. The
observed data in (b) were extrapolated using a poisson model to even longer
read lengths. To provide additional support for the extrapolation, we directly
approximated the proportion of reads containing at least one heterozygous site
based on the GENCODE annotation and the genotype of this individual, assuming
that the true gene expression levels are the same as those observed in L262,
that all mRNA isoforms of a gene are uniformly expressed, and that each
starting location of a transcript is equally likely to be included in the
library. We only considered reads that mapped to autosomal chromosomes. Our
approximation agreed well with both the observed data and the extrapolation.
Figure 7:  Longer reads enable more effective detection of allele-specific
patterns. To make a fair comparison, L75 was randomly subsampled ($n=10$) down
to the same number of bases (L75B) and the same number of reads (L75R) as
L262. We ran the subsampled libraries and L262 through the pipeline for
detecting allele-specific gene expression and allele-specific exon inclusion-
exclusion patterns. The figures above show the number of significant cases we
called at FDR = 10% using the Benjamini and Hochberg procedure [25].  Figure
8:  Examples of allele-specific isoform distribution detected with L262. Shown
are the IGV [30] coverage plots of (a) genes corresponding to the top two
significant cases (raw $p$-value $<10^{-12}$) and (b) two significant cases
(FDR = 10%) with a strong prior support by having a heterozygous site at a
known sQTL from a large population-based study. The implicated sQTL and the
exon that is known to be affected are noted by a red line and a black box,
respectively. Red asterisks above each plot denote blocks with significant
differential inclusion-exclusion pattern (FDR = 10%).

##  High-resolution transcriptome analysis with long-read RNA sequencing
Supplementary Material

Cho, H. et al

Supplementary Figure 1:  Estimated fragment length distribution. We used the
reads that mapped to genes with exactly one transcript to estimate the
fragment length distributions for L262 and L75. The read alignments produced
by GSNAP were used.  Supplementary Figure 2:  Histogram of number of reference
bases covered by each read. As a measure of information content, we counted
the number of bases in the reference genome covered by each read pair. For
overlapping reads, for instance, this will be less than double the read
length. From the histograms, we can see that the vast majority of reads in L75
are non-overlapping while reads in L262 show high degrees of overlap. Overall,
however, L262 reads still have more information content than L75 reads.
Supplementary Figure 3:  Estimated per-base sequencing error rate for each
aligned dataset. Sequencing error rate is estimated as the total number of
mismatches in the alignments divided by the total number of aligned bases.
Supplementary Figure 4: Comparison of read depth distributions. Genotyped
sites were stratified into the same four read depth bins as used in the
genotype concordance analysis. This plot compares the distribution of read
depths within each category across mappers and read length. The distributions
within each read depth category are highly similar across read lengths.
Supplementary Figure 5: Illustration of differential exon usage test. We
tested whether each exon block exhibits a differential usage pattern using a
chi-square test on a 2-by-2 contingency table where the rows indicate whether
a read is maternal or paternal and the columns indicate whether a read
overlaps with the exon block of interest. Each example read (shown as single-
end for simplicity) is marked with the cell in the contingency table it is
counted towards.  Supplementary Figure 6:  Distribution of number of splice
junctions spanned by each read.  For each read pair, we counted the number of
splice events represented in the alignment from GSNAP. The bar graph shows the
distribution of this number in both L262 and L75. We can see that L262 reads
tend to span more splice junctions than L75.  Supplementary Figure 7:
Breakdown of reads by gene category. The fraction of reads mapped to each of
the five major gene categories in GENCODE annotation (version 15) is shown.
L262 displayed a relatively high proportion of reads originating from
pseudogenes compared to L75.  Supplementary Figure 8:  Expected indel size
from DNA-seq compared to most frequently observed size from RNA-seq for L262
and L75. Points along the diagonal (gray line) in each plot indicate a
positive correlation between most frequently observed indel size and expected
size based on genotyping. In general, STAR and GSNAP perform reasonably well
in mapping indels as large as 20 bp in length. Tophat, on the other hand, is
very limited in its ability to call indels, maxing out in its ability to call
an indel at 2 bp in length. STAR appears to perform equally well in mapping
insertions and deletions, performing best with the long reads. GSNAP appears
to have a bias in mapping deletions. The L75 dataset was subsampled to match
the number of bases in L262 for this analysis.  Supplementary Figure 9:  Genes
quantified by longer reads display higher mapping retention rate. Mapping
retention rate is a simulation-based metric that is inversely correlated with
the severity of allelic bias in allele-specific quantification. It is computed
for each gene as the proportion of reads that map to the same location after
flipping all bases that align with heterozygous sites in the reference so that
the new bases are the alleles from the other parental haplotype. The above
histogram compares the distribution of mapping retention rates in L75 versus
L262.  Supplementary Figure 10:  Comparison of read mapping tools and read
length for ASE detection. A. Presents the total number of heterozygous sites
tested per mapper (Bowtie, STAR, GSNAP, and Tophat) and read length (75bp and
262bp) combination. The 75bp read length data was subsampled to mirror the
number of bases (red) and reads (green) in the 262bp data. The number of
tested sites is very similar between STAR and GSNAP regardless of read length,
while the number of tested sites for Tophat and Bowtie is significantly lower
(on the order of a few thousand sites). This result is in keeping with the
number of mapped reads from STAR and GSNAP versus Bowtie and Tophat. B.
Percent of significant ASE sites (pvalue ¡ .05) for each mapper and read
length pair. In general, the long read samples, regardless of mapper, showed a
higher percentage of significant sites.  Supplementary Figure 11:  Comparison
of allele-specific pattern detection with matching monetary costs. L262 was
subsampled down to 45.45% of the original data to match its monetary cost with
that of L75. This severely reduces the amount of information contained in the
subsampled L262 library, and thus results in a significantly fewer discoveries
compared to L75.


# A Note on quasi-transitive graphs
quasi-isometric to planar (Cayley) graphs

Joseph Paul MacManus  Mathematical Institute, University of Oxford, Oxford,
OX2 6GG, UK<EMAIL_ADDRESS>
(Date: July 18, 2024)

###### Abstract.

Given a connected, locally finite, quasi-transitive graph $X$ which is quasi-
isometric to a planar graph $\Gamma$, we remark that one can upgrade $\Gamma$
to be a planar Cayley graph, answering a question raised by Esperet–Giocanti
and Hamann.

## Introduction

The topic of ‘coarse graph theory’ is a fast-growing area, where the goal is
to understand how ‘large-scale’ assumptions affect the geometry of graphs.
This point-of-view in graph theory was recently popularised by Georgakopoulos
and Papasoglu, and we direct the interested reader to their article for
further reading [8].

One particular point of interest has been in understanding those graphs which
are quasi-isometric to planar graphs. For example, it was shown by the author
in [14] that if a quasi-transitive graph is quasi-isometric to a planar graph
then it is necessarily accessible, in the sense of Thomassen–Woess [23]. By
combining this statement with a deep theorem originating in the work of Mess
[17], it is further shown that if a Cayley graph is quasi-isometric to a
planar graph, then it is quasi-isometric to a planar Cayley graph.

Recently, it was shown by Esperet and Giocanti that quasi-transitive graphs
which are ‘minor-excluded’ are quasi-isometric to planar graphs [5]. There,
they raise the question as to whether one can take the aforementioned planar
graph to be itself quasi-transitive, or possibly to be a Cayley graph. In [12]
Hamann shows that one can upgrade the planar graph in their theorem to indeed
be quasi-transitive, but raises the more general question of if given any
quasi-transitive graph which is quasi-isometric to a planar graph $\Gamma$,
can we similarly upgrade $\Gamma$? In this short note, we settle the above and
record the following fact.

alphtheoremmain Let $X$ be a connected, locally finite, quasi-transitive
graph. Suppose $X$ is quasi-isometric to a planar graph. Then $X$ is quasi-
isometric to a planar Cayley graph.

More precisely, such a graph is quasi-isometric to either $\\{1\\}$,
$\mathbf{Z}$, $\mathbf{Z}^{2}$, a cocompact Fuchsian group, or some free
product of the aforementioned groups. It is also possible to further upgrade
the planar Cayley graph in Theorem Introduction to be that of a Kleinian
function group. See Remark 2.2 below for more details.

It was a long-standing question of Woess [26, Prob. 1] as to whether every
vertex-transitive graph is quasi-isometric to a Cayley graph. A counterexample
was eventually presented in [4]. Somehow, the geometry of our setting means
that a similar counterexample is impossible here.

The proof of Theorem Introduction proceeds case-wise. First considering one-
ended graphs, we separately consider those which are hyperbolic and those
which are not, in the sense of Gromov. In the hyperbolic case, we combine the
convergence theorem of Tukia, Gabai, and Casson–Jungreis with a corollary of
the Montgomery–Zippin theorem to deduce the desired conclusion. In the non-
hyperbolic case, we appeal to an inequality of Varopoulos to deduce that our
graph has a quadratic growth function. Following this, a theorem of Woess [27]
implies that our graph is quasi-isometric to $\mathbf{Z}^{2}$. The infinite-
ended case is then deduced by combining a theorem of Hamann [11]
characterising quasi-isometries between accessible graphs, together with a
theorem of the author stating that the graphs in question are indeed
accessible [14].

### Conventions

Every graph in this note is assumed to be simplicial. We take $0$ to be a
natural number.

##  1\. One-ended graphs

Throughout this section, let $X$ be a connected, locally finite, quasi-
transitive, one-ended graph which is quasi-isometric to some planar graph. We
will assume a working knowledge of the theory of isometry groups of graphs,
viewed as topological groups. We point the reader to [3] as a good modern
reference. The survey article of Woess [26] also serves as a good introductory
text.

###  1.1. Hyperbolic graphs

In this subsection we assume that $X$ is hyperbolic, in the sense of Gromov.
To avoid confusion we will call a metric space Gromov-hyperbolic, while a
finitely generated group shall be said to be word-hyperbolic. That is, a
finitely generated group is word-hyperbolic if it admits a Gromov-hyperbolic
Cayley graph.

By [14, Thm. B], we have that $X$ is quasi-isometric to a complete Riemannian
plane. It is a standard fact that the boundary of a Gromov-hyperbolic complete
Riemannian plane is a circle. Let $G=\operatorname{Isom}(X)$ denote the
isometry group of $X$, which acts faithfully and cocompactly on $X$. We have
that $G$ is a totally disconnected, locally compact topological group.
Consider the kernel $N\leq G$ of the action of $G$ on $\partial_{\infty}X$.
This is a compact normal subgroup of $G$, and $H:=G/N$ acts faithfully and
continuously on the circle. Every quotient of a locally compact group is
locally compact, so $H$ (or an index-two subgroup of $H$) is a locally compact
subgroup of $\operatorname{Homeo}_{+}(S^{1})$. By a corollary of the
Montgomery-Zippin theorem [9, Thm. 4.7], it follows that $H$ is a Lie group.
Since $H$ is totally disconnected, we deduce that it is discrete.

Let $Y=X/N$ be the quotient graph. Note that since $N$ is a compact subgroup
of $\operatorname{Isom}(X)$ we have that $N$ acts on $X$ with bounded orbits,
and so the quotient map $\pi:X\to Y$ is a quasi-isometry. It is easy to check
that $H$ acts faithfully and continuously on $Y$. It follows that the quotient
group $H=G/N$ homeomorphically embeds as a topological group into
$\operatorname{Isom}(Y)$ [3, Prop. 5.B.6]. Since $H$ is discrete, we deduce
that the action of $H$ on $Y$ is properly discontinuous [3, Prop. 5.B.10].
Since $Y/H$ is compact, so we see by the Švarc–Milnor lemma that $H$ is quasi-
isometric to $Y$ (and thus also to $X$). Therefore $H$ is a word-hyperbolic
group with a circular boundary. By the convergence group theorem of Tukia,
Gabai, and Casson–Jungreis [24, 6, 2] we deduce that $H$ is virtually
Fuchsian. We have thus shown the following.

######  Theorem 1.1.

Let $X$ be a connected, locally finite, quasi-transitive graph. Suppose $X$ is
one-ended, Gromov-hyperbolic, and quasi-isometric to a planar graph. Then $X$
is quasi-isometric to a Fuchsian group, and thus to a planar Cayley graph.

###  1.2. Non-hyperbolic graphs

Now, we assume that $X$ is not hyperbolic. In this case it will be helpful to
assume without loss of generality that $X$ is vertex-transitive. This is
harmless, as it is a standard fact that every connected, locally finite,
quasi-transitive graph is quasi-isometric to a connected, vertex-transitive,
locally finite graph [27, Prop. 5.1].

Proceeding in the spirit of an argument presented in [19, § 3], we will make
use of an isoperimetric inequality of Varopoulos, as it appears in [22, Thm.
2.1].

Given our transitive graph $X$, let $V(n)$ denote the size of the closed
$n$-ball around a given vertex. In other words, $V:\mathbf{N}\to\mathbf{N}$ is
the growth function of $X$. Given a subset $\Omega\subset VX$, we denote by
$\partial\Omega$ the set of vertices which lie at a distance of exactly 1 from
$\Omega$.

######  Proposition 1.2.

The vertex-transitive graph $X$ satisfies the following isoperimetric
inequality. For all finite $\Omega\subset VX$, we have that

$\frac{|\Omega|}{\phi(2|\Omega|)}\leq 4|\partial\Omega|,$

where $\phi(\lambda)=\inf\\{n\in\mathbf{N}:V(n)>\lambda\\}$.

The statement above has the easy consequence that if there exists
$C,D\in\mathbf{N}$ such that $V(n)\geq Cn^{D}$ for all $n\geq 1$, then there
exists $C^{\prime}>0$ such that

$|\Omega|^{\tfrac{D}{D-1}}\leq C^{\prime}|\partial\Omega|.$

for all finite $\Omega\subset VX$. The contrapositive of this says that if we
can find ‘large’ subsets with a ‘small’ boundary, then this subsequently
bounds the growth of $X$.

By [14, Thm. B], we know that $X$ is quasi-isometric to a complete Riemannian
plane $P$ which is not Gromov-hyperbolic. Let $K$ be some triangulation of
$P$, such that the inclusion $K\hookrightarrow P$ is a quasi-isometry.111Note
that the existence of such a triangulation is non-obvious and requires proof.
However, the plane $P$ constructed in the proof of Theorem B in [14] comes
naturally equipped with such a triangulation, so we shall ignore this
technicality.

In what follows, a local $(\lambda,\varepsilon)$-quasi-geodesic is a closed
loop $\ell$ such that any subpath of length less than
$\operatorname{length}(\ell)/2$ is a $(\lambda,\varepsilon)$-quasi-geodesic. A
theorem of Papasoglu states that a connected, locally finite graph is Gromov-
hyperbolic if and only if it has uniformly thin bigons [18]. As noted in [19,
Lem. 3.8], this easily implies that $K$ contains arbitrarily long simple loops
which are uniformly locally $(\lambda,\varepsilon)$-quasi-geodesic. Let
$\ell\subset K$ be such a loop. Then by the Jordan curve theorem, $\ell$
separates $P$ into a bounded and an unbounded component. Let $U\subset P$
denote the bounded component, so $U$ is an open disk. It is an immediate
consequence of the fact that $\ell$ is a local $(\lambda,\varepsilon)$-quasi-
geodesic that $U$ must contain a vertex $v$ of $K$ which lies at a distance of
at least $c\cdot\operatorname{length}(\ell)$ from $\ell$, where $c>0$ is some
constant depending only on $\lambda$, $\varepsilon$, and $K$. Let $\psi:K\to
X$ be a quasi-isometry. Let $y=\psi(v)$. It is easy to check that there must
exist some uniform $m>0$ such that $y$ is contained in a bounded connected
component of $X-N_{m}(\psi(\ell))$ (provided $\ell$ is taken to be
sufficiently long). In other words, if we write
$n=\operatorname{length}(\ell)$ and let $\Omega\subset VX$ denote the vertex
set of $U^{\prime}$, we see that

$|\partial\Omega|\leq c_{1}n,\ \ \ V(c_{2}n)\leq|\Omega|.$

where $c_{1},c_{2}>0$ are some uniform constants. By repeating this
construction with arbitrarily long $\ell$, we deduce by Varopoulos’ inequality
that $X$ has at most quadratic growth.

Building on deep results of Gromov [10] and Bass [1], a theorem of Woess [27,
Thm. 4.1] implies our desired conclusion.

######  Theorem 1.3 (Woess).

Let $X$ be a connected, locally finite, quasi-transitive graph with quadratic
growth. Then $G=\operatorname{Isom}(X)$ contains a compact normal subgroup $N$
such that $G/N$ is a finitely generated discrete group with a finite index
subgroup isomorphic to $\mathbf{Z}^{2}$ acting properly discontinuously on the
quotient graph $X/N$.

As noted in the hyperbolic case, the quotient map $X\to X/N$ is a quasi-
isometry, so by the Švarc–Milnor lemma we have that $X$ is quasi-isometric to
$\mathbf{Z}^{2}$. We thus deduce the following.

######  Theorem 1.4.

Let $X$ be a one-ended, connected, locally finite, quasi-transitive graph
which is not Gromov-hyperbolic. If $X$ is quasi-isometric to a complete
Riemannian plane, then $X$ is quasi-isometric to $\mathbf{Z}^{2}$, and thus to
a planar Cayley graph.

##  2\. Infinite-ended graphs

We will need the following theorem of Hamann [11], which is inspired by the
analogous result for finitely generated groups due to Papasoglu and Whyte
[20].

######  Theorem 2.1 (Hamann).

Let $X$, $Y$ be connected, locally finite, infinite-ended, quasi-transitive
graphs. Let $(T,\mathcal{V})$ and $(S,\mathcal{U})$ be canonical tree
decompositions of $X$, $Y$ with bounded adhesion and tight separators, such
that each torso has at most one end. Let $\Upsilon_{X}$, $\Upsilon_{Y}$ denote
the set of quasi-isometry classes of the one-ended torsos of the
aforementioned tree decompositions, without repetitions. Then $X$ and $Y$ are
quasi-isometric if and only if $\Upsilon_{X}=\Upsilon_{Y}$.

We direct the reader to [11] for definitions and statements relating to
(canonical) tree decompositions. Our main result now follows easily from
combining all that has come before.

*

###### Proof.

If $X$ is finite or two-ended then there is nothing to prove. If $X$ is one-
ended then we are done by Theorems 1.1 and 1.4.

Suppose then that $X$ is infinite-ended. By the main theorem of [14] we have
that $X$ is accessible. More precisely, by [14, Cor. C] we have that $X$
admits a canonical tree decomposition with bounded adhesion and tight
separators, where each infinite torso is quasi-isometric to a complete
Riemannian plane. If $\Upsilon$ denotes the set of quasi-isometry classes of
these one-ended torsos, then by Theorems 1.1 and 1.4 we have that $\Upsilon$
contains at most two elements, picked from the QI-classes of the hyperbolic
plane $\mathbf{H}^{2}$ and the Euclidean plane $\mathbf{R}^{2}$. Combining
this observation with Theorem 2.1, we deduce that $X$ is quasi-isometric to
one of

$F_{2},\ \ \mathbf{Z}^{2}\ast\mathbf{Z}^{2},\ \ \Sigma\ast\Sigma,\ \
\mathbf{Z}^{2}\ast\Sigma,$

where $\Sigma$ denotes some fixed cocompact Fuchsian group. The $F_{2}$-case
arises exactly when every torso is finite. In particular, we deduce that $X$
is quasi-isometric to a planar Cayley graph. ∎

######  Remark 2.2.

We conclude by remarking that every planar Cayley graph is quasi-isometric to
a ‘Kleinian function group’.

Recall that a Kleinian group $G$ is a discrete subgroup of
$\rm{PSL}(2,\mathbf{C})$. Such a group admits a canonical action on
$\mathbf{S}^{2}$ by Möbius transformations. The unique maximal open subspace
$\Omega$ of $\mathbf{S}^{2}$ on which $G$ acts properly discontinuously is
called the domain of discontinuity. If $\Omega$ contains a $G$-invariant
connected component then we call $G$ a function group. These groups can be
completely classified, see [15] or [16] for further reading.

There is a natural combinatorial characterisation of function groups in terms
of their Cayley graphs. Given a Cayley graph $\Gamma$ of a finitely generated
group $G$, a planar embedding $\vartheta:\Gamma\hookrightarrow\mathbf{S}^{2}$
is said to be covariant if the action of $G$ on $\Gamma$ sends facial paths to
facial paths. By a theorem of Levinson–Maskit [13], the class of finitely
generated groups group admitting a covariant embedding into $\mathbf{S}^{2}$
corresponds precisely to the Kleinian function groups (see also [7]).

Whitney’s planar embedding theorem [25, 21] states that 3-connected222Recall
that a graph $\Gamma$ is said to be 3-connected if it is connected and the
removal of any vertex or pair of vertices cannot disconnect $\Gamma$. planar
graphs embed uniquely into 2-sphere, up to post-composition with
homeomorphisms of the 2-sphere. Therefore, a 3-connected planar Cayley graph
is necessarily covariantly embedded. Thus, it is sufficient to find examples
of 3-connected planar Cayley graphs in each quasi-isometry class. The only
non-trivial case to consider is infinite-ended groups which are not virtually
free.

These can be constructed by amalgamating triangle groups along certain cyclic
subgroups. For example, consider the Euclidean (2,3,6)-triangle group

$H=\langle r,s,t\ ;\ (rs)^{2},\ (st)^{3},\ (rt)^{6}\rangle.$

The Cayley graph corresponding to the above presentation is planar. The order-
three cyclic subgroup $C=\langle st\rangle\leq H$ lies on the boundary of a
single face. This implies that the Cayley graph of the double $G=H\ast_{C}H$
given by the presentation

$G=\langle a,b,c,r,s,t\ ;\ (rs)^{2},\ (st)^{3},\ (rt)^{6},\ (ab)^{2},\
(bc)^{3},\ (ac)^{6},\ bc=st\rangle$

is planar. Since $|C|\geq 3$ we have that this Cayley graph is also
3-connected. Thus $G$ is a Kleinian function group by the above discussion.

The theorem of Papasoglu and Whyte [20] implies that $G$ is indeed quasi-
isometric to $\mathbf{Z}^{2}\ast\mathbf{Z}^{2}$. A similar construction is
possible within the other infinite-ended quasi-isometry classes, which we
leave as an exercise to the interested reader.

### Acknowledgements

I am very grateful to Yves Cornulier for taking the time to share some insight
into the structure of isometry groups of hyperbolic graphs with circular
boundaries, and also to Louis Esperet, Agelos Georgakopoulos, Matthias Hamann,
and Panos Papasoglu for comments on an earlier version of this note. I also
thank Lawk Mineh for providing some very sobering counterexamples.

## References

  * [1] H. Bass.  The degree of polynomial growth of finitely generated nilpotent groups.  Proceedings of the London Mathematical Society, 3(4):603–614, 1972.
  * [2] A. Casson and D. Jungreis.  Convergence groups and Seifert fibered 3-manifolds.  Inventiones mathematicae, 118(1):441–456, 1994.
  * [3] Y. Cornulier and P. de La Harpe.  Metric geometry of locally compact groups.  arXiv preprint arXiv:1403.3796, 2014.
  * [4] A. Eskin, D. Fisher, and K. Whyte.  Coarse differentiation of quasi-isometries I: Spaces not quasi-isometric to cayley graphs.  Annals of Mathematics, pages 221–260, 2012.
  * [5] L. Esperet and U. Giocanti.  Coarse geometry of quasi-transitive graphs beyond planarity.  arXiv preprint arXiv:2312.08902, 2023.
  * [6] D. Gabai.  Convergence groups are Fuchsian groups.  Annals of Mathematics, 136(3):447–510, 1992.
  * [7] A. Georgakopoulos.  On planar Cayley graphs and Kleinian groups.  Transactions of the American Mathematical Society, 373(7):4649–4684, 2020.
  * [8] A. Georgakopoulos and P. Papasoglu.  Graph minors and metric spaces.  arXiv preprint arXiv:2305.07456, 2023.
  * [9] É. Ghys.  Groups acting on the circle.  Enseignement Mathematique, 47(3/4):329–408, 2001.
  * [10] M. Gromov.  Groups of polynomial growth and expanding maps (with an appendix by Jacques Tits).  Publications Mathématiques de l’IHÉS, 53:53–78, 1981.
  * [11] M. Hamann.  Tree amalgamations and quasi-isometries.  arXiv preprint arXiv:1812.04987, 2018.
  * [12] M. Hamann.  Quasi-transitive $K_{\infty}$-minor free graphs.  arXiv preprint arXiv:2405.17218, 2024.
  * [13] H. Levinson and B. Maskit.  Special embeddings of Cayley diagrams.  Journal of Combinatorial Theory, Series B, 18(1):12–17, 1975.
  * [14] J. MacManus.  Accessibility, planar graphs, and quasi-isometries.  arXiv preprint arXiv:2310.15242, 2023.
  * [15] A. Marden.  The geometry of finitely generated Kleinian groups.  Annals of Mathematics, 99(3):383–462, 1974.
  * [16] B. Maskit.  Kleinian groups, volume 287.  Springer Science & Business Media, 2012.
  * [17] G. Mess.  The Seifert conjecture and groups which are coarse quasiisometric to planes.  preprint, 1988.
  * [18] P. Papasoglu.  Strongly geodesically automatic groups are hyperbolic.  Inventiones mathematicae, 121:323–334, 1995.
  * [19] P. Papasoglu.  Quasi-isometry invariance of group splittings.  Annals of mathematics, pages 759–830, 2005.
  * [20] P. Papasoglu and K. Whyte.  Quasi-isometries between groups with infinitely many ends.  Commentarii Mathematici Helvetici, 77(1):133–144, 2002.
  * [21] R. Richter and C. Thomassen.  3-connected planar spaces uniquely embed in the sphere.  Transactions of the American Mathematical Society, 354(11):4585–4595, 2002.
  * [22] L. Saloff-Coste.  Isoperimetric inequalities and decay of iterated kernels for almost-transitive Markov chains.  Combinatorics, Probability and Computing, 4(4):419–442, 1995.
  * [23] C. Thomassen and W. Woess.  Vertex-transitive graphs and accessibility.  Journal of Combinatorial Theory, Series B, 58(2):248–268, 1993.
  * [24] P. Tukia.  Homeomorphic conjugates of Fuchsian groups.  Journal für die reine und angewandte Mathematik, 390:1–54, 1988.
  * [25] H. Whitney.  Congruent graphs and the connectivity of graphs.  Hassler Whitney Collected Papers, pages 61–79, 1992.
  * [26] W. Woess.  Topological groups and infinite graphs.  Discrete mathematics, 95(1-3):373–384, 1991.
  * [27] W. Woess.  Topological groups and recurrence of quasi transitive graphs.  Rendiconti del Seminario Matematico e Fisico di Milano, 64(1):185–213, 1994.


* [517] A. Nichol and P. Dhariwal, “Improved Denoising Diffusion Probabilistic Models,” Feb. 2021.  arXiv:2102.09672 [cs, stat].
  * [518] D. Rezende and S. Mohamed, “Variational Inference with Normalizing Flows,” in Proceedings of the 32nd International Conference on Machine Learning, pp. 1530–1538, PMLR, June 2015.
  * [519] I. Kobyzev, S. J. D. Prince, and M. A. Brubaker, “Normalizing Flows: An Introduction and Review of Current Methods,” IEEE Transactions on Pattern Analysis and Machine Intelligence, vol. 43, pp. 3964–3979, Nov. 2021\.  arXiv:1908.09257 [cs, stat].
  * [520] L. Mosser, O. Dubrule, and M. J. Blunt, “Reconstruction of three-dimensional porous media using generative adversarial neural networks,” Physical Review E, vol. 96, p. 043309, Oct. 2017.
  * [521] J. Feng, X. He, Q. Teng, C. Ren, H. Chen, and Y. Li, “Reconstruction of porous media from extremely limited information using conditional generative adversarial networks,” Physical Review E, vol. 100, p. 033308, Sept. 2019\.
  * [522] R. Shams, M. Masihi, R. B. Boozarjomehry, and M. J. Blunt, “Coupled generative adversarial and auto-encoder neural networks to reconstruct three-dimensional multi-scale porous media,” Journal of Petroleum Science and Engineering, vol. 186, p. 106794, Mar. 2020.
  * [523] P. Xia, H. Bai, and T. Zhang, “Multi-scale reconstruction of porous media based on progressively growing generative adversarial networks,” Stochastic Environmental Research and Risk Assessment, vol. 36, pp. 3685–3705, Nov. 2022.
  * [524] A. Henkes and H. Wessels, “Three-dimensional microstructure generation using generative adversarial neural networks in the context of continuum micromechanics,” Computer Methods in Applied Mechanics and Engineering, vol. 400, p. 115497, Oct. 2022.
  * [525] S. Rawat and M. H. H. Shen, “A novel topology design approach using an integrated deep learning network architecture,” Jan. 2019.
  * [526] K. Yaji, S. Yamasaki, and K. Fujita, “Data-driven multifidelity topology design using a deep generative model: Application to forced convection heat transfer problems,” Computer Methods in Applied Mechanics and Engineering, vol. 388, p. 114284, Jan. 2022.
  * [527] K.-H. Lee and G. J. Yun, “Microstructure reconstruction using diffusion-based generative models,” Jan. 2023.  arXiv:2211.10949 [cond-mat, physics:physics].
  * [528] C. Düreth, P. Seibert, D. Rücker, S. Handford, M. Kästner, and M. Gude, “Conditional diffusion-based microstructure reconstruction,” Materials Today Communications, vol. 35, p. 105608, June 2023.
  * [529] N. N. Vlassis and W. Sun, “Denoising diffusion algorithm for inverse design of microstructures with fine-tuned nonlinear material properties,” Computer Methods in Applied Mechanics and Engineering, vol. 413, p. 116126, Aug. 2023.
  * [530] J. Feng, Q. Teng, B. Li, X. He, H. Chen, and Y. Li, “An end-to-end three-dimensional reconstruction framework of porous media from a single two-dimensional image based on deep learning,” Computer Methods in Applied Mechanics and Engineering, vol. 368, p. 113043, Aug. 2020.
  * [531] S. Kench and S. J. Cooper, “Generating three-dimensional structures from a two-dimensional slice with generative adversarial network-based dimensionality expansion,” Nature Machine Intelligence, vol. 3, pp. 299–305, Apr. 2021.
  * [532] Y. Li, P. Jian, and G. Han, “Cascaded Progressive Generative Adversarial Networks for Reconstructing Three-Dimensional Grayscale Core Images From a Single Two-Dimensional Image,” Frontiers in Physics, vol. 10, 2022.
  * [533] F. Zhang, X. He, Q. Teng, X. Wu, and X. Dong, “3D-PMRNN: Reconstructing three-dimensional porous media from the two-dimensional image with recurrent neural network,” Journal of Petroleum Science and Engineering, vol. 208, p. 109652, Jan. 2022.
  * [534] Q. Zheng and D. Zhang, “RockGPT: reconstructing three-dimensional digital rocks from single two-dimensional slice with deep learning,” Computational Geosciences, vol. 26, pp. 677–696, June 2022.
  * [535] J. Phan, L. Ruspini, G. Kiss, and F. Lindseth, “Size-invariant 3D generation from a single 2D rock image,” Journal of Petroleum Science and Engineering, vol. 215, p. 110648, Aug. 2022.
  * [536] F. Zhang, Q. Teng, H. Chen, X. He, and X. Dong, “Slice-to-voxel stochastic reconstructions on porous media with hybrid deep generative model,” Computational Materials Science, vol. 186, p. 110018, Jan. 2021.
  * [537] S. Rawat and M. H. Shen, “Application of Adversarial Networks for 3D Structural Topology Optimization,” pp. 2019–01–0829, Apr. 2019.
  * [538] S. Rawat and M.-H. H. Shen, “A Novel Topology Optimization Approach using Conditional Deep Learning,” Jan. 2019.
  * [539] M.-H. H. Shen and L. Chen, “A New CGAN Technique for Constrained Topology Design Optimization,” Apr. 2019.
  * [540] H. Wessels, C. Böhm, F. Aldakheel, M. Hüpgen, M. Haist, L. Lohaus, and P. Wriggers, “Computational Homogenization Using Convolutional Neural Networks,” in Current Trends and Open Problems in Computational Mechanics (F. Aldakheel, B. Hudobivnik, M. Soleimani, H. Wessels, C. Weißenfels, and M. Marino, eds.), pp. 569–579, Cham: Springer International Publishing, 2022.
  * [541] L. Mosser, O. Dubrule, and M. J. Blunt, “Stochastic Seismic Waveform Inversion Using Generative Adversarial Networks as a Geological Prior,” Mathematical Geosciences, vol. 52, pp. 53–79, Jan. 2020.
  * [542] T. Guo, D. J. Lohan, R. Cang, M. Y. Ren, and J. T. Allison, “An Indirect Design Representation for Topology Optimization Using Variational Autoencoder and Style Transfer,” in 2018 AIAA/ASCE/AHS/ASC Structures, Structural Dynamics, and Materials Conference, AIAA SciTech Forum, American Institute of Aeronautics and Astronautics, Jan. 2018.
  * [543] P. S. Vulimiri, H. Deng, F. Dugast, X. Zhang, and A. C. To, “Integrating Geometric Data into Topology Optimization via Neural Style Transfer,” Materials, vol. 14, p. 4551, Jan. 2021.
  * [544] L. Gatys, A. Ecker, and M. Bethge, “A Neural Algorithm of Artistic Style,” Journal of Vision, vol. 16, p. 326, Sept. 2016.
  * [545] S. Khan, K. Goucher-Lambert, K. Kostas, and P. Kaklis, “ShipHullGAN: A generic parametric modeller for ship hull design using deep convolutional generative model,” Computer Methods in Applied Mechanics and Engineering, vol. 411, p. 116051, June 2023.
  * [546] Q. Chen, J. Wang, P. Pope, W. (Wayne) Chen, and M. Fuge, “Inverse Design of Two-Dimensional Airfoils Using Conditional Generative Models and Surrogate Log-Likelihoods,” Journal of Mechanical Design, vol. 144, p. 021712, Feb. 2022.
  * [547] W. Chen and M. Fuge, “B\’ezierGAN: Automatic Generation of Smooth Curves from Interpretable Low-Dimensional Parameters,” Jan. 2021.  arXiv:1808.08871 [cs, stat].
  * [548] W. Chen and F. Ahmed, “MO-PaDGAN: Reparameterizing Engineering Designs for augmented multi-objective optimization,” Applied Soft Computing, vol. 113, p. 107909, Dec. 2021.
  * [549] A. Richardson, “Generative Adversarial Networks for Model Order Reduction in Seismic Full-Waveform Inversion,” June 2018.  arXiv:1806.00828 [physics].
  * [550] Y. Zhang, P. Seibert, A. Otto, A. Raßloff, M. Ambati, and M. Kästner, “DA-VEGAN: Differentiably Augmenting VAE-GAN for microstructure reconstruction from extremely small data sets,” Feb. 2023.  arXiv:2303.03403 [cs].
  * [551] W. Chen and F. Ahmed, “PaDGAN: Learning to Generate High-Quality Novel Designs,” Journal of Mechanical Design, vol. 143, p. 031703, Mar. 2021.
  * [552] A. Kulesza and B. Taskar, “Determinantal point processes for machine learning,” Foundations and Trends® in Machine Learning, vol. 5, no. 2-3, pp. 123–286, 2012.  arXiv:1207.6083 [cs, stat].
  * [553] S. J. Bates, J. Sienz, and D. S. Langley, “Formulation of the Audze–Eglais Uniform Latin Hypercube design of experiments,” Advances in Engineering Software, vol. 34, pp. 493–506, Aug. 2003.
  * [554] A. Heyrani Nobari, M. F. Rashad, and F. Ahmed, “CreativeGAN: Editing Generative Adversarial Networks for Creative Design Synthesis,” in Volume 3A: 47th Design Automation Conference (DAC), (Virtual, Online), p. V03AT03A002, American Society of Mechanical Engineers, Aug. 2021.
  * [555] D. Bau, S. Liu, T. Wang, J.-Y. Zhu, and A. Torralba, “Rewriting a Deep Generative Model,” July 2020.  arXiv:2007.15646 [cs].
  * [556] A. Elgammal, B. Liu, M. Elhoseiny, and M. Mazzone, “CAN: Creative Adversarial Networks, Generating ”Art” by Learning About Styles and Deviating from Style Norms,” June 2017.  arXiv:1706.07068 [cs].
  * [557] S. Oh, Y. Jung, S. Kim, I. Lee, and N. Kang, “Deep Generative Design: Integration of Topology Optimization and Generative Models,” Journal of Mechanical Design, vol. 141, p. 111405, Nov. 2019.
  * [558] M. Greminger, “Generative Adversarial Networks With Synthetic Training Data for Enforcing Manufacturing Constraints on Topology Optimization,” in Volume 11A: 46th Design Automation Conference (DAC), (Virtual, Online), p. V11AT11A005, American Society of Mechanical Engineers, Aug. 2020.
  * [559] S. Yoo, S. Lee, S. Kim, K. H. Hwang, J. H. Park, and N. Kang, “Integrating deep learning into CAD/CAE system: generative design and evaluation of 3D conceptual wheel,” Structural and Multidisciplinary Optimization, vol. 64, pp. 2725–2747, Oct. 2021.
  * [560] W. Zhang, Y. Wang, Z. Du, C. Liu, S.-K. Youn, and X. Guo, “Machine-learning assisted topology optimization for architectural design with artistic flavor,” Computer Methods in Applied Mechanics and Engineering, vol. 413, p. 116041, Aug. 2023.
  * [561] M. P. Bendsøe and O. Sigmund, Topology optimization: theory, methods, and applications.  Berlin ; New York: Springer, 2003.
  * [562] F. Yang and J. Ma, “FWIGAN: Full‐Waveform Inversion via a Physics‐Informed Generative Adversarial Network,” Journal of Geophysical Research: Solid Earth, vol. 128, p. e2022JB025493, Apr. 2023.
  * [563] S. Radhakrishnan, V. Bharadwaj, V. Manjunath, and R. Srinath, “Creative Intelligence – Automating Car Design Studio with Generative Adversarial Networks (GAN),” in Machine Learning and Knowledge Extraction (A. Holzinger, P. Kieseberg, A. M. Tjoa, and E. Weippl, eds.), Lecture Notes in Computer Science, (Cham), pp. 160–175, Springer International Publishing, 2018.
  * [564] W. Chen and M. Fuge, “Synthesizing Designs With Interpart Dependencies Using Hierarchical Generative Adversarial Networks,” Journal of Mechanical Design, vol. 141, p. 111403, Nov. 2019.
  * [565] Z. Nie, T. Lin, H. Jiang, and L. B. Kara, “TopologyGAN: Topology Optimization Using Generative Adversarial Networks Based on Physical Fields Over the Initial Domain,” Mar. 2020.
  * [566] N. Hertlein, P. R. Buskohl, A. Gillman, K. Vemaganti, and S. Anand, “Generative adversarial network for early-stage design flexibility in topology optimization for additive manufacturing,” Journal of Manufacturing Systems, vol. 59, pp. 675–685, Apr. 2021.
  * [567] A. Heyrani Nobari, W. W. Chen, and F. Ahmed, “RANGE-GAN: Design Synthesis Under Constraints Using Conditional Generative Adversarial Networks,” Journal of Mechanical Design, pp. 1–16, Sept. 2021.
  * [568] J. Wang, W. W. Chen, D. Da, M. Fuge, and R. Rai, “IH-GAN: A conditional generative model for implicit surface-based inverse design of cellular structures,” Computer Methods in Applied Mechanics and Engineering, vol. 396, p. 115060, June 2022.
  * [569] L. Duque, G. Gutiérrez, C. Arias, A. Rüger, and H. Jaramillo, “Automated Velocity Estimation by Deep Learning Based Seismic-to-Velocity Mapping,” vol. 2019, pp. 1–5, European Association of Geoscientists & Engineers, June 2019.
  * [570] Y.-Q. Wang, Q. Wang, W.-K. Lu, Q. Ge, and X.-F. Yan, “Seismic impedance inversion based on cycle-consistent generative adversarial network,” Petroleum Science, vol. 19, pp. 147–161, Feb. 2022.
  * [571] J.-Y. Zhu, T. Park, P. Isola, and A. A. Efros, “Unpaired Image-to-Image Translation using Cycle-Consistent Adversarial Networks,” Aug. 2020\.  arXiv:1703.10593 [cs].
  * [572] B. Li, C. Huang, X. Li, S. Zheng, and J. Hong, “Non-iterative structural topology optimization using deep learning,” Computer-Aided Design, vol. 115, pp. 172–180, Oct. 2019.
  * [573] Y. Xie, E. Franz, M. Chu, and N. Thuerey, “tempoGAN: a temporally coherent, volumetric GAN for super-resolution fluid flow,” ACM Transactions on Graphics, vol. 37, pp. 95:1–95:15, July 2018.
  * [574] G. Pang, C. Shen, L. Cao, and A. V. D. Hengel, “Deep Learning for Anomaly Detection: A Review,” ACM Computing Surveys, vol. 54, pp. 1–38, Mar. 2022.
  * [575] S. Hawkins, H. He, G. Williams, and R. Baxter, “Outlier Detection Using Replicator Neural Networks,” in Data Warehousing and Knowledge Discovery (Y. Kambayashi, W. Winiwarter, and M. Arikawa, eds.), Lecture Notes in Computer Science, (Berlin, Heidelberg), pp. 170–180, Springer, 2002.
  * [576] T. Schlegl, P. Seeböck, S. M. Waldstein, U. Schmidt-Erfurth, and G. Langs, “Unsupervised Anomaly Detection with Generative Adversarial Networks to Guide Marker Discovery,” in Information Processing in Medical Imaging (M. Niethammer, M. Styner, S. Aylward, H. Zhu, I. Oguz, P.-T. Yap, and D. Shen, eds.), Lecture Notes in Computer Science, (Cham), pp. 146–157, Springer International Publishing, 2017.
  * [577] H. Zenati, C. S. Foo, B. Lecouat, G. Manek, and V. R. Chandrasekhar, “Efficient GAN-Based Anomaly Detection,” May 2019.  arXiv:1802.06222 [cs, stat].
  * [578] T. Schlegl, P. Seeböck, S. M. Waldstein, G. Langs, and U. Schmidt-Erfurth, “f-AnoGAN: Fast unsupervised anomaly detection with generative adversarial networks,” Medical Image Analysis, vol. 54, pp. 30–44, May 2019.
  * [579] A. Henkes, L. Herrmann, H. Wessels, and S. Kollmannsberger, “GAN enables Microstructure Monitoring for Additive Manufacturing of Complex Structures,” in preparation, 2023.
  * [580] V. Mnih, K. Kavukcuoglu, D. Silver, A. A. Rusu, J. Veness, M. G. Bellemare, A. Graves, M. Riedmiller, A. K. Fidjeland, G. Ostrovski, S. Petersen, C. Beattie, A. Sadik, I. Antonoglou, H. King, D. Kumaran, D. Wierstra, S. Legg, and D. Hassabis, “Human-level control through deep reinforcement learning,” Nature, vol. 518, pp. 529–533, Feb. 2015.
  * [581] D. Silver, J. Schrittwieser, K. Simonyan, I. Antonoglou, A. Huang, A. Guez, T. Hubert, L. Baker, M. Lai, A. Bolton, Y. Chen, T. Lillicrap, F. Hui, L. Sifre, G. van den Driessche, T. Graepel, and D. Hassabis, “Mastering the game of Go without human knowledge,” Nature, vol. 550, pp. 354–359, Oct. 2017.
  * [582] O. Vinyals, I. Babuschkin, W. M. Czarnecki, M. Mathieu, A. Dudzik, J. Chung, D. H. Choi, R. Powell, T. Ewalds, P. Georgiev, J. Oh, D. Horgan, M. Kroiss, I. Danihelka, A. Huang, L. Sifre, T. Cai, J. P. Agapiou, M. Jaderberg, A. S. Vezhnevets, R. Leblond, T. Pohlen, V. Dalibard, D. Budden, Y. Sulsky, J. Molloy, T. L. Paine, C. Gulcehre, Z. Wang, T. Pfaff, Y. Wu, R. Ring, D. Yogatama, D. Wünsch, K. McKinney, O. Smith, T. Schaul, T. Lillicrap, K. Kavukcuoglu, D. Hassabis, C. Apps, and D. Silver, “Grandmaster level in StarCraft II using multi-agent reinforcement learning,” Nature, vol. 575, pp. 350–354, Nov. 2019.
  * [583] J. Kober, J. A. Bagnell, and J. Peters, “Reinforcement learning in robotics: A survey,” The International Journal of Robotics Research, vol. 32, pp. 1238–1274, Sept. 2013.
  * [584] H. Kim, M. Jordan, S. Sastry, and A. Ng, “Autonomous Helicopter Flight via Reinforcement Learning,” in Advances in Neural Information Processing Systems, vol. 16, MIT Press, 2003.
  * [585] P. Abbeel, A. Coates, M. Quigley, and A. Ng, “An Application of Reinforcement Learning to Aerobatic Helicopter Flight,” in Advances in Neural Information Processing Systems, vol. 19, MIT Press, 2006.
  * [586] P. Abbeel, A. Coates, and A. Y. Ng, “Autonomous Helicopter Aerobatics through Apprenticeship Learning,” The International Journal of Robotics Research, vol. 29, pp. 1608–1639, Nov. 2010.
  * [587] R. S. Sutton and A. G. Barto, Reinforcement learning: an introduction.  Adaptive computation and machine learning series, Cambridge, Massachusetts: The MIT Press, second edition ed., 2018.
  * [588] G. Novati, S. Verma, D. Alexeev, D. Rossinelli, W. M. van Rees, and P. Koumoutsakos, “Synchronised Swimming of Two Fish,” Bioinspiration & Biomimetics, vol. 12, p. 036001, Mar. 2017.  arXiv:1610.04248 [physics].
  * [589] S. Verma, G. Novati, and P. Koumoutsakos, “Efficient collective swimming by harnessing vortices through deep reinforcement learning,” Proceedings of the National Academy of Sciences, vol. 115, pp. 5849–5854, June 2018.
  * [590] P. Ma, Y. Tian, Z. Pan, B. Ren, and D. Manocha, “Fluid directed rigid body control using deep reinforcement learning,” ACM Transactions on Graphics, vol. 37, pp. 96:1–96:11, July 2018.
  * [591] J. Rabault, M. Kuchta, A. Jensen, U. Réglade, and N. Cerardi, “Artificial neural networks trained through deep reinforcement learning discover control strategies for active flow control,” Journal of Fluid Mechanics, vol. 865, pp. 281–302, Apr. 2019.
  * [592] D. Fan, L. Yang, Z. Wang, M. S. Triantafyllou, and G. E. Karniadakis, “Reinforcement learning for bluff body active flow control in experiments and simulations,” Proceedings of the National Academy of Sciences, vol. 117, pp. 26091–26098, Oct. 2020.
  * [593] J. Xu, T. Du, M. Foshey, B. Li, B. Zhu, A. Schulz, and W. Matusik, “Learning to fly: computational controller design for hybrid UAVs with reinforcement learning,” ACM Transactions on Graphics, vol. 38, pp. 42:1–42:12, July 2019.
  * [594] X. Y. Lee, A. Balu, D. Stoecklein, B. Ganapathysubramanian, and S. Sarkar, “Flow Shape Design for Microfluidic Devices Using Deep Reinforcement Learning,” Nov. 2018.  arXiv:1811.12444 [cs, stat].
  * [595] K. Wang and W. Sun, “Meta-modeling game for deriving theory-consistent, microstructure-based traction–separation laws via deep reinforcement learning,” Computer Methods in Applied Mechanics and Engineering, vol. 346, pp. 216–241, Apr. 2019.
  * [596] M. P. Bendsøe, “Optimal shape design as a material distribution problem,” Structural optimization, vol. 1, pp. 193–202, Dec. 1989.
  * [597] M. P. Bendsøe and O. Sigmund, Topology Optimization.  Berlin, Heidelberg: Springer, 2004.
  * [598] K. Hayashi and M. Ohsaki, “Reinforcement Learning and Graph Embedding for Binary Truss Topology Optimization Under Stress and Displacement Constraints,” Frontiers in Built Environment, vol. 6, 2020\.
  * [599] S. Zhu, M. Ohsaki, K. Hayashi, and X. Guo, “Machine-specified ground structures for topology optimization of binary trusses using graph embedding policy network,” Advances in Engineering Software, vol. 159, p. 103032, Sept. 2021.
  * [600] H. Sun and L. Ma, “Generative Design by Using Exploration Approaches of Reinforcement Learning in Density-Based Structural Topology Optimization,” Designs, vol. 4, p. 10, June 2020.
  * [601] S. Jang, S. Yoo, and N. Kang, “Generative Design by Reinforcement Learning: Enhancing the Diversity of Topology Optimization Designs,” Computer-Aided Design, vol. 146, p. 103225, May 2022.
  * [602] J. Han, A. Jentzen, and W. E, “Solving high-dimensional partial differential equations using deep learning,” Proceedings of the National Academy of Sciences, vol. 115, pp. 8505–8510, Aug. 2018.
  * [603] W. E and B. Yu, “The Deep Ritz method: A deep learning-based numerical algorithm for solving variational problems,” arXiv:1710.00211 [cs, stat], Sept. 2017.  arXiv: 1710.00211.
  * [604] J. Yang, T. Dzanic, B. Petersen, J. Kudo, K. Mittal, V. Tomov, J.-S. Camier, T. Zhao, H. Zha, T. Kolev, R. Anderson, and D. Faissol, “Reinforcement Learning for Adaptive Mesh Refinement,” in Proceedings of The 26th International Conference on Artificial Intelligence and Statistics, pp. 5997–6014, PMLR, Apr. 2023.
  * [605] J. Rabault and A. Kuhnle, “Accelerating Deep Reinforcement Learning strategies of Flow Control through a multi-environment approach,” Physics of Fluids, vol. 31, p. 094105, Sept. 2019.  arXiv:1906.10382 [physics].
  * [606] G. Novati, H. L. de Laroussilhe, and P. Koumoutsakos, “Automating Turbulence Modeling by Multi-Agent Reinforcement Learning,” arXiv:2005.09023 [physics], Oct. 2020.  arXiv: 2005.09023.
  * [607] X.-Y. Liu and J.-X. Wang, “Physics-informed Dyna-style model-based deep reinforcement learning for dynamic control,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 477, p. 20210618, Nov. 2021.
  * [608] H. Shi, Y. Zhou, K. Wu, S. Chen, B. Ran, and Q. Nie, “Physics-informed deep reinforcement learning-based integrated two-dimensional car-following control strategy for connected automated vehicles,” Knowledge-Based Systems, vol. 269, p. 110485, June 2023.
  * [609] A. Ramesh and B. Ravindran, “Physics-Informed Model-Based Reinforcement Learning,” May 2023.  arXiv:2212.02179 [cs].
  * [610] C. Rodwell and P. Tallapragada, “Physics-informed reinforcement learning for motion control of a fish-like swimming robot,” Scientific Reports, vol. 13, p. 10754, July 2023.
  * [611] R. S. Sutton, “Dyna, an integrated architecture for learning, planning, and reacting,” ACM SIGART Bulletin, vol. 2, pp. 160–163, July 1991.
  * [612] M. Janner, J. Fu, M. Zhang, and S. Levine, “When to trust your model: model-based policy optimization,” in Proceedings of the 33rd International Conference on Neural Information Processing Systems, no. 1122, pp. 12519–12530, Red Hook, NY, USA: Curran Associates Inc., Dec. 2019.
  * [613] L. Kaiser, M. Babaeizadeh, P. Milos, B. Osinski, R. H. Campbell, K. Czechowski, D. Erhan, C. Finn, P. Kozakowski, S. Levine, A. Mohiuddin, R. Sepassi, G. Tucker, and H. Michalewski, “Model-Based Reinforcement Learning for Atari,” Feb. 2020.  arXiv:1903.00374 [cs, stat].
  * [614] Y. Luo, H. Xu, Y. Li, Y. Tian, T. Darrell, and T. Ma, “Algorithmic Framework for Model-based Deep Reinforcement Learning with Theoretical Guarantees,” Feb. 2021.  arXiv:1807.03858 [cs, stat].
  * [615] M. P. Deisenroth and C. E. Rasmussen, “PILCO: a model-based and data-efficient approach to policy search,” in Proceedings of the 28th International Conference on International Conference on Machine Learning, ICML’11, (Madison, WI, USA), pp. 465–472, Omnipress, June 2011\.
  * [616] S. Levine and P. Abbeel, “Learning Neural Network Policies with Guided Policy Search under Unknown Dynamics,” in Advances in Neural Information Processing Systems, vol. 27, Curran Associates, Inc., 2014\.
  * [617] N. Heess, G. Wayne, D. Silver, T. Lillicrap, T. Erez, and Y. Tassa, “Learning Continuous Control Policies by Stochastic Value Gradients,” in Advances in Neural Information Processing Systems, vol. 28, Curran Associates, Inc., 2015.
  * [618] I. Clavera, V. Fu, and P. Abbeel, “Model-Augmented Actor-Critic: Backpropagating through Paths,” May 2020.  arXiv:2005.08068 [cs, stat].
  * [619] D. Hafner, T. Lillicrap, J. Ba, and M. Norouzi, “Dream to Control: Learning Behaviors by Latent Imagination,” Mar. 2020.  arXiv:1912.01603 [cs].
  * [620] D. Hafner, T. Lillicrap, M. Norouzi, and J. Ba, “Mastering Atari with Discrete World Models,” Feb. 2022.  arXiv:2010.02193 [cs, stat].
  * [621] R. J. Williams, “Simple statistical gradient-following algorithms for connectionist reinforcement learning,” Machine Learning, vol. 8, pp. 229–256, May 1992.
  * [622] R. S. Sutton, D. McAllester, S. Singh, and Y. Mansour, “Policy Gradient Methods for Reinforcement Learning with Function Approximation,” in Advances in Neural Information Processing Systems, vol. 12, MIT Press, 1999.
  * [623] S. M. Kakade, “A Natural Policy Gradient,” in Advances in Neural Information Processing Systems, vol. 14, MIT Press, 2001.
  * [624] D. Silver, G. Lever, N. Heess, T. Degris, D. Wierstra, and M. Riedmiller, “Deterministic Policy Gradient Algorithms,” in Proceedings of the 31st International Conference on Machine Learning, pp. 387–395, PMLR, Jan. 2014.
  * [625] J. Schulman, S. Levine, P. Abbeel, M. Jordan, and P. Moritz, “Trust Region Policy Optimization,” in Proceedings of the 32nd International Conference on Machine Learning, pp. 1889–1897, PMLR, June 2015.
  * [626] C. J. C. H. Watkins and P. Dayan, “Q-learning,” Machine Learning, vol. 8, pp. 279–292, May 1992.
  * [627] H. v. Hasselt, A. Guez, and D. Silver, “Deep reinforcement learning with double Q-Learning,” in Proceedings of the Thirtieth AAAI Conference on Artificial Intelligence, AAAI’16, (Phoenix, Arizona), pp. 2094–2100, AAAI Press, Feb. 2016.
  * [628] Z. Wang, T. Schaul, M. Hessel, H. Van Hasselt, M. Lanctot, and N. De Freitas, “Dueling network architectures for deep reinforcement learning,” in Proceedings of the 33rd International Conference on International Conference on Machine Learning - Volume 48, ICML’16, (New York, NY, USA), pp. 1995–2003, JMLR.org, June 2016.
  * [629] J. Schulman, F. Wolski, P. Dhariwal, A. Radford, and O. Klimov, “Proximal Policy Optimization Algorithms,” Aug. 2017.  arXiv:1707.06347 [cs].
  * [630] R. Bellman, “A Markovian Decision Process,” Journal of Mathematics and Mechanics, vol. 6, no. 5, pp. 679–684, 1957.
  * [631] I. C. Dolcetta and H. Ishii, “Approximate solutions of the bellman equation of deterministic control theory,” Applied Mathematics & Optimization, vol. 11, pp. 161–181, Feb. 1984.
  * [632] R. S. Sutton, “Learning to predict by the methods of temporal differences,” Machine Learning, vol. 3, pp. 9–44, Aug. 1988.
  * [633] S. J. Bradtke and A. G. Barto, “Linear Least-Squares algorithms for temporal difference learning,” Machine Learning, vol. 22, pp. 33–57, Mar. 1996.