WO2023157933A1 - 細胞の系統解析を行う方法 - Google Patents

細胞の系統解析を行う方法 Download PDF

Info

Publication number
WO2023157933A1
WO2023157933A1 PCT/JP2023/005534 JP2023005534W WO2023157933A1 WO 2023157933 A1 WO2023157933 A1 WO 2023157933A1 JP 2023005534 W JP2023005534 W JP 2023005534W WO 2023157933 A1 WO2023157933 A1 WO 2023157933A1
Authority
WO
WIPO (PCT)
Prior art keywords
cell
cells
mutation
sites
somatic
Prior art date
Application number
PCT/JP2023/005534
Other languages
English (en)
French (fr)
Inventor
聡史 太田
Original Assignee
国立研究開発法人理化学研究所
Priority date (The priority date is an assumption and is not a legal conclusion. Google has not performed a legal analysis and makes no representation as to the accuracy of the date listed.)
Filing date
Publication date
Application filed by 国立研究開発法人理化学研究所 filed Critical 国立研究開発法人理化学研究所
Publication of WO2023157933A1 publication Critical patent/WO2023157933A1/ja

Links

Classifications

    • G—PHYSICS
    • G16—INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR SPECIFIC APPLICATION FIELDS
    • G16B—BIOINFORMATICS, i.e. INFORMATION AND COMMUNICATION TECHNOLOGY [ICT] SPECIALLY ADAPTED FOR GENETIC OR PROTEIN-RELATED DATA PROCESSING IN COMPUTATIONAL MOLECULAR BIOLOGY
    • G16B10/00—ICT specially adapted for evolutionary bioinformatics, e.g. phylogenetic tree construction or analysis

Definitions

  • the present invention relates to a method, system and program for performing lineage analysis of cells.
  • Somatic mutations accumulate during development and aging, causing somatic genome mosaicism. Somatic mutations are associated with various diseases such as dementia, cardiovascular disease, and cancer (Non-Patent Documents 1 and 2) and aging (Non-Patent Document 3). Furthermore, somatic mutations are also known to be involved in normal development, such as the acquisition of neuronal population diversity in the central nervous system (Non-Patent Documents 4 and 5). In addition, it is known from RNA sequence data analysis that somatic mutations spread clonally in normal tissues (Non-Patent Document 6).
  • Non-Patent Documents 7, 8 2.8 (95% confidence interval, 2.4-3.3) somatic mutations occur per cell replication. This is slightly higher than the incidence of mutations in germ cells (germline mutations) (Non-Patent Documents 7, 8). However, later in development, 6.4 ⁇ 10 ⁇ 10 to 7.8 ⁇ 10 ⁇ 10 somatic mutations occur per base pair per cell division in various cell types (Non-Patent Document 8). , 9). These can be up to 10 times more germline mutations (Non-Patent Documents 10, 11). As a result, even monozygotic twins may become genetically diverse due to somatic mutation (Non-Patent Document 12).
  • somatic mutation is an event that occurs in various places and times in various organisms, including humans (Non-Patent Document 13). Somatic mutations can be considered as counterparts to germline mutations (Non-Patent Document 14), whose dynamics are subject to evolutionary processes (Non-Patent Document 15). Like germ cells, somatic cells have an 'evolutionary record' as 'scars' in the genome over time, and their history can theoretically be restored by retrospective cell tracking methods (Non-Patent Documents 16, 17). . In other words, somatic mutations have cell lineage temporal information.
  • Non-Patent Document 18 For example, especially in the field of oncology, somatic mutations have been extensively studied in relation to cancer evolution (Non-Patent Document 18).
  • researchers have devised mathematical models to elucidate the dynamics of somatic mutations from various biological perspectives (Non-Patent Documents 19 and 20).
  • Non-Patent Documents 32, 33 t-distributed stochastic neighborhood embedding
  • UMAP uniform manifold estimation and projection for dimension reduction
  • the present inventors focused on somatic mutations identified by comparing transcriptome sequence data with reference genome sequence data, and found that The inventors found that phylogenetic analysis was performed, and conceived of the present invention.
  • the present invention focuses on phylogenetic signatures of single-cell transcriptomes. For example, 60 somatic cells have about 10 96 possible lineages, which is more than the number of atoms in the universe [31].
  • the method of the present invention takes advantage of the coarse phylogenetic signal carried in RNA-sequence data to narrow the window covering true lineages.
  • sources of RNA sequencing have available hundreds or thousands of copies of the RNA molecule within the cell. False positives due to biological and technical factors such as RNA editing, sequencing errors (e.g. random errors that occur during reverse transcription and PCR), and potential sampling errors should be considered [13, 25]. ).
  • single-cell transcriptome data are overwhelmingly abundant and compensate for these disadvantages of RNA-sequencing data (26).
  • a method for performing lineage analysis of cells comprising: (1) For a plurality of single cells derived from the same individual, compare the RNA sequence data of the trunk liptome of each single cell with the corresponding genome sequence data, and screen for sites where the RNA sequence and the genome sequence differ; (2) Of the sites identified in the screening in step (1), a site where three or more single nucleotide polymorphisms are detected is screened as a site of somatic mutation; and (3) creating a cell lineage tree based on the nucleotide sequence information of the site of somatic mutation obtained in step (2);
  • the method as described above comprising: [2] (4) The method of [1], further comprising estimating the cell type of each single cell from the cell lineage tree.
  • [3] The method of [1] or [2], further comprising comparing the information of the single cell type deduced from the gene expression profile of each single cell with the information of the cell lineage tree.
  • [4] The method according to any one of [1] to [3], wherein the single nucleotide polymorphism is a single nucleotide substitution.
  • step (1) in 50% or more of all cells, sites for which trunk cryptome RNA sequence data cannot be obtained are presumed to have no mutation and are excluded from screening [1]- The method according to any one of [4].
  • [6] The method according to any one of [1] to [5], wherein the somatic mutation screened in step (2) is a rare somatic mutation.
  • [7] The method of any one of [1]-[6] for detecting a somatic mutation associated with a disease or condition.
  • a recording medium recording a program for executing the method according to any one of [1] to [8].
  • a system comprising a processor and a memory storing a program that, when executed by the processor, performs the method of any one of [1]-[8].
  • the method of the present invention utilizes latent somatic mutations detected from single-cell transcriptome data to provide a new framework for detecting cell lineage trajectories. Unlike pseudo-time course analysis, in which heterogeneously differentiated cells are analyzed in a single time snapshot, the present invention can trace cell lineages back to ancestral cells using somatic mutations.
  • FIG. 1 shows the pathway for single-cell somatic mutation pattern analysis.
  • FASTQ (78), BAM (4) and VCF (79) files were generated for each single cell.
  • FASTQ, BAM, and VCF are data formats for describing genetic information.
  • FIG. 2 shows detection of somatic mutations by mutants with multiple mutation sites.
  • the mutational context of this cell is distinguishable from (a) and (b) by a third type of nucleotide G.
  • FIG. 3 shows an evaluation of the degree of monophyleticity (DoM) at each node as a percentage of clusters in the subtree. Map the gene expression profiles (table of clusters) to a tree topology. If the DoM and gene expression profiles are in perfect agreement with each other, a set of subtrees can be observed that are monophyletic in terms of clusters. If the DoM and gene expression profiles do not match perfectly, paraphyletic subtrees or multiphyletic subtrees will be observed in terms of clusters.
  • DoM degree of monophyleticity
  • FIG. 6 is a mapping of gene expression profiles to a cell lineage tree (SRP090944 data; batch 2, 23 cells).
  • FIG. 7 is a mapping of gene expression profiles re-analyzed by t-SNE to a cell lineage tree (SRP090944 data; batch 1, 54 cells).
  • FIG. 8 shows the results of a pseudo-time course analysis of SRP090944 data; Batch 1.
  • FIG. 9 relates to a model system linking the function of Notch/Wnt signaling to cytotrophoblast and extravillous trophoblast differentiation (65).
  • vCTB villous cytotrophoblast
  • CCT column trophoblast
  • EVT extravillous trophoblast
  • N1 ICD Notch1 intercellular domain
  • TEAD4 transcription enhancer factor TEF-3 gene
  • IRF6 interferon regulatory factor 6
  • Notch2 Notch receptor 2.
  • the present invention includes the following aspects. Unless defined otherwise herein, technical and scientific terms used herein have the same meaning as commonly understood by one of ordinary skill in the art.
  • the materials, materials and examples disclosed herein are exemplary only and are not intended to be limiting. When the specification refers to "in one embodiment,” it is meant to be non-limiting, ie, non-limiting, to that embodiment.
  • the present invention relates to methods of performing lineage analysis of cells.
  • the method includes: (1) For a plurality of single cells derived from the same individual, compare the RNA sequence data of the trunk liptome of each single cell with the corresponding genome sequence data, and screen for sites where the RNA sequence and the genome sequence differ; (2) Of the sites identified in the screening in step (1), a site where three or more single nucleotide polymorphisms are detected is screened as a site of somatic mutation; and (3) creating a cell lineage tree based on the nucleotide sequence information of the site of somatic mutation obtained in step (2); Including.
  • the method of the invention utilizes multiple single cells derived from the same individual.
  • the type of individual is not particularly limited. Non-limiting examples include humans, non-human primates (monkeys, gorillas, chimpanzees, etc.), mice, rats, and the like.
  • multiple single cells from the same tissue or from adjacent tissues of the same individual are utilized. Adjacent tissues include, but are not limited to, tissues that are physically close to each other in vivo (e.g., small intestine and large intestine, placenta and fetus, etc.), functionally related tissues (e.g., digestive system, respiratory system, central nervous system, etc.). system, etc.).
  • RNA sequence data of the trunk liptome of each single cell is not particularly limited.
  • data published in data banks such as DDBJ, EMBL-Bank/EBI, GenBank/NCBI, etc.
  • data published in documents such as Non-Patent Document 36 can be used.
  • the source of the corresponding genome sequence data to be compared with the trunk cryptome RNA sequence data is not particularly limited.
  • data published in data banks such as DDBJ, EMBL-Bank/EBI, and GenBank/NCBI can be used.
  • the step of comparing the RNA sequence data of the trunk liptome of each single cell with the corresponding genomic sequence data and screening for sites where the RNA sequence and the genomic sequence differ uses known software for comparing and mapping sequences. can be done by
  • mapping can use Burrows-Wheeler Aligner (BWA).
  • BWA is a software package for mapping low-range sequences against a large reference genome, such as the human genome, and includes two algorithms: BWA-backtrack, BWA-SW and BWA-MEM.
  • Other software such as BarraCUDA, STAR, etc. may be used to map the trunk glyptomic RNA sequence data to the corresponding genomic sequence data.
  • an adapter sequence may be added to the RNA sequence for analysis. In that case, it is preferable to remove the adapter sequences prior to mapping the RNA sequence data to the genomic sequence data.
  • RNA sequence differs from the genome sequence can be performed using publicly known software for polymorphism detection and alignment display.
  • Samtools may be used.
  • Samtools is software used for post-short-read sequence alignment processing in DNA sequencing.
  • the main functions of Samtools include polymorphism detection, alignment display, indexing, data extraction, and file format conversion.
  • software such as elPrep and Picard may be used for screening sites where the RNA sequence differs from the genomic sequence.
  • step (2) among the sites identified in the screening in step (1), sites where three or more single nucleotide polymorphisms are detected are screened as sites of somatic mutation.
  • Non-Patent Document 8 the in vivo tissue-specific somatic mutation probability in normal lymphocytes, precursors of chronic lymphocytic leukemia (CLL), is 7.6 per base per cell division. ⁇ 10 ⁇ 10 ⁇ 1.1 ⁇ 10 ⁇ 10 (SE). Therefore, stochastic somatic changes mean that an average of 3 mutations occur per cell division in the human genome, excluding repetitive regions (Non-Patent Documents 70, 71). The abundance of such somatic cell mutations makes it possible to create cell lineages that trace back to the past even in normal tissues.
  • Three single nucleotide polymorphisms means three types of nucleotide polymorphisms. If a locus (here, a coordinate on the genome) is homozygous, it can be determined whether a mutation occurred at that locus by comparison with the standard genomic sequence. However, if a locus is heterozygous, it cannot be immediately determined whether it is the result of a mutation or is originally heterozygous. That is, when there are two kinds of mutations, it is difficult to determine whether or not mutation occurs at the gene locus. However, in the case of three mutations, this observation cannot be explained without assuming at least one mutation. In the present invention, based on this logic, "sites where three or more single nucleotide polymorphisms are detected" are targeted for mutation detection.
  • the number of polymorphisms in the "site where three or more single nucleotide polymorphisms are detected there is no particular upper limit on the number of polymorphisms in the "site where three or more single nucleotide polymorphisms are detected". For convenience of analysis, it is preferable not to have too many polymorphisms. Without limitation, the number of polymorphisms is 15 or less, 12 or less, 10 or less, 8 or less, 6 or less, 5 or less, 4 or less. In one aspect, the number of polymorphisms is three.
  • the "single nucleotide polymorphism" used in the present invention means a general mutation to another base (group) that occurs at the single base level, and may be any of base (group) substitution, insertion (addition), or deletion. .
  • the term "mutation occurring at the single base level” includes not only single base mutation (point mutation) but also mutation occurring at multiple base units.
  • the single nucleotide polymorphism is a single nucleotide mutation (point mutation).
  • the single nucleotide polymorphism is a single nucleotide insertion and/or substitution.
  • the single nucleotide polymorphism is a single nucleotide substitution.
  • step (1) sites for which no RNA sequence data of the trunk liptome can be obtained in 50% or more of all cells are presumed to have no mutation and are excluded from screening.
  • the method utilizes RNA sequence data of the trunk cryptome. Unlike a static genome, the amount of transcripts (transcripts) fluctuates from cell to cell and over time. Therefore, depending on the cell, there may be cases (sites) in which RNA of a specific transcript is not detected and RNA sequence data cannot be obtained.
  • trunk cryptome RNA sequence data is obtained in 50% or more of all cells, that is, transcript sequences mapped to 50% or more cells are used for mutation detection. do.
  • step (1) sites for which RNA sequence data of the trunk liptome cannot be obtained in 50% or more of all cells are presumed to have no mutation and are excluded from screening. "Assuming no mutation occurred” means substituting the corresponding genomic sequence. "50% or more" of "in 50% or more of all cells” is, without limitation, 55% or more, 60% or more, 65% or more, 70% or more, 75% or more, 80% or more may
  • a threshold may be selected for excluding potential mutations in the screening of the present invention.
  • a threshold of 80% was used as an example solution. Specifically, multiple mutation sites observed in 80% or more of single cells tested were screened. Without limitation, higher thresholds, such as 90% or greater, or lower thresholds, such as 60% or greater, 70% or greater, may be employed.
  • the somatic mutation screened in step (2) is a rare somatic mutation.
  • Infrequent somatic mutations refer to less frequent somatic mutations in transcripts observed in cell populations.
  • Allyl originally means allele, but in the present specification, it may be used to extend the meaning to single-nucleotide mutation.
  • the concept of allele originally belongs to the genome sequence, but assuming that there is some correlation between the allele frequency according to the transcript sequence and the allele frequency according to the genome sequence, alleles with lower frequency in the genome sequence, that is, Minor alleles are candidates for somatic mutations to be screened.
  • Non-Patent Documents 74, 75 The assumption that there is some degree of correlation between allele frequencies from transcript and genome sequences is based on the empirical expectation that allele frequencies of genome and transcriptome reads are correlated.
  • Non-Patent Documents 74, 75 In order to reduce the error, it is preferable to further apply a more quantitative approach. For example, statistical haplotyping phasing [76] is believed to be effective for improving estimates.
  • the frequency of mutations found in alleles can be estimated, for example, by comparing the reference genome and the lead sequence.
  • the meaning of "low frequency” is non-limiting, for example, within 50% of the low frequency of all mutations found in alleles, preferably within 40%, within 30%, within 20%, within 10 % or less means within 5%.
  • the meaning of "low frequency” means within 10th, within 8th, within 5th, within 3rd of all mutations found in the allele.
  • a rare somatic mutation refers to the lowest frequency somatic mutation among all mutations found in an allele.
  • step (3) a cell lineage tree is created based on the base sequence information of the somatic mutation sites obtained in step (2).
  • Phylip For creating a cell phylogenetic tree, known software can be used to create a phylogenetic tree from sequence information. Some MEGA-X can be used. Alternatively, software such as Phylip may be used.
  • the maximum parsimony method may be applied in the creation of the cell lineage tree.
  • the maximum parsimony method is a kind of method for creating a phylogenetic tree by finding a mutation pattern that minimizes the number of mutations. Similar analysis can also be performed using other methods such as the distance matrix method, maximum likelihood method, Bayesian method, and the like.
  • the method may further include estimating the cell type of each single cell from the cell lineage tree.
  • Estimation (identification) of the cell type of each cell from the cell lineage tree can be performed using a known method.
  • known tools such as principal component analysis, t-SNE, and UMAP may be used.
  • a gene expression matrix consisting of single cell types and gene types may be subjected to linear or nonlinear dimensionality reduction to cluster cells, or information on known marker genes may be used.
  • Non-Patent Document 45 is a tool that annotates predicted mutations and provides information on their impact. It is possible to obtain information on the number of each mutation, positional information/type (SNP, insertion, deletion, etc.), and degree of influence (substitution to amino acid).
  • the mutation site can be annotated based on the annotation assigned to the genome, depending on where in the genome the mutation site is located. However, since the data used in the above method may be transcript sequences, there may be some uncertainty.
  • the SnpEff software can take into account the circumstances and adopt annotations that are considered appropriate. Furthermore, these annotations allow biological and evolutionary interpretation of each mutation.
  • the method may further comprise comparing single cell type information deduced from the gene expression profile of each single cell with information in the cell lineage tree.
  • the method can identify single cell types that could not be identified in gene expression profiles.
  • gene expression profiles can review and correct misidentified single cell types.
  • the method can be used to detect somatic mutations associated with a disease or condition.
  • the type of disease or symptom is not particularly limited.
  • the disease or condition is a disease or condition caused by genetic mutation.
  • the disease or condition is selected from the group consisting of cancer, dementia, cardiovascular disease, aging, autoimmune disease, neurodegenerative disease, and psychiatric disease.
  • the present invention provides a new framework for estimating cell lineage trees using somatic mutations detected from single-cell transcriptomes of the cover. In doing so, we focused on systematic features of single-cell transcriptomes rather than individual mutations. The results showed that it was possible to reconstruct a cell lineage tree consistent with known biological knowledge.
  • the significance of the present invention is to infer a global phylogenetic tree of a single cell and provide supporting information for interpreting gene expression profiles. In gene expression analysis, dimensionality reduction methods are often used. However, how to extract biologically relevant findings from the results obtained at the single-cell level remains an open question (Non-Patent Document 69). This framework allows us to interpret gene expression profiles from different angles of single-cell mosaicism.
  • the real-time course of the present approach can take advantage of somatic mutations to trace cell lineages back to progeny cells. . That is, the method of the present invention makes it possible to infer the time course of cells.
  • “Time course speculation” is speculation about unobservable past events, etc., including, for example, stem cell self-renewal stages (Figs. 5c and 6a).
  • the real-time course of the present invention can yield results that closely match those of the pseudo-time-course analysis, but with a much different kind of data, namely mapped read sequences and their depth, which allows for more granularity. It is possible to obtain results in a granular way.
  • programs, Systems, Etc. for Performing Cell Lineage Analysis In order to efficiently execute the cell lineage analysis method of the present invention, it is preferable to create a program for automatically executing the method by a computer. . Also, this program may be recorded on a recording medium so as to be read by a computer. Furthermore, it may be a dedicated system in which a program is recorded in memory so that it can be executed by a processor of a computer.
  • the invention relates to a program for carrying out the method of the invention.
  • the present invention is a program for executing a method for lineage analysis of cells, the method comprising: (1) A step of comparing the RNA sequence data of the trunk liptome of each single cell with the corresponding genome sequence data for a plurality of single cells derived from the same individual, and screening for sites where the RNA sequence differs from the genome sequence; (2) A step of screening a site where three or more single nucleotide polymorphisms are detected among the sites identified in the screening of (1) as a site of somatic mutation; and (3) The above program, which includes the step of creating a cell lineage tree based on the nucleotide sequence information of the site of somatic mutation obtained in (1).
  • the present invention relates to a recording medium recording a program for executing the method of the present invention.
  • the invention relates to a system comprising a processor and a memory having stored thereon a program which, when executed by said processor, performs the method according to the method of the invention.
  • the present invention provides a system for performing lineage analysis of cells, comprising: (1) A means of screening for sites where the RNA sequence differs from the genome sequence by comparing the RNA sequence data of the trunk liptome of each single cell with the corresponding genome sequence data for a plurality of single cells derived from the same individual; (2) Means for screening the site where three or more single nucleotide polymorphisms are detected among the sites identified in the screening of (1) as the site of somatic mutation; and (3) The above-described system, which includes means for creating a cell lineage tree based on the base sequence information of the somatic mutation site obtained in (1).
  • Lineage analysis of cells "(1) For multiple single cells derived from the same individual, compare the RNA sequence data of the trunk liptome of each single cell with the corresponding genome sequence data, screen for different sites", "(2) Among the sites identified in (1), the site where three or more single nucleotide polymorphisms are detected is screened as the site of somatic mutation ", "(3) Create a cell lineage tree based on the nucleotide sequence information of the site of somatic mutation obtained in (1)” ” is as described in the item. All of the items described in “1. Method for Performing Lineage Analysis of Cells” also apply to systems, programs, and recording media.
  • Example 1 Mapping of transcriptome sequence data and detection of mutations (1-1) Transcriptome data
  • SRP090944 batch 1 54 cells
  • batch 2 33 cells
  • Non-Patent Document 36 placental data were analyzed from the perspective of cellular communication networks between two semiallogenic individuals, a mother and a fetus.
  • the intercellular interactome (Non-Patent Document 37) was estimated from the gene expression of receptor-ligand pairs across cell types.
  • DDBJ Sequence Read Archive (DRA) as SRS1732266 (SRX2225269)-SRS1732319 (SRX2225328).
  • Pavlicev et al. analyzed single-cell data by gene expression profiles and principal component analysis (PCA) on 300 marker genes into cytotrophoblast (CYT) 1, CYT2, CYT3, extravillous trophoblast ( EVT), divided into 5 clusters (cell types) of maternal decidual cells (DC)
  • PCA principal component analysis
  • CYT cytotrophoblast
  • EVT extravillous trophoblast
  • DC maternal decidual cells
  • Pavlicev et al. reported the combined expression of /CD4+/CD83+/CD86+ markers and determined to be uterine dendritic cells, but some of the genes characteristic of DCs, such as CLEC4C, THBD, No expression of CD1C, CD80, IL10, IL12B, etc. was observed.
  • FIG. 1 shows an overview of the data analysis pipeline of this example.
  • Single-cell transcriptome sequence data were mapped to the human genome (GRCh38) using Burrows-Wheeler Aligner (BWA) 40 after removal of adapter sequences with trimomatic (41). did.
  • BWA Burrows-Wheeler Aligner
  • FIG. 2 shows a conceptual diagram of one embodiment of the method of the present invention. The method of the present invention includes the following steps.
  • step (1) For a plurality of single cells, compare the trunk liptome RNA sequence data of each single cell with the corresponding genome sequence data, and screen for sites where the RNA sequence and the genome sequence differ; (2) The data of the sites identified in the screening of step (1) are compared between cells, and sites where three or more single nucleotide polymorphisms are detected are screened as sites of somatic cell mutation.
  • Non-Patent Document 45 a VCF file obtained by mapping the transcript sequence to the genome was given as input data, and the annotation was added by designating the used genome data.
  • the average coverage was 0.685% (SD: 0.231) (Fig. 4, batch 1 cells).
  • the initial screen detected 1,965,629 sites that differed from the reference genome.
  • 89 polymutation sites where three or more types of nucleotides were mutated were observed in all 54 single cells.
  • Tables 1 and 2 include alternative annotations such as nested intron genes (Non-Patent Document 63). Some mutation numbers overlap between categories. For example, 2083 and 574 mutations were detected from Batch 1 and Batch 2 data, respectively, on an 80% basis. However, the SnpEff software estimates 1903 and 1398 mutations from the Batch 1 and Batch 2 data, respectively, with the default parameter set. The fact that the orders of the numbers of both mutation sites estimated by independent software are almost the same indicates the validity of the annotation used in this analysis.
  • Example 2 Phylogenetic Analysis of Cells and Creation of Cell Phylogenetic Tree (2-1) Phylogenetic Analysis of Cells All observable mutation sites obtained in Example 1 were linked to create a sequence alignment. "Mutation site ligation" means joining the sequence of the coding region flanking the mutation with the mutation site to create a 3-base codon to create a codon sequence containing the mutation. Cell lineage trees were reconstructed with default parameters using the maximum parsimony method (46, 47) implemented in MEGA X (48).
  • the created multiple alignment was given to MEGA X as input data, and the phylogenetic tree was reconstructed by selecting the maximum parsimony method from the GUI.
  • the results were output in the Newick tree format (Non-Patent Document 49) and used for subsequent processing.
  • Non-Patent Documents 32, 33 t-SNE (Non-Patent Documents 32, 33) and UMAP (Non-Patent Documents 34, 35) were applied to perform nonlinear dimensionality reduction (NDR).
  • the Louvain method (Non-Patent Document 55) was used for clustering. Specifically, after performing global linear dimensional compression by principal component analysis, gene expression matrices were applied to each software as input data.
  • Mathematica is mathematical formula manipulation software.
  • Physicals for Mathematica is a library for mathematical formula manipulation software.
  • AssignCluster2Cell is the name of the program developed to integrate the cell lineage tree and the gene expression profile in this example.
  • AssignCluster2Cell reads the reconstructed cell lineage tree in Newick format, compares it with the pre-clustered cell population, and displays the degree of agreement on each node as a pie chart.
  • the library of functions provided by Phylogetics for Mathematica and Phylogetics is used to load and visualize the tree.
  • AssignCluster2Cell itself is a function written in Mathematica and optimized to run on Mathematica Notebook.
  • Non-Patent Document 62 Standard data format for expressing phylogenetic trees
  • DoM degree of monophyleticity
  • the pie chart in FIG. 3 represents the DoM of each node. For example, if there is only one pie chart for node n, subtree n is completely monophyletic, and if there are two pie charts for node n', subtree n' is polyphyletic (subtree B in FIG. 3). I understand. Since cell type 2 predominates in subtree B, we speculate that cell type 1 was derived from cell type 2. Thus, we delineated cell type associations based on gene expression profiles and cell lineages.
  • the root of the phylogenetic tree represents the zygotic cell (fertilized egg), and its instance represents the zygotic genome of the fertilized egg.
  • the root of the deduced cell lineage tree may represent the progenitor cells of the observed cell population. Junctional cells lie somewhere between the root of the phylogenetic tree and the reference genome.
  • Cytotrophoblast (CYT: CYT1, CYT2, CYT3 in FIGS. 3 and 4) and extravillous trophoblast (EVT) cell lineages have been successfully modeled in the context of Notch/Wnt signaling (65).
  • EVT differentiates from CYT via cell column trophoblasts (CCT) (Fig. 9).
  • the results of this example are consistent with the model of Notch/Wnt signaling described in Non-Patent Document 65.
  • Pavlicev et al. state that their data include maternal decidual cells (DC).
  • the results of this example suggest that the putative DC cells are differentiated from fetal CYT or CYT stem cells (FIGS. 5 and 6).
  • reanalysis of transcriptome data using t-SNE (Non-Patent Documents 32, 33) also supported the correctness of the results of this example (Fig. 7).
  • Non-Patent Documents 50, 51 the maximum likelihood method suggests that overall somatic mutations are subject to purging selection, suggesting that the majority of detected somatic mutations are at least far from positive false.
  • Example 3 Evaluation of the dN/dS Ratio in the Coding Region
  • the dN/dS ratio in the coding region was evaluated.
  • Non-Patent Documents 50, 51 In normal tissues, genuine somatic mutations are expected to undergo purification selection (Non-Patent Documents 50, 51).
  • the selective pressure of the mutants detected in Example 1 was evaluated in terms of the dN/dS ratio.
  • the codon sequences containing the detected variants were assembled to generate a codon alignment including the exon variants.
  • the overall dN/dS ratio was calculated using Codeml from the Paml package [52].
  • Comparative Example 1 Pseudo Time Course Analysis
  • monocle3 version 1.0.0
  • Non-Patent Document 56-58 version 4.1.2
  • R version 4.1.2
  • SRP090944 A pseudo-time course analysis of data batch 1 was performed. Specifically, a dimensionally compressed (26-dimensional) single-cell gene expression matrix was provided as input data, and each single cell was assigned a relative pseudo-time expected to reflect differentiation aspects.
  • Cellular pseudo-time is represented by contour lines generated with interpolation order 3 by Mathematica's ListContourPlot function [59].
  • the results were fairly consistent with the results using the method of the present invention, except for the direction of pseudo-time.
  • the number of cells with a pseudo time greater than 2.5 is 16 568, SRR4371569, SRR4371570, SRR4371571, SRR4371572, SRR4371575 and SRR4371577), Met.
  • RNA-seq data can be used to reconstruct the somatic tree of human placental tissue.
  • This phylogenetic tree is consistent with the known placental cell lineages for four cell types: cytotrophoblast (CYT) I, CYT II, CYT III, and extravillous trophoblast (EVT).
  • CYT cytotrophoblast
  • EVT extravillous trophoblast
  • the method of the present invention uses somatic mutations to trace cell lineages back to progeny. Inferred progeny cells are therefore represented as internal nodes (vertices) of the reconstructed celler trajectory tree.
  • the quality of the detected somatic mutations was evaluated in an evolutionary framework. It was confirmed that the detected somatic mutations underwent purifying selection (dN/dS ⁇ 1) as a whole, but also showed signs of adaptive evolution (dN/dS>1), at least at certain sites. Ta. These evolutionary insights support moderate confidence in our results, at least in the coding region. Furthermore, the results of the Examples were consistent with the results of pseudo-time course analysis using a completely different type of data: sequence and depth of mapped reads.
  • RNA sequence data is a "by-product", and polymorphism (mosaic) information between single cells has often been ignored.
  • the present invention has shown that it is possible to extract significant genotype information by a gleaning method.
  • the present invention's somatic mutational analysis of single-cell transcriptome data enables the interpretation of high-dimensional gene expression data. Moreover, the biological significance of somatic mutation sheds light on a new perspective of 'evolution' within individuals.

Landscapes

  • Physics & Mathematics (AREA)
  • Health & Medical Sciences (AREA)
  • Life Sciences & Earth Sciences (AREA)
  • Animal Behavior & Ethology (AREA)
  • Biophysics (AREA)
  • Physiology (AREA)
  • Engineering & Computer Science (AREA)
  • Bioinformatics & Cheminformatics (AREA)
  • Bioinformatics & Computational Biology (AREA)
  • Biotechnology (AREA)
  • Evolutionary Biology (AREA)
  • General Health & Medical Sciences (AREA)
  • Medical Informatics (AREA)
  • Spectroscopy & Molecular Physics (AREA)
  • Theoretical Computer Science (AREA)
  • Measuring Or Testing Involving Enzymes Or Micro-Organisms (AREA)

Abstract

本発明は、細胞の系統解析を行う方法に関する。本発明の方法は、（１）同一個体由来の複数の単一細胞について、各単一細胞のトランクリプトームのＲＮＡ配列データを、対応するゲノム配列データと比較し、ＲＮＡ配列とゲノム配列が異なる部位をスクリーニングし；（２）工程（１）のスクリーニングで特定された部位のうち、３つ若しくはそれ以上の１塩基多型が検出される部位を、体細胞変異の部位としてスクリーニングし；そして、　（３）工程（２）で得られた体細胞変異の部位の塩基配列情報に基づき、細胞系統樹を作成する、ことを含む。

Description

細胞の系統解析を行う方法
　本発明は、細胞の系統解析を行うための方法、システム及びプログラムに関する。
　体細胞変異（体細胞突然変異）は、発生や加齢の過程で蓄積され、体細胞ゲノムのモザイク化を引き起こす。体細胞変異は、認知症、循環器疾患、癌などの様々な疾患（非特許文献１、２）や、加齢と関連している（非特許文献３）。さらに、体細胞変異は、中枢神経系における神経細胞集団の多様性の獲得など、正常な発生にも関与していることも知られている（非特許文献４、５）。また、ＲＮＡ配列データ解析により、体細胞変異が正常組織でクローン的に拡大していることが知られている（非特許文献６）。
　ヒトを含む哺乳動物の発生初期には、１回の細胞複製あたり２．８（９５％信頼区間、２．４－３．３）個の体細胞変異が発生する。これは、生殖細胞における突然変異（生殖細胞変異）の発生率よりもわずかに多い（非特許文献７、８）。しかし、発生後期には、様々な細胞種において、細胞分裂あたり塩基対あたり６．４×１０－１０から７．８×１０－１０までの体細胞突然変異が起こるようになる（非特許文献８、９）。これらは、生殖細胞変異の１０倍に達することもある（非特許文献１０、１１）。その結果、一卵性双生児であっても、体細胞変異によって遺伝的に多様になることがある（非特許文献１２）。
　このように、体細胞変異は、ヒトを含む様々な生物の至る場所・時期において生じる事象である（非特許文献１３）。体細胞変異は、生殖細胞変異（非特許文献１４）と対をなすものと考えることができ、その動態は進化の過程に左右される（非特許文献１５）。生殖細胞同様、体細胞は時間とともにゲノムに「傷跡」として「進化の記録」を持ち、理論的には遡及的な細胞追跡法によってその歴史を復元することができる（非特許文献１６、１７）。すなわち、体細胞変異は細胞系統の時間情報を持つ。
　例えば、特に腫瘍学の分野では、体細胞変異は癌の進化との関連で広く研究されている（非特許文献１８）。研究者は、体細胞変異の動態を様々な生物学的観点から解明するために、数理モデルを考案している（非特許文献１９、２０）。
　しかし、癌の進化の研究において、ｄｅ　ｎｏｖｏ変異を検出するためには困難な複数の問題を解消する必要がある。次世代シーケンス（ＮＧＳ）技術のエラー率は非常に高く、稀な変異を検出することは困難である。この難問を解決するためには、通常、精度を上げるために、超深度シーケンシング（非特許文献２１）が必要である。また、真の体細胞変異と生殖細胞変異を区別する必要がある（非特許文献２２）。このような、エラー率の高さ、体細胞変異を生殖細胞変異と区別して検出する必要性、などの課題により、子孫細胞における祖先の変異を推定するための「系統的」特徴（ｓｉｇｎａｔｕｒｅ）にノイズを発生させる可能性がある。
　一方、シングルセルシーケンス（ＳＣＳ）の登場により、希少細胞種、未培養微生物、体組織のモザイク化など、様々な分野での解析が可能になった（非特許文献２３）。
　体細胞変異はゲノムの変化に起因するため、体細胞ゲノム配列データを用いて、接合体参照ゲノムに対する希少変異体を検出することは、明快なアプローチである。しかし、ＳＣＳでは、１つの細胞には各ゲノムＤＮＡ分子のコピーが２つしかないため、カバー範囲の不均一性、対立遺伝子脱落（ＡＤＯ）事象、偽陽性（ＦＰ）エラー、偽陰性（ＦＮ）エラーなど、さまざまな技術的問題を引き起こすという欠点がある（非特許文献２４）。
　ｔ分布型確率的近傍埋め込み法（ｔ－ｄｉｓｔｒｉｂｕｔｅｄ　ｓｔｏｃｈａｓｔｉｃ　ｎｅｉｇｈｂｏｒ　ｅｍｂｅｄｄｉｎｇ：ｔ－ＳＮＥ）（非特許文献３２、３３）及び低次元化のためのＵＭＡＰ（ｕｎｉｆｏｒｍ　ｍａｎｉｆｏｌｄ　ａｐｐｒｏｘｉｍａｔｉｏｎ　ａｎｄ　ｐｒｏｊｅｃｔｉｏｎ　ｆｏｒ　Ｄｉｍｅｎｓｉｏｎ　Ｒｅｄｕｃｔｉｏｎ）（非特許文献３４、３５）など、機械学習アプローチによる細胞種の分類や疑似時間経過解析が開発されている。しかしながら、これらの手法では、低次元データのクラスタリングを含め、生物学的に適切な解釈や結果の再現性に欠けることがある。
Kennedy, S. R., Loeb, L. A. & Herr, A. J. Somatic mutations in aging, cancer and neurodegeneration. Mech Ageing Dev 133, 118-126, doi:10.1016/j.mad.2011.10.009 (2012). Morley, A. A. The somatic mutation theory of ageing. Mutat Res 338, 19-23 (1995). Kelly, D. P. Ageing theories unified. Nature 470, 342, doi:10.1038/nature09896 (2011). Abeliovich, A. et al. On somatic recombination in the central nervous system of transgenic mice. Science 257, 404-410 (1992). McConnell, M. J. et al. Failed clearance of aneuploid embryonic neural progenitor cells leads to excess aneuploidy in the Atm-deficient but not the Trp53-deficient adult cerebral cortex. J Neurosci 24, 8090-8096, doi:10.1523/JNEUROSCI.2263-04.2004 (2004). Yizhak, K. et al. RNA sequence analysis reveals macroscopic somatic clonal expansion across normal tissues. Science 364, eaaw0726, doi:10.1126/science.aaw0726 (2019). Ju, Y. S. et al. Somatic mutations reveal asymmetric cellular dynamics in the early human embryo. Nature 543, 714-718, doi:10.1038/nature21703 (2017). Tomasetti, C., Vogelstein, B. & Parmigiani, G. Half or more of the somatic mutations in cancers of self-renewing tissues originate prior to tumor initiation. Proceedings of the National Academy of Sciences 110, 1999, doi:10.1073/pnas.1221068110 (2013). Welch, John S. et al. The Origin and Evolution of Mutations in Acute Myeloid Leukemia. Cell 150, 264-278, doi:https://doi.org/10.1016/j.cell.2012.06.023 (2012). Lynch, M. Rate, molecular spectrum, and consequences of human mutation. Proceedings of the National Academy of Sciences 107, 961, doi:10.1073/pnas.0912629107 (2010). Van Horebeek, L., Dubois, B. & Goris, A. Somatic Variants: New Kids on the Block in Human Immunogenetics. Trends in Genetics 35, 935-947, doi:https://doi.org/10.1016/j.tig.2019.09.005 (2019). Jonsson, H. et al. Differences between germline genomes of monozygotic twins. Nature Genetics 53, 27-34, doi:10.1038/s41588-020-00755-1 (2021). Garcia-Nieto, P. E., Morrison, A. J. & Fraser, H. B. The somatic mutation landscape of the human body. Genome Biology 20, 298, doi:10.1186/s13059-019-1919-5 (2019). Milholland, B. et al. Differences between germline and somatic mutation rates in humans and mice. Nat Commun 8, 15183, doi:10.1038/ncomms15183 (2017). Rozhok, A. I. & DeGregori, J. Toward an evolutionary model of cancer: Considering the mechanisms that govern the fate of somatic mutations. Proc Natl Acad Sci U S A 112, 8914-8921, doi:10.1073/pnas.1501713112 (2015). Woodworth, M. B., Girskis, K. M. & Walsh, C. A. Building a lineage from single cells: genetic techniques for cell lineage tracking. Nature Reviews Genetics 18, 230, doi:10.1038/nrg.2016.159 (2017). Oota, S. Somatic mutations - Evolution within the individual. Methods 176, 91-98, doi:https://doi.org/10.1016/j.ymeth.2019.11.002 (2020). McGranahan, N. & Swanton, C. Clonal Heterogeneity and Tumor Evolution: Past, Present, and the Future. Cell 168, 613-628, doi:10.1016/j.cell.2017.01.018 (2017). Beerenwinkel, N., Schwarz, R. F., Gerstung, M. & Markowetz, F. Cancer Evolution: Mathematical Models and Computational Inference. Systematic Biology 64, e1-e25, doi:10.1093/sysbio/syu081 (2015). Altrock, P., Liu, L. & Michor, F. The mathematics of cancer: Integrating quantitative models. Nature Reviews Cancer 15, 730-745, doi:10.1038/nrc4029 (2015). Rheinbay, E. et al. Recurrent and functional regulatory mutations in breast cancer. Nature 547, 55-60, doi:10.1038/nature22992 (2017). Sun, J. X. et al. A computational approach to distinguish somatic vs. germline origin of genomic alterations from deep sequencing of cancer specimens without a matched normal. PLoS Comput Biol 14, e1005965, doi:10.1371/journal.pcbi.1005965 (2018). Method of the year 2013. Nat Methods 11, 1, doi:10.1038/nmeth.2801 (2014). Wang, Y. & Navin, N. E. Advances and applications of single-cell sequencing technologies. Mol Cell 58, 598-609, doi:10.1016/j.molcel.2015.05.005 (2015). Sheng, Q., Zhao, S., Li, C. I., Shyr, Y. & Guo, Y. Practicability of detecting somatic point mutation from RNA high throughput sequencing data. Genomics 107, 163-169, doi:10.1016/j.ygeno.2016.03.006 (2016). Tam, P. P. L. & Ho, J. W. K. Cellular diversity and lineage trajectory: insights from mouse single cell transcriptomes. Development 147, dev179788, doi:10.1242/dev.179788 (2020). Ji, Z. & Ji, H. TSCAN: Pseudo-time reconstruction and evaluation in single-cell RNA-seq analysis. Nucleic acids research 44, e117-e117, doi:10.1093/nar/gkw430 (2016). Hou, W. et al. A statistical framework for differential pseudotime analysis with multiple single-cell RNA-seq samples. bioRxiv : the preprint server for biology, 2021.2007.2010.451910, doi:10.1101/2021.07.10.451910 (2021). Campbell, K. R. & Yau, C. Uncovering pseudotemporal trajectories with covariates from single cell and bulk expression data. Nature communications 9, 2442-2442, doi:10.1038/s41467-018-04696-6 (2018). Felsenstein, J. The Number of Evolutionary Trees. Systematic Biology 27, 27-33, doi:10.2307/2412810 (1978). Gott, J. R., III et al. A Map of the Universe. The Astrophysical Journal 624, 463-484, doi:10.1086/428890 (2005). Hinton、G. & Roweis, S. Stochastic Neighbor Embedding. Advances in neural information processing systems 15, 833--840 (2003). Maaten, L. v. d. & Hinton, G. Visualizing Data using t-SNE. Journal of Machine Learning Research 9, 2579-2605 (2008). McInnes, L., Healy, J. & Melville, J. UMAP: Uniform Manifold Approximation and Projection for Dimension Reduction. arXiv, 1802.03426 (2020). McInnes, L., Healy, J., Saul, N. & Grosberger, L. UMAP: Uniform Manifold Approximation and Projection. Journal of Open Source Software 3, 861 (2018). Pavlicev, M. et al. Single-cell transcriptomics of the human placenta: inferring the cell communication network of the maternal-fetal interface. Genome Res 27, 349-361, doi:10.1101/gr.207597.116 (2017). Sanchez, C., et al. Grasping at molecular interactions and genetic networks in Drosophila melanogaster using FlyNets, an Internet database. Nucleic acids research 27, 89-94, doi:10.1093/nar/27.1.89 (1999). Schneider, V. A. et al. Evaluation of GRCh38 and de novo haploid genome assemblies demonstrates the enduring quality of the reference assembly. bioRxiv, 072116, doi:10.1101/072116 (2016). Church, D. M. et al. Modernizing reference genome assemblies. PLoS Biol 9, e1001091, doi:10.1371/journal.pbio.1001091 (2011). Li, H. & Durbin, R. Fast and accurate short read alignment with Burrows-Wheeler transform. Bioinformatics 25, 1754-1760, doi:10.1093/bioinformatics/btp324 (2009). Bolger, A. M., Lohse, M. & Usadel, B. Trimmomatic: a flexible trimmer for Illumina sequence data. Bioinformatics 30, 2114-2120, doi:10.1093/bioinformatics/btu170 (2014). Li, H. et al. The Sequence Alignment/Map format and SAMtools. Bioinformatics 25, 2078-2079, doi:10.1093/bioinformatics/btp352 (2009). Kimura, M. The number of heterozygous nucleotide sites maintained in a finite population due to steady flux of mutations. Genetics 61, 893-903 (1969). Tajima, F. Infinite-allele model and infinite-site model in population genetics. Journal of Genetics 75, 27, doi:10.1007/BF02931749 (1996). Cingolani, P. et al. A program for annotating and predicting the effects of single nucleotide polymorphisms, SnpEff: SNPs in the genome of Drosophila melanogaster strain w1118; iso-2; iso-3. Fly (Austin) 6, 80-92, doi:10.4161/fly.19695 (2012). Farris, J. S. Methods for Computing Wagner Trees. Systematic Biology 19, 83-92, doi:10.1093/sysbio/19.1.83 (1970). Fitch, W. M. Toward Defining the Course of Evolution: Minimum Change for a Specific Tree Topology. Systematic Zoology 20, 406-416, doi:10.2307/2412116 (1971). Tamura, K., Stecher, G. & Kumar, S. MEGA11: Molecular Evolutionary Genetics Analysis Version 11. Molecular Biology and Evolution 38, 3022-3027, doi:10.1093/molbev/msab120 (2021). Felsenstein, J. PHYLIP - Phylogeny Inference Package (Version 3.2). Cladistics 5, 164-166 (1989). Yadav, V. K., DeGregori, J. & De, S. The landscape of somatic mutations in protein coding genes in apparently benign human tissues carries signatures of relaxed purifying selection. Nucleic Acids Res 44, 2075-2084, doi:10.1093/nar/gkw086 (2016). Persi, E., Wolf, Y. I., Leiserson, M. D. M., Koonin, E. V. & Ruppin, E. Criticality in tumor evolution and clinical outcome. Proceedings of the National Academy of Sciences 115, E11101-E11110, doi:10.1073/pnas.1807256115 (2018). Yang, Z. PAML 4: phylogenetic analysis by maximum likelihood. Mol Biol Evol 24, 1586-1591, doi:10.1093/molbev/msm088 (2007). Nei, M. & Gojobori, T. Simple methods for estimating the numbers of synonymous and nonsynonymous nucleotide substitutions. Mol Biol Evol 3, 418-426, doi:10.1093/oxfordjournals.molbev.a040410 (1986). R Core Team. R: A Language and Environment for Statistical Computing (Vienna, Austria, 2016). Blondel, V. D., Guillaume, J.-L., Lambiotte, R. & Lefebvre, E. Fast unfolding of communities in large networks. Journal of Statistical Mechanics: Theory and Experiment 2008, P10008, doi:10.1088/1742-5468/2008/10/p10008 (2008). Trapnell, C. et al. The dynamics and regulators of cell fate decisions are revealed by pseudotemporal ordering of single cells. Nature Biotechnology 32, 381-386, doi:10.1038/nbt.2859 (2014). Qiu, X. et al. Reversed graph embedding resolves complex single-cell trajectories. Nature Methods 14, 979-982, doi:10.1038/nmeth.4402 (2017). Cao, J. et al. The single-cell transcriptional landscape of mammalian organogenesis. Nature 566, 496-502, doi:10.1038/s41586-019-0969-x (2019). Mathematica (Wolfram Research, Inc., Champaign, Illinois, 2020). Phylogenetics for Mathematica. Version 6.5 (Indiana University: Bloomington, Indiana, Department of Earth and Atmospheric Sciences, 2019). Zachar, I. IstvanZachar/Phylogenetics, <https://github.com/IstvanZachar/Phylogenetics/releases/tag/1.1.0> (2017). Archie, J. et al. The Newick tree format, <https://evolution.genetics.washington.edu/phylip/newicktree.html> (1986). Kumar, A. An overview of nested genes in eukaryotic genomes. Eukaryot Cell 8, 1321-1329, doi:10.1128/EC.00143-09 (2009). Long, C. A. Sokal, Robert R., and Peter H. A. Sneath. Principles of Numerical Taxonomy. W. H. Freeman and Co., San Francisco and London. Pp. xvi + 359, illus. 1963. Price $8.50. Journal of Mammalogy 46, 111-112, doi:10.2307/1377831 (1965). Knofler, M. et al. Human placenta and trophoblast development: key molecular mechanisms and model systems. Cellular and Molecular Life Sciences 76, 3479-3496, doi:10.1007/s00018-019-03104-6 (2019). Yang, Z., Wong, W. S. W. & Nielsen, R. Bayes Empirical Bayes Inference of Amino Acid Sites Under Positive Selection. Molecular Biology and Evolution 22, 1107-1118, doi:10.1093/molbev/msi097 (2005). Zhang, T., Periz, G., Lu, Y. N. & Wang, J. USP7 regulates ALS-associated proteotoxicity and quality control through the NEDD4L-SMAD pathway. Proc Natl Acad Sci U S A 117, 28114-28125, doi:10.1073/pnas.2014349117 (2020). Bittoni, A. et al. Retrospective Cohort Study of Caveolin-1 Expression as Prognostic Factor in Unresectable Locally Advanced or Metastatic Pancreatic Cancer Patients. Curr Oncol 28, 3525-3536, doi:10.3390/curroncol28050303 (2021). Kobak, D. & Berens, P. The art of using t-SNE for single-cell transcriptomics. Nature Communications 10, 5416, doi:10.1038/s41467-019-13056-x (2019). Ortega, M. A. et al. Using single-cell multiple omics approaches to resolve tumor heterogeneity. Clinical and Translational Medicine 6, 46, doi:10.1186/s40169-017-0177-y (2017). Araten, D. J. et al. A Quantitative Measurement of the Human Somatic Mutation Rate. Cancer Research 65, 8111-8117, doi:10.1158/0008-5472.Can-04-1198 (2005). Dou, Y., Gold, H. D., Luquette, L. J. & Park, P. J. Detecting Somatic Mutations in Normal Cells. Trends Genet 34, 545-557, doi:10.1016/j.tig.2018.04.003 (2018). Rhee, J.-K., Lee, S., Park, W.-Y., Kim, Y.-H. & Kim, T.-M. Allelic imbalance of somatic mutations in cancer genomes and transcriptomes. Scientific Reports 7, 1653, doi:10.1038/s41598-017-01966-z (2017). Ju, Y. S. et al. Extensive genomic and transcriptional diversity identified through massively parallel DNA and RNA sequencing of eighteen Korean individuals. Nature Genetics 43, 745-752, doi:10.1038/ng.872 (2011). Neums, L. et al. VaDiR: an integrated approach to Variant Detection in RNA. Gigascience 7, 1-13, doi:10.1093/gigascience/gix122 (2018). Browning, S. R. & Browning, B. L. Haplotype phasing: existing methods and new developments. Nature reviews. Genetics 12, 703-714, doi:10.1038/nrg3054 (2011). Wasik, K. et al. Comparing low-pass sequencing and genotyping for trait mapping in pharmacogenetics. BMC Genomics 22, 197, doi:10.1186/s12864-021-07508-2 (2021). Cock, P. J. A., Fields, C. J., Goto, N., Heuer, M. L. & Rice, P. M. The Sanger FASTQ file format for sequences with quality scores, and the Solexa/Illumina FASTQ variants. Nucleic acids research 38, 1767-1771, doi:10.1093/nar/gkp1137 (2010). Danecek, P. et al. The variant call format and VCFtools. Bioinformatics (Oxford, England) 27, 2156-2158, doi:10.1093/bioinformatics/btr330 (2011).
　本発明者らは、上記問題解決のため鋭意研究に努めた結果、トランスクリプトームの配列データを参照ゲノム配列のデータと比較することにより同定される、体細胞変異に着目することにより、細胞の系統解析を行うことを見出し、本発明を想到した。
　本発明は、単一細胞のトランスクリプトームが持つ系統学的なシグネチャーに着目するものである。例えば、６０個の体細胞は約１０９６の可能な系統を持ち、これは宇宙に存在する原子の数（非特許文献３１）より多い。本発明の方法は、ＲＮＡ配列データに保持されている粗い系統学的シグナルを利用して、真の系統をカバーするウィンドウを絞り込むものである。シングルセルシーケンス（ＳＣＳ）と異なり、ＲＮＡ配列のソースでは、細胞内のＲＮＡ分子の数百、数千のコピーを利用することができる。ＲＮＡ編集、配列決定エラー（例えば逆転写やＰＣＲの際に生じるランダムエラー）、潜在的なサンプリングエラーなど、生物学的および技術的要因による偽陽性を考慮する必要がある（非特許文献１３、２５）。しかし、単一細胞のトランスクリプトームデータは圧倒的に豊富で、ＲＮＡシーケンスデータのこれらのデメリットを補うものである（非特許文献２６）。
　限定されるわけではないが、本発明は、以下の態様を含む。
［１］
　細胞の系統解析を行う方法であって、
　（１）同一個体由来の複数の単一細胞について、各単一細胞のトランクリプトームのＲＮＡ配列データを、対応するゲノム配列データと比較し、ＲＮＡ配列とゲノム配列が異なる部位をスクリーニングし；
　（２）工程（１）のスクリーニングで特定された部位のうち、３つ若しくはそれ以上の１塩基多型が検出される部位を、体細胞変異の部位としてスクリーニングし；そして、
　（３）工程（２）で得られた体細胞変異の部位の塩基配列情報に基づき、細胞系統樹を作成する、
ことを含む、前記方法。
［２］
　（４）前記細胞系統樹から、各単一細胞の細胞の種類を推定する、ことをさらに含む、［１］に記載の方法。
［３］
　各単一細胞の遺伝子発現プロファイルから推定される単一細胞の種類の情報を、前記細胞系統樹の情報と比較する、ことをさらに含む、［１］又は［２］の方法。
［４］
　１塩基多型が、１塩基置換である、［１］－［３］のいずれか１項に記載の方法。
［５］
　工程（１）において、全細胞のうち５０％以上の細胞において、トランクリプトームのＲＮＡ配列データが得られない部位については、変異は生じていないと推定し、スクリーニングから除外する、［１］－［４］のいずれか１項に記載の方法。
［６］
　工程（２）でスクリーニングする体細胞変異が、頻度の少ない体細胞変異である、［１］－［５］のいずれか１項に記載の方法。
［７］
　疾患又は症状と関連する体細胞変異を検出するための、［１］－［６］のいずれか１項に記載の方法。
［８］
　疾患又は症状が、癌、認知症、循環器疾患、加齢、自己免疫疾患、神経変性疾患、及び精神疾患からなる群から選択される、［７］に記載の方法。
［９］
　［１］－［８］のいずれか１項に記載の方法を実行するためのプログラム。
［１０］
　［１］－［８］のいずれか１項に記載の方法を実行するためのプログラムを記録した記録媒体。
［１１］
　プロセッサと、前記プロセッサで実行されたときに、［１］－［８］のいずれか１項に記載の方法を実行するプログラムを記録したメモリーを有するシステム。
　本発明の方法は、単一細胞のトランスクリプトームデータから検出される潜在的な体細胞変異を利用して、細胞系列の軌跡を検出する新しい枠組みを提供する。不均質に分化した細胞を１回のタイムスナップショットで解析する擬似時間経過解析とは異なり、本発明では体細胞変異を用いて祖先細胞まで細胞系譜を遡ることができる。
図１は、単一細胞の体細胞変異パターン分析のための経路を示す。ＦＡＳＴＱ（非特許文献７８）、ＢＡＭ（非特許文献４）及びＶＣＦ（非特許文献７９）ファイルを各単一細胞について作成した。ＦＡＳＴＱ、ＢＡＭ、及びＶＣＦは、遺伝情報を記述するためのデータ形式である。 図２は、多変異部位を有する変異体による体細胞変異の検出を示す。（ａ）生殖細胞変異及び体細胞変異を有しない祖先ヘテロ接合部位由来の見掛け上の変異部位。（ｂ）体細胞変異を有する祖先ヘテロ接合部位由来の見掛け上の変異部位。この細胞の変異の状況を（ａ）と区別することはできない。（ｃ）体細胞変異を有する祖先ヘテロ接合部位由来の変異部位。この細胞の変異の状況は、第３のタイプのヌクレオチドＧにより、（ａ）、（ｂ）と区別可能である。（ｄ）生殖細胞変異を有する祖先ヘテロ接合部位由来の見掛け上の変異部位。この細胞の変異の状況を（ｈ）と区別することはできない。（ｅ）生殖細胞変異、次いで、体細胞変異を有する祖先ヘテロ接合部位由来の変異部位。この細胞の変異の状況は、第３のタイプのヌクレオチドＴにより、（ｄ）、（ｈ）と区別可能である。（ｆ）ホモ接合部位由来のホモ部位。（ｇ）体細胞変異を有する祖先ヘテロ接合部位由来の変異部位。この細胞の変異の状況を（ａ）、（ｂ）と区別することはできない。（ｈ）生殖細胞変異、次いで、体細胞変異を有する祖先ヘテロ接合部位由来の見掛け上の変異部位。この細胞の変異の状況を（ｄ）と区別することはできない。図２では、簡便のため参照ゲノムはハプロタイプフェーズとした。実線は、ハプロタイプフェーズの参照部位を示す。点線は、ホモ接合体参照部位を示す。本図では、例示的に、祖先細胞由来の典型的な２細胞を表現している。 図３は、各ノードにおける単一系統性（ｍｏｎｏｐｈｙｌｅｔｉｃｉｔｙ）の程度（ＤｏＭ）の評価を、サブツリーにおけるクラスターの割合で示したものである。遺伝子発現プロファイル（クラスターの表）をツリートポロジーにマップする。ＤｏＭと遺伝子発現プロファイルが互いに完全に一致する場合、クラスターの観点で単一系統である、サブツリーのセットが観察され得る。ＤｏＭと遺伝子発現プロファイルが完全には一致しない場合、クラスターの観点で、側系統のサブツリー又は多系統のサブツリーが観察されるであろう。このようにしてサブツリーの総数に対する、単一系統のサブツリーのサイズの一致性を評価することができる。Ｍａｔｈｅｍａｔｉｃａ（非特許文献５９）のコード、ＡｓｓｉｇｎＣｌｕｓｔｅｒ２Ｃｅｌｌにより、クラスターに割り当てられた細胞の数の割合を算出し、それによって、サブツリーのＤｏＭの単一系統性を評価することができる。Ａ：祖先幹細胞；Ｂ、Ｃ及びＤ：幹細胞由来細胞；Ｅ、Ｆ、Ｇ、及びＩ：観察された分化した単一細胞。 図４は、参照ゲノム（ＧＲＣｈ３８）上にマップされたトランスクリプトームのカバー割合を示す。カバー範囲については、各細胞によって大きな相違があった。横軸：参照ゲノムにおけるカバー割合；縦軸：２名の妊婦からサンプリングされたバッチ１及びバッチ２の各単一細胞の名称。 図５は、遺伝子発現プロファイルの細胞系統樹へのマッピングである（ＳＲＰ０９０９４４データ；バッチ１、５４細胞）。Ｐａｖｌｉｃｅｖらの研究（非特許文献３６）に基づき、遺伝子発現プロファイルと３００のマーカー遺伝子に関する主成分分析（ｐｒｉｎｃｉｐａｌ　ｃｏｍｐｏｎｅｎｔ　ａｎａｌｙｓｉｓ（ＰＣＡ）により、単一細胞データを、細胞栄養芽層（ＣＹＴ）１、ＣＹＴ２、ＣＹＴ３、絨毛外栄養膜（ＥＶＴ）、母体脱落膜細胞（ＤＣ）の５つのクラスター（細胞の種類）に分類し、ツリートポロジーにマッピングした。（ａ）以前のＰａｖｌｉｃｅｖらの研究（非特許文献３６）では母体脱落膜細胞（ＤＣ）として報告されているが、胎盤の細胞；（ｂ）注釈されていない胎盤細胞（非特許文献３６）；（ｃ）自己再生段階の推定幹細胞。スケールバー：変異の期待値。円グラフ中の各数値は、遺伝子発現プロファイルに基づいて指定されたクラスターを示す。 図６は、遺伝子発現プロファイルの細胞系統樹へのマッピングである（ＳＲＰ０９０９４４データ；バッチ２、２３細胞）。（ａ）自己再生段階の推定幹細胞。その他の表記は、図５と同様である。 図７は、ｔ－ＳＮＥにより再分析された遺伝子発現プロファイルの細胞系統樹（ＳＲＰ０９０９４４データ；バッチ１、５４細胞）へのマッピングである。 図８は、ＳＲＰ０９０９４４データ；バッチ１の擬似時間経過解析の結果を示す。等高線は、補間次数３で、ＭａｔｈｅｍａｔｉｃａのＬｉｓｔＣｏｎｔｏｕｒＰｌｏｔ関数（非特許文献５９）によって作成された。図は、擬似時間を色で示している。Ｍａｔｈｅｍａｔｉｃａアルゴリズムのため、本図では全ての細胞が標識されているわけではない。 図９は、Ｎｏｔｃｈ／Ｗｎｔシグナルの機能を細胞栄養芽層（ｃｙｔｏｔｒｏｐｈｏｂｌａｓｔ）と絨毛外栄養膜（ｅｘｔｒａｖｉｌｌｏｕｓ　ｔｒｏｐｈｏｂｌａｓｔ）の分化（非特許文献６５）に結び付けるモデルシステムに関する。ｖＣＴＢ：絨毛細胞栄養芽層；ＣＣＴ：細胞柱栄養膜；ＥＶＴ：絨毛外栄養膜；Ｎ１ＩＣＤ：Ｎｏｔｃｈ１細胞間ドメイン；ＴＥＡＤ４：転写エンハンサー因子ＴＥＦ－３遺伝子；ｐ６３：腫瘍タンパク質ｐ６３；ＴＣＦ１：転写因子Ｔ細胞因子１；ＴＣＦ４：転写因子Ｔ細胞因子４；ＩＲＦ６：インターフェロン制御因子６；Ｎｏｔｃｈ２：Ｎｏｔｃｈ受容体２。
　非限定的に本発明は、以下の態様を含む。本明細書において他に断りがない限り、本明細書で使用される技術及び科学用語は、当業者が通常理解している意味と同じ意味を有する。本明細書に開示された物質、材料及び例は単なる例示であり、制限することを意図していない。本明細書において「一態様において」と言及する場合は、その態様に限定されない、即ち、非限定的であることを意味する。
　１．細胞の系統解析を行う方法
　一態様において、本発明は、細胞の系統解析を行う方法に関する。
　前記方法は、
　（１）同一個体由来の複数の単一細胞について、各単一細胞のトランクリプトームのＲＮＡ配列データを、対応するゲノム配列データと比較し、ＲＮＡ配列とゲノム配列が異なる部位をスクリーニングし；
　（２）工程（１）のスクリーニングで特定された部位のうち、３つ若しくはそれ以上の１塩基多型が検出される部位を、体細胞変異の部位としてスクリーニングし；そして、
　（３）工程（２）で得られた体細胞変異の部位の塩基配列情報に基づき、細胞系統樹を作成する、
ことを含む。
　前記発明の方法では、同一個体由来の複数の単一細胞を利用する。個体の種類は特に限定されない。非限定的に、ヒト、非ヒト霊長類（サル、ゴリラ、チンパンジー等）、マウス、ラット等を含む。一態様において、同一個体の、同一組織又は近接する組織由来の複数の単一細胞を利用する。近傍する組織とは、非限定的に、生体内において物理的に近い組織（例えば、小腸と大腸、胎盤と胎児等）、機能が関連する組織（例えば、消化器系、呼吸器系、中枢神経系等）を意味する。
　各単一細胞のトランクリプトームのＲＮＡ配列データの取得先は特に限定されない。例えば、ＤＤＢＪ、ＥＭＢＬ－Ｂａｎｋ／ＥＢＩ、ＧｅｎＢａｎｋ／ＮＣＢＩ等のデータバンクで公開されているデータ、非特許文献３６等文献に公開されているデータを利用することができる。
　トランクリプトームのＲＮＡ配列データと比較する、対応するゲノム配列データの取得先は特に限定されない。例えば、ＤＤＢＪ、ＥＭＢＬ－Ｂａｎｋ／ＥＢＩ、ＧｅｎＢａｎｋ／ＮＣＢＩ等のデータバンクで公開されているデータを利用することができる。
　各単一細胞のトランクリプトームのＲＮＡ配列データを、対応するゲノム配列データと比較し、ＲＮＡ配列とゲノム配列が異なる部位をスクリーニングする工程は、配列を比較、マッピングするための公知のソフトウェアを使用して行うことができる。
　例えば、マッピングは、Ｂｕｒｒｏｗｓ－Ｗｈｅｅｌｅｒ　Ａｌｉｇｎｅｒ（ＢＷＡ）を用いることができる。ＢＷＡは、ヒトゲノムのような長大な参照ゲノムに対して、低域配列をマッピングするためのソフトウェアパッケージであり、ＢＷＡ－バックトラック、ＢＷＡ－ＳＷ及びＢＷＡ－ＭＥＭの２種類のアルゴリズムを含む。トランクリプトームのＲＮＡ配列データの対応するゲノム配列データへのマッピングには、その他に、ＢａｒｒａＣＵＤＡ、ＳＴＡＲなどのソフトウェアを使用してもよい。
　一態様において、解析のために、ＲＮＡ配列にアダプター配列を付加してもよい。その場合、ＲＮＡ配列データのゲノム配列データへのマッピングの前に、アダプター配列を削除することが好ましい。
　ＲＮＡ配列とゲノム配列が異なる部位のスクリーニングは、多型の検出やアライメントの表示のための公知のソフトウェアを使用して行うことができる。例えば、Ｓａｍｔｏｏｌｓを使用してもよい。Ｓａｍｔｏｏｌｓは、ＤＮＡ配列決定のショートリード配列アラインメント後の処理に用いられるソフトウェアである。Ｓａｍｔｏｏｌｓの主な機能としては、多型の検出やアライメントの表示、インデクシング、データ抽出、ファイル形式の変換などがある。ＲＮＡ配列とゲノム配列が異なる部位のスクリーニングには、その他に、ｅｌＰｒｅｐ、Ｐｉｃａｒｄなどのソフトウェアを使用してもよい。
　工程（２）において、工程（１）のスクリーニングで特定された部位のうち、３つ若しくはそれ以上の１塩基多型が検出される部位を、体細胞変異の部位としてスクリーニングする。
　前記発明では、体細胞変異の部位をスクリーニングする。例えば、Ｔｏｍａｓｅｔｔｉらによると（非特許文献８）、慢性リンパ性白血病（ＣＬＬ）の前駆体である正常リンパ球におけるｉｎ　ｖｉｖｏ組織特異的体細胞変異確率は、１細胞分裂あたり１塩基あたり７．６×１０－１０±１．１×１０－１０（ＳＥ）であった。したがって、確率的な体細胞変化は、反復領域を除いたヒトゲノムにおいて、１細胞分裂あたり平均３個の突然変異が起こっていることになる（非特許文献７０、７１）。このような体細胞変異の多さが、正常組織であっても過去に遡る細胞系譜の作成を可能にする。
　「３つの１塩基多型」とは、３種類の塩基多型という意味である。仮にある遺伝子座（ここではゲノム上のある座標）がホモ接合型であった場合、標準ゲノム配列との比較によりその遺伝子座で変異が生じたかどうかを判断することができる。しかし、ある遺伝子座がヘテロ接合型であった場合は、それが変異の結果なのかもともとヘテロ接合型であったのかを即座に判断することができない。即ち、変異が２種類の場合は、その遺伝子座で変異が起きているかどうかを判断することが困難である。しかし、変異が３種類の場合は、少なくとも１回の突然変異を仮定しなければ、この観察結果を説明することができない。本発明では、この論理に基づき、「３つ若しくはそれ以上の１塩基多型が検出される部位」を変異検出の対象とする。
　「３つ若しくはそれ以上の１塩基多型が検出される部位」の多型の数の上限は特に限定されない。解析の便宜上は、多型の数はあまり多くなりすぎない方が好ましい。非限定的に、多型の数は、１５以下、１２以下、１０以下、８以下、６以下、５以下、４以下である。一態様において、多型の数は３である。
　本発明で使用する「１塩基多型」は、１塩基レベルで生じる他の塩基（群）への変異一般を意味し、塩基（群）の置換、挿入（付加）、欠失のいずれでもよい。「１塩基レベルで生じる変異」とは、一塩基の変異（点変異）の他に、複数塩基単位で生じる変異も含む。一態様において、一塩基多型は、一塩基の変異（点変異）である。一態様において、一塩基多型は、一塩基の挿入及び／又は置換である。一態様において、一塩基多型は一塩基の置換である。
　一態様において、工程（１）において、全細胞のうち５０％以上の細胞において、トランクリプトームのＲＮＡ配列データが得られない部位については、変異は生じていないと推定し、スクリーニングから除外する。前記方法では、トランクリプトームのＲＮＡ配列データを利用する。静的なゲノムと異なり、細胞毎に、また、時間経過によりトランスクリプト（転写産物）の量は変動する。よって、細胞によっては特定の転写産物のＲＮＡが検出されずＲＮＡ配列データが得られない場合（部位）が生じうる。本発明の一態様において、全細胞のうち５０％以上の細胞において、トランクリプトームのＲＮＡ配列データが得られた場合、即ち、５０％以上の細胞にマップされたトランスクリプト配列を変異検出に利用する。この場合、工程（１）において、全細胞のうち５０％以上の細胞において、トランクリプトームのＲＮＡ配列データが得られない部位については、変異は生じていないと推定し、スクリーニングから除外する。「変異は生じていないと推定し」とは、対応するゲノム配列を代わりに用いることを意味する。「全細胞のうち５０％以上の細胞において」の「５０％以上」は、非限定的に、５５％以上、６０％以上、６５％以上、７０％以上、７５％以上、８０％以上であってもよい。
　トランスクリプト（転写産物）のカバー範囲が低いため、一部の単一細胞では最初のスクリーニングで潜在的な変異部位がフィルタリングされる可能性がある。本発明のスクリーニングにおいて、潜在的な変異体を除外するための閾値を選択してもよい。本明細書の実施例では、解決策の一例として８０％の閾値を使用した。具体的には、試験した単一細胞のうち８０％以上の単一細胞において観察される多変異部位をスクリーニングした。非限定的に、より高い閾値、例えば、９０％以上、あるいはより低い閾値、例えば、６０％以上、７０％以上を採用してもよい。
　前記方法の一態様において、工程（２）でスクリーニングする体細胞変異が、頻度の少ない体細胞変異である。
　頻度の少ない体細胞変異（マイナーアリル）とは、細胞集団において観察されたトランスクリプトにおいて、頻度の低い方の体細胞変異を意味する。「アリル」とは元来対立遺伝子の意味であるが、本明細書においては、一塩基変異に意味を拡張して用いる場合がある。アリルという概念は本来ゲノム配列に帰属するものであるが、トランスクリプト配列によるアリル頻度とゲノム配列によるアリル頻度の間にある程度相関があると仮定し、ゲノム配列において頻度の低い方のアリル、即ち、マイナーアリルを、スクリーニングする体細胞変異の候補とする。
　「トランスクリプト配列によるアリル頻度とゲノム配列によるアリル頻度の間にある程度相関がある」との仮定は、ゲノムとトランスクリプトームのリードの対立遺伝子頻度が相関しているという経験的な予想に基づいている（非特許文献７４，７５）。しかし、誤差を少なくするためには、より定量的なアプローチ方法をさらに適用することが好ましい。例えば、統計的なハプロタイピングの位相合わせ（非特許文献７６）は、推定値を向上させるために有効であると考えられる。
　対立遺伝子において見出される変異の頻度は、例えば、参照ゲノムとリード配列の比較により推定することができる。「頻度の少ない」の意義は、非限定的に、例えば、対立遺伝子において見出される全変異について、頻度が低い方から５０％以内、好ましくは、４０％以内、３０％以内、２０％以内、１０％以内、５％以内であることを意味する。一態様において、「頻度の少ない」の意義は、対立遺伝子において見出される全変異について、頻度が低い方から１０番以内、８番以内、５番以内、３番以内であることを意味する。一態様において、頻度の少ない体細胞変異（マイナーアリル）とは、対立遺伝子において見出される全変異について、最も頻度が低い体細胞変異を意味する。
　頻度の少ない変異（マイナーアリル）をスクリーニングすることの意味は、ある細胞において新しく生じた体細胞変異は細胞集団に拡散する途上にあるので、その観察頻度は相対的に低いはずだという仮定に基づく。前記方法において、このような頻度の低い変異を選択することにより、トランスクリプトのアリル頻度とＤＮＡ配列のアリル頻度のある程度の相関を踏まえ、複数のアリルから新しく生じた体細胞変異を効率よく選ぶことができる。もって、細胞の系統解析をより有効に行うことが可能となる。
　工程（３）において、工程（２）で得られた体細胞変異の部位の塩基配列情報に基づき、細胞系統樹を作成する。
　細胞系統樹の作成は、配列情報から系統樹を作成するための、公知のソフトウェアを使用することができる、例えば、遺伝子配列情報から系統塾を作成する（過去を推定する）パッケージソフトウェアのひとつである、ＭＥＧＡ－Ｘを使用することができる。その他にも、Ｐｈｙｌｉｐなどのソフトウェアを使用してもよい。
　非限定的に、細胞系統樹の作成において、最大節約法（ｍａｘｉｍｕｍ　ｐａｒｓｉｍｏｎｙ　ｍｅｔｈｏｄ）（非特許文献４６、４７）を適用してもよい。最大節約法とは、変異数が最小になるような変異パターンを見つけることで、系統樹を作成する方法の一種である。同様の解析は、その他の方法、例えば、距離行列法、最尤法、ベイズ法等を使用して行うことも可能である。
　前記方法において、前記細胞系統樹から、各単一細胞の細胞の種類を推定する、ことをさらに含んでもよい。
　細胞系統樹からの各細胞の細胞の種類の推定（特定）は、公知の方法を用いて行うことができる。例えば、主成分分析、ｔ－ＳＮＥ、ＵＭＡＰなどの公知のツールを用いてもよい。具体的には、単一細胞の種類と遺伝子の種類からなる遺伝子発現マトリックスを線形もしくは非線形次元削減し、細胞のクラスタリングを行なったり、既知のマーカー遺伝子の情報を用いてもよい。
　ｓｎｐＥｆｆ（非特許文献４５）とは、予測された変異にアノテーションやその影響に関する情報を付与するツールである。各変異の数や位置情報・タイプ（ＳＮＰ、挿入、欠失など）や影響度（アミノ酸への置換）の情報を得ることができる。変異部位がゲノムのどの箇所に落ちているかによって、ゲノムに割り当てられた注釈（ａｎｎｏｔａｔｉｏｎ）をもとに変異部位の注釈を行うことができる。しなしながら、前記方法で使用するデータはトランスクリプト配列であることもあり、多少の不確実性が存在する場合がある。ＳｎｐＥｆｆソフトウェアはそのあたりの事情も考慮し、適切と思われる注釈を採用することができる。さらに、これらの注釈によって各変異の生物学的・進化学的解釈が可能になる。
　一態様において、前記方法は、各単一細胞の遺伝子発現プロファイルから推定される単一細胞の種類の情報を、前記細胞系統樹の情報と比較する、ことをさらに含んでもよい。
　各単一細胞の遺伝子発現プロファイルの取得先は特に限定されない。例えばＥｘｐｒｅｓｓｉｏｎ　Ａｔｌａｓ、Ｇｅｎｅ　Ｅｘｐｒｅｓｓｉｏｎ　Ｏｍｎｉｂｕｓ等のデータバンクで公開されているデータ、非特許文献３６等文献に公開されているデータを利用することができる。
　一態様において、前記方法により、遺伝子発現プロファイルでは、特定できなかった単一細胞の種類を特定することができる。一態様において、遺伝子発現プロファイルでは、誤って特定されていた単一細胞の種類を、再検討したり、正しいものに訂正したりすることができる。
　一態様において、前記方法は、疾患又は症状と関連する体細胞変異を検出するために使用することができる。疾患又は症状の種類は特に限定されない。一態様において、疾患又は症状は、遺伝子変異に起因する疾患又は症状である。非限定的に、疾患又は症状は、癌、認知症、循環器疾患、加齢、自己免疫疾患、神経変性疾患、及び精神疾患からなる群から選択される。
　本発明の意義
　本発明の理解のため、非限定的に、本発明の意義について説明する。本発明により、カバーの単一細胞トランスクリプトームから検出された体細胞変異を用いて、細胞系統樹を推定する新しいフレームワークが提供された。その際、個々の変異ではなく、単一細胞トランスクリプトームが持つ系統的な特徴に注目した。その結果、既知の生物学的知識と一致する細胞系統樹を再構築することが可能であることを示した。本発明の意義は、単一細胞の総体的な系統樹を推測し、遺伝子発現プロファイルを解釈するための裏付けとなる情報を提供することにある。遺伝子発現解析では、次元削減法がしばしば用いられる。しかし、単一細胞レベルで得られた結果から、生物学的に適切な知見をどのように抽出するかは、未解決の問題である（非特許文献６９）。このフレームワークは、単一細胞のモザイク性という異なる角度から遺伝子発現プロファイルを解釈することを可能にする。
　擬似時間経過解析では、単一の時間スナップショットで異種分化した細胞を使用するが、本発明のアプローチであるリアルタイムコースでは、体細胞変異を利用して、子孫細胞まで細胞系譜を遡ることができる。即ち、本発明の方法は、細胞の時間経過を推測することが可能である。「時間経過の推測」とは、例えば、幹細胞の自己複製段階を含む、観察不可能な過去のイベントなどについての推測である（図５ｃおよび図６ａ）。実際、本発明のリアルタイムコースは、擬似時間経過解析の結果とほぼ一致する結果を得ることができるが、全く異なる種類のデータ、すなわちマッピングされた読み取り配列とその深さを用いることで、より細粒度な方法で結果を得ることが可能である。
　２．細胞の系統解析を行うためのプログラム、システム等
　本発明の細胞の系統解析を行う方法を、効率的に実行するためには、当該方法を計算機によって自動的に実行するプログラムを作成するのが好ましい。また、このプログラムを計算機に読み込ませるために記録媒体に記録してもよい。さらには、計算機のプロセッサで実行できるように、プログラムをメモリーに記録した専用のシステムとしてもよい。
　一態様において、本発明は、本発明の方法を実行するためのプログラムに関する。
　一態様において、本発明は、細胞の系統解析を行う方法を実行させるためのプログラムであって、前記方法は、
　（１）同一個体由来の複数の単一細胞について、各単一細胞のトランクリプトームのＲＮＡ配列データを、対応するゲノム配列データと比較し、ＲＮＡ配列とゲノム配列が異なる部位をスクリーニングする工程；
　（２）（１）のスクリーニングで特定された部位のうち、３つ若しくはそれ以上の１塩基多型が検出される部位を、体細胞変異の部位としてスクリーニングする工程；及び、
　（３）（１）で得られた体細胞変異の部位の塩基配列情報に基づき、細胞系統樹を作成する工程
を含む、前記プログラム、に関する。
　一態様において、本発明は、本発明の方法を実行するためのプログラムを記録した記録媒体に関する。
　一態様において、本発明は、プロセッサと、前記プロセッサで実行されたときに、本発明の方法に記載の方法を実行するプログラムを記録したメモリーを有する、システムに関する。
　一態様において、本発明は、細胞の系統解析を行うためのシステムであって、
　（１）同一個体由来の複数の単一細胞について、各単一細胞のトランクリプトームのＲＮＡ配列データを、対応するゲノム配列データと比較し、ＲＮＡ配列とゲノム配列が異なる部位をスクリーニングする手段；
　（２）（１）のスクリーニングで特定された部位のうち、３つ若しくはそれ以上の１塩基多型が検出される部位を、体細胞変異の部位としてスクリーニングする手段；及び、
　（３）（１）で得られた体細胞変異の部位の塩基配列情報に基づき、細胞系統樹を作成する手段
を含む、前記システム、に関する。
　「細胞の系統解析」、「（１）同一個体由来の複数の単一細胞について、各単一細胞のトランクリプトームのＲＮＡ配列データを、対応するゲノム配列データと比較し、ＲＮＡ配列とゲノム配列が異なる部位をスクリーニングする」、「（２）（１）のスクリーニングで特定された部位のうち、３つ若しくはそれ以上の１塩基多型が検出される部位を、体細胞変異の部位としてスクリーニングする」、「（３）（１）で得られた体細胞変異の部位の塩基配列情報に基づき、細胞系統樹を作成する」の意義、態様などについては、「１．細胞の系統解析を行う方法」の項目に記載した通りである。「１．細胞の系統解析を行う方法」に記載した事項はすべて、システム、プログラム、記録媒体にも適用される。
　以下、実施例に基づいて本発明を詳細に説明するが、本発明はこれらの実施例に限定されるものではない。当業者は本明細書の記載に基づいて容易に本発明に修飾・変更を加えることができ、それらは本発明の技術的範囲に含まれる。
　実施例１　トランスクリプトーム配列データのマッピングと変異の検出
　（１－１）トランスクリプトームデータ
　トランスクリプトームデータとして、正常（胎盤）組織から得られた２種類の公開トランスクリプトームデータ：ＳＲＰ０９０９４４バッチ１（５４細胞）およびバッチ２（３３細胞）を使用した（非特許文献３６）。
　非特許文献３６（Ｐａｖｌｉｃｅｖら）では、胎盤のデータを、母親と胎児という２つの半同種異系（ｓｅｍｉａｌｌｏｇｅｎｉｃ）個体間の細胞コミュニケーションネットワークという観点から解析した。当該文献では、細胞の種類を超えた受容体－リガンドペアの遺伝子発現から、細胞間インタラクトーム（非特許文献３７）を推定した。その結果、Ｇタンパク質共役型受容体の細胞型特異的な発現を見出し、リガンド－受容体プロファイルが細胞型識別のための信頼できるツールになる可能性が示唆された。データはＤＤＢＪ　Ｓｅｑｕｅｎｃｅ　Ｒｅａｄ　Ａｒｃｈｉｖｅ　（ＤＲＡ）にＳＲＳ１７３２２６６（ＳＲＸ２２２５２６９）－ＳＲＳ１７３２３１９（ＳＲＸ２２２５３２８）として登録されている。
　Ｐａｖｌｉｃｅｖらは、遺伝子発現プロファイルと３００のマーカー遺伝子に関する主成分分析（ｐｒｉｎｃｉｐａｌ　ｃｏｍｐｏｎｅｎｔ　ａｎａｌｙｓｉｓ（ＰＣＡ）により、単一細胞データを、細胞栄養芽層（ＣＹＴ）１、ＣＹＴ２、ＣＹＴ３、絨毛外栄養芽層（ＥＶＴ）、母体脱落膜細胞（ＤＣ）の５クラスター（細胞の種類）に分類している。ＤＣについては、Ｐａｖｌｉｃｅｖらは、これらの細胞は、ＣＤ１９、ＣＤ２０９、ＣＤ１６３の発現を欠き、ＩＴＧＡＸ＋／ＣＤ１４＋／ＣＤ４＋／ＣＤ８３＋／ＣＤ８６＋のマーカーを組み合わせて発現していることを報告し、子宮樹状細胞であると判断している。しかし、ＤＣの特徴的な遺伝子のいくつか、例えば、ＣＬＥＣ４Ｃ、ＴＨＢＤ、ＣＤ１Ｃ、ＣＤ８０、ＩＬ１０、ＩＬ１２Ｂなど。の発現は見られなかった。
　（１－２）トランスクリプトーム配列データのマッピング
　本実施例のデータ解析パイプラインの全体像を図１に示す。ｔｒｉｍｏｍａｔｉｃ（非特許文献４１）でアダプター配列を削除した後に、単一細胞トランスクリプトーム配列データを、Ｂｕｒｒｏｗｓ－Ｗｈｅｅｌｅｒ　Ａｌｉｇｎｅｒ（ＢＷＡ）４０を用いてヒトゲノム（ＧＲＣｈ３８）（非特許文献３８、３９）にマッピングした。具体的には、これらのツールを高速計算機に実装し、バッチ処理を実施するためのスクリプトを作成することで中間ファイルを介した連続処理を実施した。
　バッチ１とバッチ２のデータは、ｄｅ　ｄｅｖｏ変異を共有していない、と仮定した。すなわち、２人の匿名患者は、共通の（生殖細胞系列の）突然変異を共有していない、すなわち、彼らは互いに関連していない、と仮定した。Ｓａｍｔｏｏｌｓを使用して変異を検出した（非特許文献４２）。具体的には、参照ゲノムとリード配列の比較により、細胞にまたがる多重アラインメントを作成し、変異の位置と数を推定した。
　本明細書の実施例では、検出されたすべてのｉｎｄｅｌイベントを除いて、一塩基の変異のみを使用した。また、全細胞で５０％以上の不完全な部位データが検出された場合、即ち、全細胞のうち５０％以上の細胞において、トランクリプトームのＲＮＡ配列データが得られない部位については、変異は生じていないと推定し、該当する部位をスクリーニングから除外した。
　（１－３）単一細胞のトランスクリプトーム変異の検出
　明らかな変異（マッピングされた転写産物における参照ゲノムからの逸脱）には、生殖細胞変異と体細胞変異の両方が含まれる可能性がある。ハプロタイプの位相に関する情報は提供されていないため、細胞集団の多反復部位における体細胞変異を検出するための簡便な方法を採用した。１つの系統に最大でも１つの体細胞変異が発生すると仮定し、細胞全体で起こりうる突然変異パターンを分類した（非特許文献４３，４４）
　図２に、本発明の方法の一態様の概念図を示す。本発明の方法は、以下の工程を含む。
　（１）複数の単一細胞について、各単一細胞のトランクリプトームのＲＮＡ配列データを、対応するゲノム配列データと比較し、ＲＮＡ配列とゲノム配列が異なる部位をスクリーニングする；
　（２）工程（１）のスクリーニングで特定された部位のデータを細胞間で比較し、３つ若しくはそれ以上の１塩基多型が検出される部位を、体細胞変異の部位としてスクリーニングする。
　２つの変異部位だけでは、その変異が体細胞変異又は生殖細胞変異のいずれが生じたのか、あるいはその部位がヘテロ接合なのかを知ることができない（図２ａ、図２ｂ、図２ｄ、図２ｇ、図２ｈ）。しかし、３つ若しくはそれ以上の多発部位が観察された場合、少なくとも１つの体細胞変異が細胞系列のどこかで起こった、ということが説明できる（図２ｃ、図２ｅ）。
　一方、どの系統で突然変異が起こったかは、そのままでは不明である。一つのアイデアとして、観察されたヌクレオチドが参照部位を共有している場合、参照ゲノムの状態を祖先の状態として用いることできる。しかし、我々の枠組みでは、多遺伝子部位は常にヘテロ接合であり、参照ゲノムデータはハプロタイプフェーズされていない。ここでは、マイナーアリルが細胞集団の中で新たに派生したヌクレオチドであると仮定し、単純に派生変異として「マイナーアリル」を選択することにした。選択したマイナーアリルは、本実施例において対立遺伝子において見出される変異のうち、最も頻度が低かった体細胞変異である。
　検出された変異部位は、ＳｎｐＥｆｆソフトウェア（非特許文献４５）を用いてアノテーションを行った。具体的には、トランスクリプト配列をゲノムにマップすることで取得したＶＣＦファイルを入力データとして与え、使用したゲノムデータを指定することでアノテーションを付与した。
　（１－４）結果
　ＳＲＰ０９０９４４のデータ（バッチ１、５４細胞）のマッピング
　ＳＲＰ０９０９４４の３，０８８，２８６，４０１ｂｐの転写産物のデータ（バッチ１：５４細胞）を、参照ゲノム（ＧＲＣｈ３８）（非特許文献３８，３９）にマッピングした。
　平均カバー率は０．６８５％（ＳＤ：０．２３１）であった（図４、バッチ１の細胞）。最初のスクリーニングで参照ゲノムと異なる１，９６５，６２９個の部位を検出した。また、５４個の単一細胞のすべてで観察される、３種類以上のヌクレオチドが変異している多変異部位は８９箇所であった。また、平均４３．２個の単一細胞（５４個の単一細胞の８０％）において観察される多変異部位は、２，０８３個存在した。データはすべて品質管理されたものである。
　ＳＲＰ０９０９４４のデータ（バッチ２、３３細胞）
　平均カバー率は０．４７７％（ＳＤ：０．２４３）であった（図４、バッチ２の細胞）。１回目のスクリーニングで参照ゲノムと異なる８３０，９０５部位を検出した。３３個の単一細胞のすべてで観察される多変異部位は５３箇所であった。また、平均２６．４個の単一細胞（３３個の単一細胞の８０％）において観察される多変異部位は５７４箇所存在した。
　変異の注釈（ａｎｎｏｔａｔｉｏｎ）
　表１および表２は、ＳｎｐＥｆｆソフトウェア（非特許文献４５）によってアノテーションされた推定変異を示す。
　表１　ＳｎｐＥｆｆソフトウェア（ＳＲＰ０９０９４４のデータ、バッチ１、５４細胞）による推定変異型
Figure JPOXMLDOC01-appb-T000001
 注：変異の数は、各カテゴリーで排他的ではない。
　表２　　ＳｎｐＥｆｆソフトウェア（ＳＲＰ０９０９４４のデータ、バッチ２、３３細胞）による推定変異型
Figure JPOXMLDOC01-appb-T000002
 注：変異の数は、各カテゴリーで排他的ではない。
　表１、表２の結果には、例えば、ネストしたイントロン遺伝子（非特許文献６３）などの代替的な注釈が含まれている。一部の変異数はカテゴリー間で重複している。例えば、８０％基準でバッチ１とバッチ２のデータからそれぞれ２０８３箇所と５７４箇所の変異を検出した。しかしながら、ＳｎｐＥｆｆソフトウェアはデフォルトパラメータセットでバッチ１とバッチ２のデータからそれぞれ１９０３箇所と１３９８箇所の変異を推定している。独立したソフトウェアで推定した両変異箇所数のオーダーがほぼ揃っていることは、本解析で用いたアノテーションの妥当性を示している。
　ＳｎｐＥｆｆの結果、バッチ１とバッチ２のデータから、それぞれ５５０個と２１６個のミスセンス変異、１９９個と１３５個の同義変異があることが判明した。
　実施例２　細胞の系統解析及び細胞系統樹の作成
　（２－１）細胞の系統解析
　実施例１で取得した、観測可能な変異部位をすべて連結し、配列アライメントを作成した。「変異部位の連結」とは、変異を含むコドン配列を作成するため、変異の近傍のコーディング領域の配列を変異部位と結合し３塩基のコドンを作成することを意味する。ＭＥＧＡ　Ｘ（非特許文献４８）に実装されている最大節約法（ｍａｘｉｍｕｍ　ｐａｒｓｉｍｏｎｙ　ｍｅｔｈｏｄ）（非特許文献４６、４７）を用いて、デフォルトのパラメータで細胞系統樹を再構築した。
　具体的には、作成した多重アラインメントをＭＥＧＡ　Ｘに入力データとして与え、ＧＵＩから最大節約法を選択することで系統樹の再構築を行なった。結果をＮｅｗｉｃｋ　ｔｒｅｅ形式（非特許文献４９）で出力し、以降の処理に使用した。
　（２－２）遺伝子発現の再解析
　線形次元削減（ＬＤＲ）のために、Ｒ（バージョン３．６．２）を用いて遺伝子発現パターンの主成分分析（ＰＣＡ）を行った（非特許文献５４）。具体的には、単一細胞の遺伝子発現パターンは単一細胞の種類と遺伝子の種類からなる遺伝子発現マトリックスで表現されるが、このうち遺伝子の種類方向の成分を低次元空間に投射して、遺伝子発現マトリックスの次元を削減した。
　その後、ｔ－ＳＮＥ（非特許文献３２、３３）とＵＭＡＰ（非特許文献３４、３５）を適用し、非線形次元削減（ＮＤＲ）を行った。クラスタリングにはＬｏｕｖａｉｎ法（非特許文献５５）を用いた。具体的には、主成分分析により大域的な線形次元圧縮を実施した後、遺伝子発現マトリックスをそれぞれのソフトウェアに入力データとして与えることで行った。
　（２－３）変異パターンと遺伝子発現プロファイルの比較解析
　クラスター化した細胞を、細胞の遺伝子型に基づいて再構築した細胞系統樹にマッピングした。そのために，Ｐｈｙｌｏｇｅｎｅｔｉｃｓ　ｆｏｒ　Ｍａｔｈｅｍａｔｉｃａ（非特許文献６０）とＰｈｙｌｏｇｅｎｅｔｉｃｓ（非特許文献６１）を適用して、Ｍａｔｈｅｍａｔｉｃａ（非特許文献５９）のコード、ＡｓｓｉｇｎＣｌｕｓｔｅｒ２Ｃｅｌｌを開発した。
　「Ｍａｔｈｅｍａｔｉｃａ」は数式処理ソフトウェアである。「Ｐｈｙｌｏｇｅｎｅｔｉｃｓ　ｆｏｒ　Ｍａｔｈｅｍａｔｉｃａ」は、数式処理ソフトウェア用のライブラリーである。「ＡｓｓｉｇｎＣｌｕｓｔｅｒ２Ｃｅｌｌ」は、本実施例において、細胞系統樹と遺伝子発現プロファイルを統合するために開発したプログラムの名称である。
　ＡｓｓｉｇｎＣｌｕｓｔｅｒ２Ｃｅｌｌは、再構築された細胞系統樹をＮｅｗｉｃｋ形式で読み込み、これと事前にクラスタリングされた細胞集団を比較することで、それらの一致度を円グラフとして各ノード上に表示する。樹の読み込みや視覚化にはＰｈｙｌｏｇｅｎｅｔｉｃｓ　ｆｏｒ　ＭａｔｈｅｍａｔｉｃａとＰｈｙｌｏｇｅｎｅｔｉｃｓで提供されている関数群をライブラリーとして用いている。ＡｓｓｉｇｎＣｌｕｓｔｅｒ２Ｃｅｌｌ自身もＭａｔｈｅｍａｔｉｃａで記述された関数であり、Ｍａｔｈｅｍａｔｉｃａ　Ｎｏｔｅｂｏｏｋ上で動作するよう最適化されている。
　ＡｓｓｉｇｎＣｌｕｓｔｅｒ２Ｃｅｌｌを用いて、Ｎｅｗｉｃｋ（系統樹を表現するための標準的なデータ型式）（非特許文献６２）形式の樹形ファイルと遺伝子発現プロファイルのクラスターテーブルを読み込み、クラスターＩＤと樹形図を対応付けし、各ノードの単一系統性度（ＤｏＭ）をそのサブツリーのクラスター数の割合で評価した。
　結果を図３に示す。図３の円グラフは各ノードのＤｏＭを表している。例えば、ノードｎの円グラフが１枚しかない場合、サブツリーｎは完全に単系統、ノードｎ’の円グラフが２枚ある場合、サブツリーｎ’は多系統（図３ではサブツリーＢ）であることがわかる。サブツリーＢでは細胞型２が優勢であるため、細胞型１は細胞型２から派生したと推測される。このように、遺伝子発現プロファイルと細胞系譜をもとに、細胞型の関連性を描き出した。
　理論的には、細胞系統樹の根は接合細胞（受精卵）を表し、そのインスタンスは受精卵の接合体ゲノムを表す。しかし、実施例では、推定された細胞系統樹の根は、観察された細胞集団の前駆細胞を表している可能性がある。接合細胞は、細胞系統樹の根と参照ゲノムとの間のどこかに存在する。
　（２－４）細胞系統樹
　２つの異なる個体の胎盤組織の単一細胞の系統樹を再構築した（図５、図６）。ＭＥＧＡ－Ｘ（非特許文献４８）に実行させた最大節約法（非特許文献４６、４７）を使用したため、枝の長さは枝上で発生した変異の期待値を表している。なお、各葉（ｏｐｅｒａｔｉｏｎａｌ　ｔａｘｏｎｏｍｉｃ　ｕｎｉｔ，ＯＴＵ（非特許文献６４））は、各組織からサンプリングされた細胞を表している。
　細胞栄養芽層（ＣＹＴ：図３、４のＣＹＴ１、ＣＹＴ２、ＣＹＴ３）と絨毛外栄養膜（ＥＶＴ）の細胞系列は、Ｎｏｔｃｈ／Ｗｎｔシグナル（非特許文献６５）との関連で上手くモデル化されている。すなわち、ＥＶＴは細胞柱栄養細胞（ＣＣＴ）を介してＣＹＴから分化する（図９）。本実施例の結果は、非特許文献６５に記載のＮｏｔｃｈ／Ｗｎｔシグナルのモデルと一致する。一方，Ｐａｖｌｉｃｅｖらは、自分たちのデータには母体脱落膜細胞（ＤＣ）が含まれていると述べている．しかし、本実施例の結果は、推定されるＤＣ細胞は胎児ＣＹＴまたはＣＹＴの幹細胞から分化したものであることを示唆している（図５及び図６）。また、ｔ－ＳＮＥ（非特許文献３２、３３）を用いたトランスクリプトームデータの再解析でも、本実施例の結果が正しいことが支持された（図７）。
　体細胞変異は、基本的に細胞集団の中で生き残るために細胞にマイナスの影響を与えるので、正常な組織では本物の体細胞変異は純化選択を受けていると仮定することは合理的である（非特許文献５０，５１）。例えば、最尤法によれば、全体として体細胞変異は浄化選択を受けており、検出された体細胞変異の大部分は、少なくとも正の偽とは程遠いことが示唆される。
　一方、本実施例において、異なる個体から採取した２つのサンプルには、少数の正の選択部位も検出された。この結果は、体細胞変異が遺伝的修飾因子にも関与し、重要な役割を担っていることを示唆している可能性がある。
　実施例３　コーディング領域におけるｄＮ／ｄＳ比の評価
　本実施例では、コーディング領域におけるｄＮ／ｄＳ比を評価した。
　正常な組織では、本物の体細胞変異は純化選択を受けていることが予想される（非特許文献５０，５１）。本発明の方法の検出の信頼性を評価するために、実施例１で検出された変異体の選択圧をｄＮ／ｄＳ比の観点から評価した。検出された変異体が入るコドン配列を組み立て、エキソン変異体を含むコドンアラインメントを作成した。Ｐａｍｌパッケージ（非特許文献５２）のＣｏｄｅｍｌを用いて、全体のｄＮ／ｄＳ比を計算した。
　具体的には、再構成した細胞系統樹をガイドツリーとして用いて、ＡｎａｃｏｎｄａによりＰａｌｍ環境を生成し、そこに必要なパッケージを実装した後、Ｊｕｐｙｔｅｒ　Ｎｏｔｅｂｏｏｋ上でＰｙｔｈｏｎコードを用いてＰａｌｍのＣｏｄｅｍｌモジュールを呼び出してコドンから作成した多重アラインメントを処理した。また、Ｎｅｉ＆Ｇｏｊｏｂｏｒｉ法（非特許文献５３）を用いて、細胞の遺伝子型間の一対のｄＮ／ｄＳ比を計算した。具体的には、上記Ｃｏｄｅｍｌモジュールを呼び出してコドンから作成した多重アラインメントを処理した。
　ＳＲＰ０９０９４４データのバッチ１およびバッチ２の全体のｄＮ／ｄＳ比は、それぞれ０．８６５および０．５５６であった。この結果は、２つの組織における体細胞変異が、純化選択を受けていたことを示唆している。
　一方、ＳＲＰ０９０９４４データバッチ１では、Ｎａｉｖｅ　Ｅｍｐｉｒｉｃａｌ　Ｂａｙｅｓ（ＮＥＢ）解析とＢａｙｅｓ　Ｅｍｐｉｒｉｃａｌ　Ｂａｙｅｓ（ＢＥＢ）解析の両方により、特定の変異部位でのｄＮ／ｄＳ比をωとすると、Ｐｒ（ω＞１）＝１という正の選択を受けている変異体部位を検出した（非特許文献６６）。この特定の変異部位は、ＮＭ＿００１１４４９６４．１，Ｈｏｍｏ　ｓａｐｉｅｎｓ　ＮＥＤＤ４　ｌｉｋｅ　Ｅ３　ｕｂｉｑｕｉｔｉｎ　ｐｒｏｔｅｉｎ　ｌｉｇａｓｅ　（ＮＥＤＤ４Ｌ），ｔｒａｎｓｃｒｉｐｔ　ｖａｒｉａｎｔ　ｂにある（非特許文献６７）。
　また、ＳＲＰ０９０９４４のデータバッチ２では、ＮＥＢ解析でＰｒ（ω＞１）＝０．９２８、ＢＥＢ解析で０．５６７の変異部位が見つかった。この変異部位は、ＮＭ＿　００１１７２８９５．１，　Ｈｏｍｏ　ｓａｐｉｅｎｓ　ｃａｖｅｏｌｉｎ　１　（ＣＡＶ１），　ｔｒａｎｓｃｒｉｐｔ　ｖａｒｉａｎｔ　２（非特許文献６８）の遺伝子座にある。
　ＳＲＰ０９０９４４データバッチ１とバッチ２の細胞のすべての可能なペア１４８５個と５６１個のうち、ｄＮ／ｄＳ比が１より大きい系統は、それぞれ４８２個と７１個であった。したがって、細胞間の系統の６７．５％と８７．３％は、それぞれ浄化選択または中立の影響を受けていることが示唆された。
　比較例１　擬似時間経過解析
　本比較例では、Ｒ（ｖｅｒｓｉｏｎ　４．１．２）（非特許文献５４）のｍｏｎｏｃｌｅ３（ｖｅｒｓｉｏｎ　１．０．０）（非特許文献５６－５８）を用いて、ＳＲＰ０９０９４４データバッチ１の擬似時間経過解析を行った。具体的には、次元圧縮（２６次元）した単一細胞の遺伝子発現マトリックスを入力データとして与え、各単一細胞に分化の様相を反映すると期待される相対的な擬似時間を割り当てた。
　結果の可視化には、ＭａｔｈｅｍａｔｉｃａのＬｉｓｔＣｏｎｔｏｕｒＰｌｏｔ関数を用いた（非特許文献５９）。具体的には、２次元平面上に投射された単一細胞を、割り当てられた擬似時間で３次補完することによって相対時間的な等高線を得た。
　結果を図８に示す。細胞の擬似時間は、ＭａｔｈｅｍａｔｉｃａのＬｉｓｔＣｏｎｔｏｕｒＰｌｏｔ関数（非特許文献５９）によって補間次数３で作成された等高線で表されている。この結果は，擬似時間の方向を除いて，本発明の方法を用いた結果とかなり一致していた。例えば、擬似時間が２．５より大きい細胞の数は１６個（ＳＲＲ４３７１５２７，ＳＲＲ４３７１５３１，ＳＲＲ４３７１５３２，ＳＲＲ４３７１５３３，ＳＲＲ４３７１５３６，ＳＲＲ４３７１５４２，ＳＲＲ４３７１５４７，ＳＲＲ４３７１５６３，ＳＲＲ４３７１５６６，ＳＲＲ４３７１５６８，ＳＲＲ４３７１５６９，ＳＲＲ４３７１５７０，ＳＲＲ４３７１５７１，ＳＲＲ４３７１５７２，ＳＲＲ４３７１５７５及びＳＲＲ４３７１５７７），であった。そのうちの１０個は、我々の解析によれば、ＣＹＴ細胞として推定幹細胞の自己再生段階から直接由来している（図５ｃ）。
　実施例の評価
　実施例において、低域通過の単一細胞ＲＮＡ配列データを用いて、ヒト胎盤組織の体細胞の系統樹を再構築できることが示された。この系統樹は、細胞栄養芽層（ＣＹＴ）Ｉ、ＣＹＴ　ＩＩ、ＣＹＴ　ＩＩＩ、絨毛外栄養膜（ＥＶＴ）の４種類の細胞に関して、既知の胎盤の細胞系統と一致している。異種分化した細胞を１回のタイムスナップで使用する擬似時間経過解析とは異なり、本発明の方法は、体細胞変異を使用して、子孫細胞まで遡って細胞系譜を追跡するものである。したがって、推定された子孫細胞は、再構成された細胞軌跡樹（ｃｅｌｌｅｒ　ｔｒａｊｅｃｔｏｒｙ　ｔｒｅｅ）の内部ノード（頂点）として表現される。
　得られた細胞の「系統的」特徴は、２個体からサンプリングした胎盤における絨毛幹細胞および絨毛外細胞分化におけるＮｏｔｃｈ／Ｗｎｔシグナルの役割を統合したモデル系と一致した。このことは、低域通過の配列データであっても、細胞系統樹を遡及的に推定することが可能であることを示唆している。個々の体細胞変異の検出は超深度シーケンスがなければ困難であると考えられていたが、本発明はそれを可能にする。また、体細胞変異を持つ細胞は予想通り純化選択を受けているが、少なくとも特定の部位では適応的進化のシグナルがあることがわかった。
　さらに、検出された体細胞変異の質を進化論的な枠組みで評価した。検出された体細胞変異は、全体として純化選択（ｄＮ／ｄＳ＜１）を受けているが、少なくとも特定の部位では適応的進化（ｄＮ／ｄＳ＞１）のサインも示していることが確認された。これらの進化的洞察は、少なくともコーディング領域における本発明の結果の中程度の信頼性を支持するものである。さらに、実施例の結果は、マップされたリードの配列とその深さという、全く異なる種類のデータを用いた擬似時間経過解析の結果と一致するものであった。
　単一細胞遺伝子発現解析において、低域通過シーケンシングデータ（トランスクリプトームのＲＮＡ配列データ）は「副産物」であり、単一細胞間の多型（モザイク）情報はしばしば無視されてきた。本発明は、落ち穂拾い的な方法で有意な遺伝子型情報を抽出することが可能であることを示した。本発明の、単一細胞トランスクリプトームデータの体細胞変異解析は、高次元の遺伝子発現データの解釈を可能にする。さらに、体細胞変異の生物学的意義は、個体内の「進化」という新たな視点に光を当てる。

Claims (11)

  1. 　細胞の系統解析を行う方法であって、
    　（１）同一個体由来の複数の単一細胞について、各単一細胞のトランクリプトームのＲＮＡ配列データを、対応するゲノム配列データと比較し、ＲＮＡ配列とゲノム配列が異なる部位をスクリーニングし；
    　（２）工程（１）のスクリーニングで特定された部位のうち、３つ若しくはそれ以上の１塩基多型が検出される部位を、体細胞変異の部位としてスクリーニングし；そして、
    　（３）工程（２）で得られた体細胞変異の部位の塩基配列情報に基づき、細胞系統樹を作成する、
    ことを含む、前記方法。
  2. 　（４）前記細胞系統樹から、各単一細胞の細胞の種類を推定する、ことをさらに含む、請求項１に記載の方法。
  3. 　各単一細胞の遺伝子発現プロファイルから推定される単一細胞の種類の情報を、前記細胞系統樹の情報と比較する、ことをさらに含む、請求項１又は２の方法。
  4. 　１塩基多型が、１塩基置換である、請求項１－３のいずれか１項に記載の方法。
  5. 　工程（１）において、全細胞のうち５０％以上の細胞において、トランクリプトームのＲＮＡ配列データが得られない部位については、変異は生じていないと推定し、スクリーニングから除外する、請求項１－４のいずれか１項に記載の方法。
  6. 　工程（２）でスクリーニングする体細胞変異が、頻度の少ない体細胞変異である、請求項１－５のいずれか１項に記載の方法。
  7. 　疾患又は症状と関連する体細胞変異を検出するための、請求項１－６のいずれか１項に記載の方法。
  8. 　疾患又は症状が、癌、認知症、循環器疾患、加齢、自己免疫疾患、神経変性疾患、及び精神疾患からなる群から選択される、請求項７に記載の方法。
  9. 　請求項１－８のいずれか１項に記載の方法を実行するためのプログラム。
  10. 　請求項１－８のいずれか１項に記載の方法を実行するためのプログラムを記録した記録媒体。
  11. 　プロセッサと、前記プロセッサで実行されたときに、請求項１－８のいずれか１項に記載の方法を実行するプログラムを記録したメモリーを有するシステム。
     
     
     
PCT/JP2023/005534 2022-02-18 2023-02-16 細胞の系統解析を行う方法 WO2023157933A1 (ja)

Applications Claiming Priority (2)

Application Number Priority Date Filing Date Title
JP2022023985 2022-02-18
JP2022-023985 2022-02-18

Publications (1)

Publication Number Publication Date
WO2023157933A1 true WO2023157933A1 (ja) 2023-08-24

Family

ID=87578654

Family Applications (1)

Application Number Title Priority Date Filing Date
PCT/JP2023/005534 WO2023157933A1 (ja) 2022-02-18 2023-02-16 細胞の系統解析を行う方法

Country Status (1)

Country Link
WO (1) WO2023157933A1 (ja)

Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012505666A (ja) * 2008-10-20 2012-03-08 ファーマテスト・サービシーズ・リミテッド Ｎａｖ３の遺伝子異常および複数遺伝子の異常発現を含む方法および使用
US20140206006A1 (en) * 2011-08-25 2014-07-24 Bgi Shenzhen Single cell classification method, gene screening method and device thereof
JP2021179867A (ja) * 2020-05-15 2021-11-18 株式会社D’isum ゲノム分析装置及び方法

Patent Citations (3)

* Cited by examiner, † Cited by third party
Publication number Priority date Publication date Assignee Title
JP2012505666A (ja) * 2008-10-20 2012-03-08 ファーマテスト・サービシーズ・リミテッド Ｎａｖ３の遺伝子異常および複数遺伝子の異常発現を含む方法および使用
US20140206006A1 (en) * 2011-08-25 2014-07-24 Bgi Shenzhen Single cell classification method, gene screening method and device thereof
JP2021179867A (ja) * 2020-05-15 2021-11-18 株式会社D’isum ゲノム分析装置及び方法

Non-Patent Citations (1)

* Cited by examiner, † Cited by third party
Title
"Panoramic View Report: Life Science and Clinical Research Field (2021). Center for Research and Development Strategy", 1 March 2021, JAPAN SCIENCE AND TECHNOLOGY AGENCY, JP, ISBN: 978-4-88890-711-8, article "Passages", pages: 230 - 520-538, XP009548948 *

Similar Documents

Publication Publication Date Title
US20200251180A1 (en) Resolving genome fractions using polymorphism counts
JP2020058393A (ja) 母体血漿の無侵襲的出生前分子核型分析
JP7009518B2 (ja) 既知又は未知の遺伝子型の複数のコントリビューターからのｄｎａ混合物の分解及び定量化のための方法並びにシステム
AU2014314151B2 (en) Haplotyping and copy number typing using polymorphic variant allelic frequencies
JP2019512122A (ja) 異常な核型の検出のための方法およびシステム
Hård et al. Conbase: a software for unsupervised discovery of clonal somatic mutations in single cells through read phasing
CN110770839A (zh) 来自未知基因型贡献者的dna混合物的精确计算分解的方法
JP7333838B2 (ja) 胚における遺伝パターンを決定するためのシステム、コンピュータプログラム及び方法
Fortier et al. Detection of CNVs in NGS data using VS-CNV
WO2023157933A1 (ja) 細胞の系統解析を行う方法
US10424395B2 (en) Computation pipeline of single-pass multiple variant calls
US20220059185A1 (en) Method and apparatus for detecting copy number variations in a genome
Carrami et al. A highly accurate platform for clone-specific mutation discovery enables the study of active mutational processes
Niehus et al. PopDel identifies medium-size deletions jointly in tens of thousands of genomes
Huang et al. The extent of linkage disequilibrium and computational challenges of single nucleotide polymorphisms in genome-wide association studies
Dharanipragada et al. Copy number variation detection workflow using next generation sequencing data
Oota et al. Real-Time Course: Reconstruction of cellular diversity and lineage trajectory based on somatic mutational patterns detected from low-pass single-cell transcriptome data
Zhang Detection copy number variants profile by multiple constrained optimization
Hård et al. Conbase: a software for discovery of clonal somatic mutations in single cells through read phasing
KEPLER Analyzing NGS Data with Machine Learning-From IBD Segments to Copy Number Variations
Hedges Bioinformatics of Human Genetic Disease Studies
Deveau et al. Clonal assessment of functional mutations in cancer based on a genotype-aware method for clonal reconstruction
Laiho Data analysis tools and methods for DNA microarray and high-throughput sequencing data
Goldstein et al. Calling Rare Variants from Genotype Data
Heinrich Aspects of Quality Control for Next Generation Sequencing Data in Medical Genetics

Legal Events

Date Code Title Description
121 Ep: the epo has been informed by wipo that ep was designated in this application

Ref document number: 23756457

Country of ref document: EP

Kind code of ref document: A1

DPE1 Request for preliminary examination filed after expiration of 19th month from priority date (pct application filed from 20040101)