R ES E A RC H

CANCER GENOMICS Our scTrio-seq (single-cell triple omics sequenc-

ing) technique (13) can assess somatic copy num-
ber alterations (SCNAs), DNA methylation, and
Single-cell multiomics sequencing transcriptome information simultaneously from
an individual cell. Here, we introduce scTrio-seq2,

and analyses of human which integrates single-cell whole-genome bisulfite

sequencing (scBS-seq) (14) and improves detec-
tion efficiencies (fig. S1 and table S1). In this
colorectal cancer study, we performed multiregional sampling and
generated scTrio-seq2 profiles for 12 CRC patients
Shuhui Bian1,2,3*, Yu Hou1,2*, Xin Zhou4*, Xianlong Li1,2*, Jun Yong1,5*,
(stage III or IV) (Fig. 1, fig. S2, and table S1). In
total, ~1900 single cells passed quality control.
Yicheng Wang1,2*, Wendong Wang4, Jia Yan1,2, Boqiang Hu1,2, Hongshan Guo1,2,
Paired primary tumors and lymphatic or dis-
Jilian Wang4, Shuai Gao1,2, Yunuo Mao1,2, Ji Dong1,2, Ping Zhu1,2,3, Dianrong Xiu4,
tant metastases were obtained from 10 patients
Liying Yan1,5, Lu Wen1,2, Jie Qiao1,3,5,6†, Fuchou Tang1,2,3†, Wei Fu4†
(table S1). For patient CRC01, we obtained 534
single cells (after quality control) from adjacent
Although genomic instability, epigenetic abnormality, and gene expression dysregulation are
normal colon (NC) tissue and 16 tumor regions,
hallmarks of colorectal cancer, these features have not been simultaneously analyzed at
including the primary tumor (PT), the lymph
single-cell resolution. Using optimized single-cell multiomics sequencing together with
node metastasis (LN), the liver metastasis (ML),
multiregional sampling of the primary tumor and lymphatic and distant metastases, we
and a posttreatment liver metastasis (MP) after
developed insights beyond intratumoral heterogeneity. Genome-wide DNA methylation levels
chemotherapy (fig. S2).

Downloaded from http://science.sciencemag.org/ on January 24, 2019

were relatively consistent within a single genetic sublineage. The genome-wide DNA
Most cancer cells from six of our study pa-
demethylation patterns of cancer cells were consistent in all 10 patients whose DNA we
tients (CRC01, CRC03, CRC04, CRC06, CRC09,
sequenced. The cancer cells’ DNA demethylation degrees clearly correlated with the densities
of the heterochromatin-associated histone modification H3K9me3 of normal tissue and 1
Beijing Advanced Innovation Center for Genomics, College
those of repetitive element long interspersed nuclear element 1. Our work demonstrates the
of Life Sciences, Department of Obstetrics and Gynecology,
feasibility of reconstructing genetic lineages and tracing their epigenomic and transcriptomic Third Hospital, Peking University, Beijing 100871, China.
dynamics with single-cell multiomics sequencing. 2
Biomedical Pioneering Innovation Center and Center for
Reproductive Medicine, Ministry of Education Key Laboratory

C
of Cell Proliferation and Differentiation, Beijing 100871,
olorectal cancer (CRC), a major cause of studies have been limited to analysis of bulk cells China. 3Peking-Tsinghua Center for Life Sciences, Academy
mortality, is characterized by heterogeneous which consist of non–tumor cells and complex for Advanced Interdisciplinary Studies, Peking University,
Beijing 100871, China. 4Department of General Surgery,
features of genomic, epigenomic, and tran- subclones and only reflect the average profiles Peking University Third Hospital, Beijing 100191, China. 5Key
scriptomic alterations (1–4), which are not of tumor samples. Single-cell genome and tran- Laboratory of Assisted Reproduction, Ministry of Education,
separate events, as multiple cellular pro- scriptome sequencing have revealed ITH in sev- Beijing 100191, China. 6Beijing Key Laboratory of
cesses may interact to promote tumorigenesis eral cancer types (8–12). However, single-cell Reproductive Endocrinology and Assisted Reproductive
Technology, Beijing 100191, China.
(5, 6). Intratumoral heterogeneity (ITH) across sequencing has been limited in the ability to *These authors contributed equally to this work.
multiple layers of molecular features is a barrier characterize multiple layers of molecular fea- †Corresponding author. Email: jie.qiao@263.net (J.Q.);
for effective diagnosis and treatment (7). However, tures of each genetic lineage. tangfuchou@pku.edu.cn (F.T.); fuwei@bjmu.edu.cn (W.F.)

Fig. 1. Reconstruction of genetic

lineages with scTrio-seq2. Global
SCNA patterns (250-kb resolution)
of CRC01. Each row represents an
individual cell. The subclonal SCNAs
used for identifying genetic sublin-
eages were marked and indexed;
for details, see fig. S6B. On the
top of the heatmap, the ampli-
fication or deletion frequency of
each genomic bin (250 kb) of
the non-hypermutated CRC
samples from the TCGA Project
and patient CRC01’s cancer cells
are shown.

Bian et al., Science 362, 1060–1063 (2018) 30 November 2018 1 of 4

R ES E A RC H | R E PO R T

and CRC11) were assigned to the CMS2 group SCNAs indicate the emergence of sublineages relatively homogenous within a genetic lineage
(fig. S3, A to C), a canonical CRC group with ab- (Fig. 1 and figs. S6 and S7). As whole-chromosome (or sublineage) but showed discrepancies among
normal activation of the WNT/b-catenin and MYC or arm-level events are more likely to be acquired different lineages (or sublineages) (Fig. 2A and
signaling pathways, frequent SCNAs, and non- independently in different lineages (18, 19), we figs. S8C and S9). Tumor lineages’ hypomethylated
hypermutation (2). Whole-genome sequencing mainly used subclonal breakpoints within chro- regions were significantly enriched in long ter-
(WGS) verified a low frequency of somatic single- mosome arms to identify genetic lineages (fig. minal repeats (LTRs), long interspersed nuclear
nucleotide variations (SNVs) in tumors from S6). For five patients for whom we had methyl- element 1 (LINE-1, L1), and heterochromatin re-
CRC01 (fig. S3D) (1) and identified an oncogenic ation data (for >90 cells), cancer cells were clas- gions (H3K9me3) (P < 0.05, Fisher’s exact test); in
mutation of NRAS and inactivating mutations of sified into several genetic sublineages (Fig. 1 and contrast, tumor lineages’ hypermethylated regions
APC and SMAD4, consistent with the features of figs. S5 to S7). For CRC01, we identified 12 sub- were enriched in CpG islands (CGIs), H3K4me3,
non-hypermutated CRC (1). lineages originating from two distinct lineages and open chromatin (P < 0.05, Fisher’s exact test)
For the 10 patients for whom we had DNA (A and B) on the basis of 21 subclonal break- (fig. S10, A and B). A representative long dif-
methylation data, we performed SCNA profiling points; each sublineage was supported by 4 to ferentially methylated region (DMR) (~34 kb)
of individual cells at 250-kilobase (kb) resolution 8 subclonal breakpoints (Fig. 1 and figs. S6 and located in the heterochromatin region of chro-
(Fig. 1 and figs. S4A and S5). Significant focal S7). The sublineage A5 of CRC01 was detected mosome 16 differed between lineages A and B of
SCNAs and probable gene targets were identified in both LN (17%) and ML (87%), indicating that CRC01 (Fig. 2B). This region was hypermeth-
(fig. S4B and table S2) (15). Additionally, the these metastases had a common origin (20). In ylated in NC but heterogeneous in cancer cells
SCNA profiles estimated from the transcriptome all five patients, the primary tumors showed more (fig. S10C).
data were consistent with those estimated from complex subclonal structures than the metastases, We further traced the dynamics of DNA meth-
DNA methylation data (fig. S4C). The scTrio-seq2 indicating that metastases tend to be clonal (fig. ylation and gene expression during metastasis
data confirmed the presence of one cancer cell S7, B and C). within one single lineage each for CRC01 and

Downloaded from http://science.sciencemag.org/ on January 24, 2019

with homozygous deletions of several whole Genome-wide DNA hypomethylation was de- CRC10. Global DNA methylation levels were rel-
chromosomes (fig. S4C). tected in individual cancer cells compared with atively stable during metastasis and accompanied
Genomic alterations in tumors provide mark- paired NC cells (Fig. 2A and figs. S8 and S9A), by changes in focal regions, such as promoters
ers for lineage tracing (16, 17). Clonal variants consistent with findings in published studies (fig. S11, A to C). We did not observe obvious
occur in early tumorigenesis, whereas subclonal (3, 4). Genome-wide DNA methylation levels were changes pre- or postmetastasis in the expression

A CRC01 B CRC01 NC
chr16: 46,385,812 - 46,419,987 (length: ~34 kb)
DNA methylation levels (1-kb tile)

0.8
Lineage A in PT

0.7
Single cells

0.6

0.5

Mean methylation levels

Lineage B in PT

0.4
NC A0 A1 A2 A3 A4 A5 A6 A7 A8 A9 B0 B1
Genetic sub-lineages

Sampling NC PT1 PT2 PT3 PT4 LN1 LN2 LN3 ML1

regions:
ML2 ML3 ML4 MP1 MP2 MP3 MP4 MP5
0.0
0.2
0.4
0.6
0.8
1.0
CpG sites: Methylated Unmethylated
Genetic sub-lineages: NC A0 A1 A2 A3 A4 B0 B1

C Correlations between gene expression and DNA methylation levels

Spearman correlation

CRC01 CRC02 CRC04 CRC09 CRC10 CRC11

0.4

0.2

0.0

−0.2

−0.4

TSS Gene Body TES TSS Gene Body TES TSS Gene Body TES TSS Gene Body TES TSS Gene Body TES TSS Gene Body TES

Fig. 2. Associations between DNA methylation and gene expression Each column shows one single CpG site. The bar plot on the right shows each
levels. (A) Single-cell global DNA methylation levels (1-kb tile) of each cell’s mean methylation levels in the region. (C) The correlations between
sublineage for patient CRC01. Colors represent sampling regions. Lines gene expression levels and DNA methylation levels across gene bodies and
represent median values. (B) One representative DMR on chromosome 16 their 15-kb flanking regions (Spearman correlations). The gray lines
between genetic lineages A and B in the PT sample of CRC01. Each row represents individual cells. The blue line represents the mean value for each
shows one individual cell sorted by the mean methylation levels in the region. patient. TSS, transcription start site; TES, transcription end site.

Bian et al., Science 362, 1060–1063 (2018) 30 November 2018 2 of 4

R ES E A RC H | R E PO R T

A CRC01 Fig. 3. DNA demethylation patterns in cancer cells.

(A) Unsupervised hierarchical clustering of relative
levels in NC
Methylation 0.85
DNA methylation levels in single cancer cells com-
0.80
0.75 pared with NC cells (10-Mb tile). In the upper panel,
0.70 Mean in NC the black points indicate the average DNA methylation
0.65 levels of each genomic bin of normal cells, and the blue
dashed line represents the mean level. (B) Box plot
showing the Pearson correlations between the DNA
Single cells sorted by sub-lineages

demethylation levels of genomic bins in individual

cancer cells and the densities of genomic features
across a range of resolutions.
Lineage A

Downloaded from http://science.sciencemag.org/ on January 24, 2019

B

Genomic bins (10-Mb)

Genetic Relative methylation levels:
sub-lineages: NC A0 A1 A2 A3 A4 A5 A6 A7 A8 A9 B0 B1
0.2
0.4
0.6
0.8
1.0

B CRC01 CRC02 CRC04 CRC09 CRC10

Pearson correlations

0.4
0.2
0.0
−0.2
−0.4 Resolutions:
−0.6 0.25-Mb
0.5-Mb
CRC11 CRC12 CRC13 CRC14 CRC15 1-Mb
2-Mb
Pearson correlations

0.4 5-Mb
0.2 10-Mb
0.0
−0.2
−0.4
−0.6

Regions: 3 L1 n 3 3 L1 n 3 3 L1 n 3 3 L1 tin me3 e3 L1 tin me3

me ati me me ati me me ati me me ma 9m ma
K9 rom K4 K9 rom K4 K9 rom 3K4 K9 hro H3K
4
3K hro H3K
4
H3 nc
h H3 H3 nc
h H3 H3 nc
h H H3 c H c
e e e en en
Op Op Op O p O p

levels or DNA methylation levels of epithelial– proliferation and migration contained DMRs (21) DNA demethylation showed positive cor-
mesenchymal transition–related genes (fig. S11, (fig. S12, C and D). relations with the densities of L1 and H3K9me3
D to G). Furthermore, we did observe molecular Cancer cells experienced variable demethyla- modifications (of NC cells) but negative correla-
associations between the methylome and tran- tion degrees across the whole genome, which tions with the densities of the H3K4me3 modifica-
scriptome of individual cells. Globally, the cor- were concordant within each genetic sublineage tions and open chromatin (of NC cells) (Fig. 3B).
relations between DNA methylation levels and but varied in degree among different sublineages Similar correlations were also observed between
gene expression levels were negative in promoter (Fig. 3A and fig. S13A). The relative demethyla- SNVs and chromatin state in the two patients for
regions but positive in gene body regions in in- tion degrees of cancer cells were correlated with whom there were WGS data (fig. S13C), con-
dividual cells (Fig. 2C). The transcriptome groups the absolute DNA methylation levels of NC cells sistent with chromatin organization influencing
were consistent with genetic lineages and DNA across the genome (fig. S13B), and this compar- regional SNV frequencies (22). Interestingly, L1,
methylation groups for CRC04 but not for CRC01, ison suggested that the regions with higher DNA which is evolutionarily younger and more active
CRC10, and CRC11 (fig. S12, A and B). Additionally, methylation levels in NC cells tended to undergo than LINE-2 (L2) (23), experienced significantly
some promoters of genes associated with tumor stronger demethylation in cancer cells. stronger DNA demethylation than L2 in cancer

Bian et al., Science 362, 1060–1063 (2018) 30 November 2018 3 of 4

R ES E A RC H | R E PO R T

A B C chrX chr1
NC: Lineage B: Ratio compared with NC r22
21 ch
chr

8
6
4
2
0
0.
0.
0.
0.
0.
0

ARID1A
Chr.X r2

FOXO41
AMER6A
KDM
Chr.17

NRAS
Chr.17 ch

EP3 K1
MA
18 ch
19

RU
22 19 r r2

P
00
ch

NX
13

PT
16 22

PR

1
8
4 Amp.

r1

T
19 16

ch
SM
3 15 SM AD H1

r17
1 AD 4
Del. ID B4
7 B

ch
3 20 AX 2 E R
IN Amp.
21 2

chr
1 9

6
TP
53

chr1
2 Del.

3
7 15 CD
H1
8 9 11
AXIN A
10 1 PIK3C

chr15
12 14
SMAD3
5 20 12 B2M
KIT
6 6 Outer: scTrioSeq2

c hr 4
6 AKT1
12 21 7 Inner: TCGA (CRC)

chr14
15 14 3 FOXA1
14 10 21
9 2 10 RB1CA2

chr1
BR AP
16 11 C

5
2

chr
3
22 8 8 3
BB
1 18 5 ER AS PR
KR
D

ch
11 5 18 M

RO
M
r1
1

AT

S1
2
20 13 13

EN 2

Downloaded from http://science.sciencemag.org/ on January 24, 2019

TC A S

r6
EG
PT F7L

CH1

ch
4

HR
17

CDKN2A
4

FR
BRA
FGFR
ch
r1

NOT
19 X

MYC
X 1

F
1
7
DNA methylation levels: Ratio: chr
10 chr
0
5
0
5
0

0
0

0.6
0.6
0.7
0.7

chr9 chr8
0.5

0.6
0.7

Fig. 4. Chromosomal patterns of DNA demethylation and SCNAs of (B) The residual ratios of chromosomal methylation levels of lineage
CRC. (A) The chromosomal DNA methylation levels of NC cells and lineage B B compared with those in NC cells of CRC01. (C) Circos plot showing the
of CRC01. Chromosomes are ranked by their mean methylation levels within amplification or deletion frequency of each genomic bin (250-kb resolution)
each lineage. The lines link the same chromosomes. The six chromosomes of the nine CRC patients in our study. The minimum scale of coordinates
with the lowest DNA methylation levels in lineage B are highlighted in red. is 18 Mb. Some oncogenes and tumor suppressor genes are labeled.

cells for all patients (P < 2.7 × 10−4, Wilcoxon tion and gene expression features. We found that 18. R. Beroukhim et al., Nature 463, 899–905 (2010).
rank-sum test) (fig. S13, D and E). This is in con- DNA methylation profiles are relatively stable 19. T. I. Zack et al., Nat. Genet. 45, 1134–1140 (2013).
trast to DNA demethylation during embryonic within a single genetic lineage or sublineage. The 20. K. Naxerova et al., Science 357, 55–60 (2017).
21. Y. Han et al., Pathology (Phila.) 46, 396–401 (2014).
development, where L1 tends to retain higher differences in DNA methylation levels between 22. B. Schuster-Böckler, B. Lehner, Nature 488, 504–507 (2012).
DNA methylation levels than L2 (24, 25). These the primary tumors and metastases could be 23. R. Cordaux, M. A. Batzer, Nat. Rev. Genet. 10, 691–703 (2009).
results suggest that during tumorigenesis and mainly caused by the differences in the composi- 24. F. Guo et al., Cell 161, 1437–1452 (2015).
progression, L1 and heterochromatin regions ex- tions of sublineages, but not de novo methylation 25. H. Guo et al., Nature 511, 606–610 (2014).
perience aberrant DNA demethylation, breaking or demethylation during metastasis. The DNA
AC KNOWLED GME NTS
the rules for normal development. demethylation patterns of individual cancer cells
We thank I. Bruce for polishing our manuscript. Funding: F.T. was
We further explored the chromosomal patterns were consistent among the 10 CRC patients whose supported by grants from the National Natural Science Foundation of
of aberrant DNA methylation and genome in- DNA we sequenced. Thus, single-cell multiomics China (81561138005, 31625018, and 81521002). This work was
stability of CRC. Chromosomal demethylation sequencing provides insights and resources for supported by the Beijing Advanced Innovation Center for Genomics
degrees were variable, with six chromosomes understanding the molecular alterations that oc- at Peking University. Author contributions: F.T., J.Q., and W.F.
conceived the project. S.B., Y.H., and F.T. wrote the manuscript with
(4, 5, 8, 13, 18, and X) exhibiting stronger DNA cur during CRC progression and metastasis. help from all authors. Y.H., X.Z., X.L., Jun.Y., Y.W., W.W., Jia.Y., H.G.,
demethylation than others (Figs. 4, A and B, S.G., and Y.M. performed the experiments. S.B., X.L., B.H., J.D., and
and fig. S14A). These six chromosomes were also RE FERENCES AND NOTES P.Z. conducted the bioinformatics analyses. Competing interests: The
significantly enriched for stronger DNA demeth- authors declare they have no competing interests. Data and materials
1. Cancer Genome Atlas Network, Nature 487, 330–337 (2012).
availability: The scTrio-seq2 data have been submitted to the NCBI
ylation in most patients (P < 0.05, hypergeo- 2. J. Guinney et al., Nat. Med. 21, 1350–1356 (2015).
Gene Expression Omnibus (GEO) under accession number GSE97693.
metric test) (fig. S14B). Additionally, population 3. T. Hinoue et al., Genome Res. 22, 271–282 (2012).
WGS data have been deposited in the European Genome-phenome
4. B. P. Berman et al., Nat. Genet. 44, 40–46 (2012).
analyses of the non-hypermutated CRC samples 5. P. A. Jones, S. B. Baylin, Cell 128, 683–692 (2007).
Archive (EGA) under the accession number EGAS00001003242.
from the TCGA Project and our study showed Analysis code is available at https://github.com/bianshuhui/CRC_code.
6. J. S. You, P. A. Jones, Cancer Cell 22, 9–20 (2012).
that three of the hypomethylated chromosomes 7. N. McGranahan, C. Swanton, Cancer Cell 27, 15–26 (2015).
8. N. E. Navin, Genome Res. 25, 1499–1507 (2015). SUPPLEMENTARY MATERIALS
recurrently generated SCNAs (chromosomes 8,
9. I. Tirosh et al., Science 352, 189–196 (2016). www.sciencemag.org/content/362/6418/1060/suppl/DC1
13, and 18) (Fig. 4C and fig. S14C). Using WGS 10. I. Tirosh et al., Nature 539, 309–313 (2016). Materials and Methods
data, we found that five of the hypomethylated 11. A. S. Venteicher et al., Science 355, eaai8478 (2017). Figs. S1 to S14
chromosomes were also significantly enriched 12. H. Li et al., Nat. Genet. 49, 708–718 (2017). Tables S1 and S2
for SNVs (chromosomes 4, 5, 8, 13, and X) (P < 13. Y. Hou et al., Cell Res. 26, 304–319 (2016). References (26–29)
14. S. A. Smallwood et al., Nat. Methods 11, 817–820 (2014).
0.05, hypergeometric test) (fig. S14D). 15. C. H. Mermel et al., Genome Biol. 12, R41 (2011). 13 July 2017; resubmitted 24 August 2018
Unambiguously assigned cancer cell lineages 16. R. Gao et al., Nat. Genet. 48, 1119–1130 (2016). Accepted 25 October 2018
enable characterization of their DNA methyla- 17. N. Navin et al., Nature 472, 90–94 (2011). 10.1126/science.aao3791

Bian et al., Science 362, 1060–1063 (2018) 30 November 2018 4 of 4

Single-cell multiomics sequencing and analyses of human colorectal cancer
Shuhui Bian, Yu Hou, Xin Zhou, Xianlong Li, Jun Yong, Yicheng Wang, Wendong Wang, Jia Yan, Boqiang Hu, Hongshan
Guo, Jilian Wang, Shuai Gao, Yunuo Mao, Ji Dong, Ping Zhu, Dianrong Xiu, Liying Yan, Lu Wen, Jie Qiao, Fuchou Tang and
Wei Fu

Science 362 (6418), 1060-1063.

DOI: 10.1126/science.aao3791

scTrio-seq identifies colon cancer lineages

To better design treatments for cancer, it is important to understand the heterogeneity in tumors and how this
contributes to metastasis. To examine this process, Bian et al. used a single-cell triple omics sequencing (scTrio-seq)

Downloaded from http://science.sciencemag.org/ on January 24, 2019

technique to examine the mutations, transcriptome, and methylome within colorectal cancer tumors and metastases from
10 individual patients. The analysis provided insights into tumor evolution, linked DNA methylation to genetic lineages,
and showed that DNA methylation levels are consistent within lineages but can differ substantially among clones.
Science, this issue p. 1060

ARTICLE TOOLS http://science.sciencemag.org/content/362/6418/1060

SUPPLEMENTARY http://science.sciencemag.org/content/suppl/2018/11/28/362.6418.1060.DC1
MATERIALS

RELATED http://stm.sciencemag.org/content/scitransmed/7/290/290ra91.full
CONTENT
http://stm.sciencemag.org/content/scitransmed/8/360/360ra135.full
http://stm.sciencemag.org/content/scitransmed/6/233/233re2.full
http://stm.sciencemag.org/content/scitransmed/8/361/361ra140.full
REFERENCES This article cites 28 articles, 4 of which you can access for free
http://science.sciencemag.org/content/362/6418/1060#BIBL

PERMISSIONS http://www.sciencemag.org/help/reprints-and-permissions

Use of this article is subject to the Terms of Service

Science (print ISSN 0036-8075; online ISSN 1095-9203) is published by the American Association for the Advancement of
Science, 1200 New York Avenue NW, Washington, DC 20005. 2017 © The Authors, some rights reserved; exclusive
licensee American Association for the Advancement of Science. No claim to original U.S. Government Works. The title
Science is a registered trademark of AAAS.