Estimating gene-level false discovery probability improves eQTL statistical fine-mapping precision.

Statistical fine-mapping prioritizes putative causal variants from a large number of candidate variants, and is widely used in expression quantitative loci (eQTLs) studies. In eQTL fine-mapping, the existence of causal variants for gene expression is not guaranteed, since the genetic heritability of gene expression explained by nearby (cis-) variants is limited. Here we introduce a refined fine-mapping algorithm, named Knockoff-Finemap combination (KFc). KFc estimates the probability that the causal variant(s) exist in the cis-window of a gene through construction of knockoff genotypes (i.e. a set of synthetic genotypes that resembles the original genotypes), and uses it to adjust the posterior inclusion probabilities (PIPs). Utilizing simulated gene expression data, we show that KFc results in calibrated PIP distribution with improved precision. When applied to gene expression data of 465 genotyped samples from the Japan COVID-19 Task Force (JCTF), KFc resulted in significant enrichment of a functional score as well as reporter assay hits in the top PIP bins. When combined with functional priors derived from an external fine-mapping study (GTEx), KFc resulted in a significantly higher proportion of hematopoietic trait putative causal variants in the top PIP bins. Our work presents improvements in the precision of a major fine-mapping algorithm.

Subclonal accumulation of immune escape mechanisms in microsatellite instability-high colorectal cancers.

BACKGROUND: Intratumor heterogeneity (ITH) in microsatellite instability-high (MSI-H) colorectal cancer (CRC) has been poorly studied. We aimed to clarify how the ITH of MSI-H CRCs is generated in cancer evolution and how immune selective pressure affects ITH. METHODS: We reanalyzed public whole-exome sequencing data on 246 MSI-H CRCs. In addition, we performed a multi-region analysis from 6 MSI-H CRCs. To verify the process of subclonal immune escape accumulation, a novel computational model of cancer evolution under immune pressure was developed. RESULTS: Our analysis presented the enrichment of functional genomic alterations in antigen-presentation machinery (APM). Associative analysis of neoantigens indicated the generation of immune escape mechanisms via HLA alterations. Multiregion analysis revealed the clonal acquisition of driver mutations and subclonal accumulation of APM defects in MSI-H CRCs. Examination of variant allele frequencies demonstrated that subclonal mutations tend to be subjected to selective sweep. Computational simulations of tumour progression with the interaction of immune cells successfully verified the subclonal accumulation of immune escape mutations and suggested the efficacy of early initiation of an immune checkpoint inhibitor (ICI) -based treatment. CONCLUSIONS: Our results demonstrate the heterogeneous acquisition of immune escape mechanisms in MSI-H CRCs by Darwinian selection, providing novel insights into ICI-based treatment strategies.

Convergent genomic diversity and novel BCAA metabolism in intrahepatic cholangiocarcinoma.

BACKGROUND: Driver alterations may represent novel candidates for driver gene-guided therapy; however, intrahepatic cholangiocarcinoma (ICC) with multiple genomic aberrations makes them intractable. Therefore, the pathogenesis and metabolic changes of ICC need to be understood to develop new treatment strategies. We aimed to unravel the evolution of ICC and identify ICC-specific metabolic characteristics to investigate the metabolic pathway associated with ICC development using multiregional sampling to encompass the intra- and inter-tumoral heterogeneity. METHODS: We performed the genomic, transcriptomic, proteomic and metabolomic analysis of 39-77 ICC tumour samples and eleven normal samples. Further, we analysed their cell proliferation and viability. RESULTS: We demonstrated that intra-tumoral heterogeneity of ICCs with distinct driver genes per case exhibited neutral evolution, regardless of their tumour stage. Upregulation of BCAT1 and BCAT2 indicated the involvement of 'Val Leu Ile degradation pathway'. ICCs exhibit the accumulation of ubiquitous metabolites, such as branched-chain amino acids including valine, leucine, and isoleucine, to negatively affect cancer prognosis. We revealed that this metabolic pathway was almost ubiquitously altered in all cases with genomic diversity and might play important roles in tumour progression and overall survival. CONCLUSIONS: We propose a novel ICC onco-metabolic pathway that could enable the development of new therapeutic interventions.

DOCK2 is involved in the host genetics and biology of severe COVID-19.

Identifying the host genetic factors underlying severe COVID-19 is an emerging challenge1-5. Here we conducted a genome-wide association study (GWAS) involving 2,393 cases of COVID-19 in a cohort of Japanese individuals collected during the initial waves of the pandemic, with 3,289 unaffected controls. We identified a variant on chromosome 5 at 5q35 (rs60200309-A), close to the dedicator of cytokinesis 2 gene (DOCK2), which was associated with severe COVID-19 in patients less than 65 years of age. This risk allele was prevalent in East Asian individuals but rare in Europeans, highlighting the value of genome-wide association studies in non-European populations. RNA-sequencing analysis of 473 bulk peripheral blood samples identified decreased expression of DOCK2 associated with the risk allele in these younger patients. DOCK2 expression was suppressed in patients with severe cases of COVID-19. Single-cell RNA-sequencing analysis (n = 61 individuals) identified cell-type-specific downregulation of DOCK2 and a COVID-19-specific decreasing effect of the risk allele on DOCK2 expression in non-classical monocytes. Immunohistochemistry of lung specimens from patients with severe COVID-19 pneumonia showed suppressed DOCK2 expression. Moreover, inhibition of DOCK2 function with CPYPP increased the severity of pneumonia in a Syrian hamster model of SARS-CoV-2 infection, characterized by weight loss, lung oedema, enhanced viral loads, impaired macrophage recruitment and dysregulated type I interferon responses. We conclude that DOCK2 has an important role in the host immune response to SARS-CoV-2 infection and the development of severe COVID-19, and could be further explored as a potential biomarker and/or therapeutic target.

The whole blood transcriptional regulation landscape in 465 COVID-19 infected samples from Japan COVID-19 Task Force.

Coronavirus disease 2019 (COVID-19) is a recently-emerged infectious disease that has caused millions of deaths, where comprehensive understanding of disease mechanisms is still unestablished. In particular, studies of gene expression dynamics and regulation landscape in COVID-19 infected individuals are limited. Here, we report on a thorough analysis of whole blood RNA-seq data from 465 genotyped samples from the Japan COVID-19 Task Force, including 359 severe and 106 non-severe COVID-19 cases. We discover 1169 putative causal expression quantitative trait loci (eQTLs) including 34 possible colocalizations with biobank fine-mapping results of hematopoietic traits in a Japanese population, 1549 putative causal splice QTLs (sQTLs; e.g. two independent sQTLs at TOR1AIP1), as well as biologically interpretable trans-eQTL examples (e.g., REST and STING1), all fine-mapped at single variant resolution. We perform differential gene expression analysis to elucidate 198 genes with increased expression in severe COVID-19 cases and enriched for innate immune-related functions. Finally, we evaluate the limited but non-zero effect of COVID-19 phenotype on eQTL discovery, and highlight the presence of COVID-19 severity-interaction eQTLs (ieQTLs; e.g., CLEC4C and MYBL2). Our study provides a comprehensive catalog of whole blood regulatory variants in Japanese, as well as a reference for transcriptional landscapes in response to COVID-19 infection.

Immunogenomic pan-cancer landscape reveals immune escape mechanisms and immunoediting histories.

Immune reactions in the tumor microenvironment are an important hallmark of cancer, and emerging immune therapies have been proven effective against several types of cancers. To investigate cancer genome-immune interactions and the role of immunoediting or immune escape mechanisms in cancer development, we analyzed 2834 whole genome and RNA sequencing datasets across 31 distinct tumor types with respect to key immunogenomic aspects and provided comprehensive immunogenomic profiles of pan-cancers. We found that selective copy number changes in immune-related genes may contribute to immune escape. Furthermore, we developed an index of the immunoediting history of each tumor sample based on the information of mutations in exonic regions and pseudogenes and evaluated the immunoediting history of each tumor. Our immuno-genomic analyses of pan-cancers have the potential to identify a subset of tumors with immunogenicity and diverse backgrounds or intrinsic pathways associated with their immune status and immunoediting history.

The Evolving Genomic Landscape of Esophageal Squamous Cell Carcinoma Under Chemoradiotherapy.

Esophageal squamous cell carcinoma (ESCC) often recurs after chemoradiotherapy, and the prognosis of ESCC after chemoradiotherapy has not improved over the past few decades. The mutation process in chemoradiotherapy-resistant clones and the functional relevance of genetic alterations remain unclear. To address these problems, we performed whole-exome sequencing of 52 tumor samples from 33 patients with ESCC who received radiotherapy combined with 5-fluorouracil/platinum. In multiregion analyses of pretreatment and locally recurrent lesions from five cases, most driver gene-altered clones remained under chemoradiotherapy selection pressure, while few driver gene alterations were acquired at recurrence. The mutation signatures of recurrent ESCC, including increased deletion frequency and platinum dose-dependent base substitution signatures, were substantially different from those of primary ESCC and reflected the iatrogenic impacts of chemoradiotherapy. Single-region analysis of 28 pretreatment tumors indicated that focal copy-number gain at the MYC locus was significantly associated with poor progression-free survival and overall survival after chemoradiotherapy. MYC gain remained throughout the chemoradiotherapy course and potentially contributes to intrinsic resistance to chemoradiotherapy. Consistent with these findings, MYC copy number and mRNA and protein levels in ESCC cell lines correlated positively with resistance to radiotherapy, and MYC knockdown improved sensitivity to radiotherapy. Overall, these data characterize the clonal evolution process induced by chemoradiotherapy and clinically relevant associations for genetic alterations in ESCC. These findings increase our understanding of therapeutic resistance and support the rationale for precision chemoradiotherapy. SIGNIFICANCE: Whole-exome sequencing reveals the genetic evolution of ESCC during chemoradiotherapy, highlighting MYC gain in pretreatment tumors as a potential marker of therapy resistance.

Deficiency of TMEM53 causes a previously unknown sclerosing bone disorder by dysregulation of BMP-SMAD signaling.

Bone formation represents a heritable trait regulated by many signals and complex mechanisms. Its abnormalities manifest themselves in various diseases, including sclerosing bone disorder (SBD). Exploration of genes that cause SBD has significantly improved our understanding of the mechanisms that regulate bone formation. Here, we discover a previously unknown type of SBD in four independent families caused by bi-allelic loss-of-function pathogenic variants in TMEM53, which encodes a nuclear envelope transmembrane protein. Tmem53-/- mice recapitulate the human skeletal phenotypes. Analyses of the molecular pathophysiology using the primary cells from the Tmem53-/- mice and the TMEM53 knock-out cell lines indicates that TMEM53 inhibits BMP signaling in osteoblast lineage cells by blocking cytoplasm-nucleus translocation of BMP2-activated Smad proteins. Pathogenic variants in the patients impair the TMEM53-mediated blocking effect, thus leading to overactivated BMP signaling that promotes bone formation and contributes to the SBD phenotype. Our results establish a previously unreported SBD entity (craniotubular dysplasia, Ikegawa type) and contribute to a better understanding of the regulation of BMP signaling and bone formation.

Impaired tumor immune response in metastatic tumors is a selective pressure for neutral evolution in CRC cases.

A Darwinian evolutionary shift occurs early in the neutral evolution of advanced colorectal carcinoma (CRC), and copy number aberrations (CNA) are essential in the transition from adenoma to carcinoma. In light of this primary evolution, we investigated the evolutionary principles of the genome that foster postoperative recurrence of CRC. CNA and neoantigens (NAG) were compared between early primary tumors with recurrence (CRCR) and early primary tumors without recurrence (precancerous and early; PCRC). We compared CNA, single nucleotide variance (SNV), RNA sequences, and T-cell receptor (TCR) repertoire between 9 primary and 10 metastatic sites from 10 CRCR cases. We found that NAG in primary sites were fewer in CRCR than in PCRC, while the arm level CNA were significantly higher in primary sites in CRCR than in PCRC. Further, a comparison of genomic aberrations of primary and metastatic conditions revealed no significant differences in CNA. The driver mutations in recurrence were the trunk of the evolutionary phylogenic tree from primary sites to recurrence sites. Notably, PD-1 and TIM3, T cell exhaustion-related molecules of the tumor immune response, were abundantly expressed in metastatic sites compared to primary sites along with the increased number of CD8 expressing cells. The postoperative recurrence-free survival period was only significantly associated with the NAG levels and TCR repertoire diversity in metastatic sites. Therefore, CNA with diminished NAG and diverse TCR repertoire in pre-metastatic sites may determine postoperative recurrence of CRC.

Neoantimon: a multifunctional R package for identification of tumor-specific neoantigens.

SUMMARY: It is known that some mutant peptides, such as those resulting from missense mutations and frameshift insertions, can bind to the major histocompatibility complex and be presented to antitumor T cells on the surface of a tumor cell. These peptides are termed neoantigen, and it is important to understand this process for cancer immunotherapy. Here, we introduce an R package termed Neoantimon that can predict a list of potential neoantigens from a variety of mutations, which include not only somatic point mutations but insertions, deletions and structural variants. Beyond the existing applications, Neoantimon is capable of attaching and reflecting several additional information, e.g. wild-type binding capability, allele specific RNA expression levels, single nucleotide polymorphism information and combinations of mutations to filter out infeasible peptides as neoantigen. AVAILABILITY AND IMPLEMENTATION: The R package is available at http://github/hase62/Neoantimon.

Association of single nucleotide polymorphisms in the NRF2 promoter with vascular stiffness with aging.

PURPOSE: Pulse wave velocity (PWV), an indicator of vascular stiffness, increases with age and is increasingly recognized as an independent risk factor for cardiovascular disease (CVD). Although many mechanical and chemical factors underlie the stiffness of the elastic artery, genetic risk factors related to age-dependent increases in PWV in apparently healthy people are largely unknown. The transcription factor nuclear factor E2 (NF-E2)-related factor 2 (Nrf2), which is activated by unidirectional vascular pulsatile shear stress or oxidative stress, regulates vascular redox homeostasis. Previous reports have shown that a SNP in the NRF2 gene regulatory region (-617C>A; hereafter called SNP-617) affects NRF2 gene expression such that the minor A allele confers lower gene expression compared to the C allele, and it is associated with various diseases, including CVD. We aimed to investigate whether SNP-617 affects vascular stiffness with aging in apparently healthy people. METHODS: Analyzing wide-ranging data obtained from a public health survey performed in Japan, we evaluated whether SNP-617 affected brachial-ankle PWV (baPWV) in never-smoking healthy subjects (n = 642). We also evaluated the effects of SNP-617 on other cardiovascular and blood test measurements. RESULTS: We have shown that not only AA carriers (n = 55) but also CA carriers (n = 247) show arterial stiffness compared to CC carriers (n = 340). Furthermore, SNP-617 also affected blood pressure indexes such as systolic blood pressure and mean arterial pressure but not the ankle brachial pressure index, an indicator of atherosclerosis. Multivariate analysis showed that SNP-617 accelerates the incremental ratio of baPWV with age. CONCLUSIONS: This study is the first to show that SNP-617 affects the age-dependent increase in vascular stiffness. Our results indicate that low NRF2 activity induces premature vascular aging and could be targeted for the prevention of cardiovascular diseases associated with aging.

A unified simulation model for understanding the diversity of cancer evolution.

Because cancer evolution underlies the therapeutic difficulties of cancer, it is clinically important to understand the evolutionary dynamics of cancer. Thus far, a number of evolutionary processes have been proposed to be working in cancer evolution. However, there exists no simulation model that can describe the different evolutionary processes in a unified manner. In this study, we constructed a unified simulation model for describing the different evolutionary processes and performed sensitivity analysis on the model to determine the conditions in which cancer growth is driven by each of the different evolutionary processes. Our sensitivity analysis has successfully provided a series of novel insights into the evolutionary dynamics of cancer. For example, we found that, while a high neutral mutation rate shapes neutral intratumor heterogeneity (ITH) characterized by a fractal-like pattern, a stem cell hierarchy can also contribute to shaping neutral ITH by apparently increasing the mutation rate. Although It has been reported that the evolutionary principle shaping ITH shifts from selection to accumulation of neutral mutations during colorectal tumorigenesis, our simulation revealed the possibility that this evolutionary shift is triggered by drastic evolutionary events that occur in a short time and confer a marked fitness increase on one or a few cells. This result helps us understand that each process works not separately but simultaneously and continuously as a series of phases of cancer evolution. Collectively, this study serves as a basis to understand in greater depth the diversity of cancer evolution.

The relationship between cigarette smoking and the tongue microbiome in an East Asian population.

Background: The oral microbiome, which consists of various habitats, has been shown to be influenced by smoking. However, differences in the tongue microbiomes of current and former smokers, as well as their resultant functional consequences, have rarely been investigated in East Asian populations. Methods: We used 16S rRNA amplicon sequencing of tongue-coating samples obtained from East Asian subjects who were current, former, or never smokers to identify differences in their tongue microbiomes and related metagenomic functions. Two sets of participants from 2016 to 2017 (n = 657 and n = 187, respectively) were analyzed separately. Results: We found significant differences between the overall microbiome compositions of current versus never smokers (p = 0.0015), but not between former versus never smokers (p = 0.43) based on the weighted UniFrac distance. Twenty-nine of 43 investigated genera showed significantly different expression levels in current versus never smokers. Neisseria and Capnocytophaga were less abundant, and Streptococcus and Megasphaera were more abundant in current smokers. Moreover, the abundances of metagenomic pathways, including those related to nitrate reduction and the tricarboxylic acid cycle, were significantly different between current and never smokers. Conclusions: The tongue microbiomes and related metagenomic pathways of current smokers differ from those of never smokers among East Asians.

Comprehensive analysis of indels in whole-genome microsatellite regions and microsatellite instability across 21 cancer types.

Microsatellites are repeats of 1- to 6-bp units, and approximately 10 million microsatellites have been identified across the human genome. Microsatellites are vulnerable to DNA mismatch errors and have thus been used to detect cancers with mismatch repair deficiency. To reveal the mutational landscape of microsatellite repeat regions at the genome level, we analyzed approximately 20.1 billion microsatellites in 2717 whole genomes of pan-cancer samples across 21 tissue types. First, we developed a new insertion and deletion caller (MIMcall) that takes into consideration the error patterns of different types of microsatellites. Among the 2717 pan-cancer samples, our analysis identified 31 samples, including colorectal, uterus, and stomach cancers, with a higher proportion of mutated microsatellite (≥0.03), which we defined as microsatellite instability (MSI) cancers of genome-wide level. Next, we found 20 highly mutated microsatellites that can be used to detect MSI cancers with high sensitivity. Third, we found that replication timing and DNA shape were significantly associated with mutation rates of microsatellites. Last, analysis of mutations in mismatch repair genes showed that somatic SNVs and short indels had larger functional impacts than germline mutations and structural variations. Our analysis provides a comprehensive picture of mutations in the microsatellite regions and reveals possible causes of mutations, as well as provides a useful marker set for MSI detection.

Metagenomic analysis of bacterial species in tongue microbiome of current and never smokers.

Cigarette smoking affects the oral microbiome, which is related to various systemic diseases. While studies that investigated the relationship between smoking and the oral microbiome by 16S rRNA amplicon sequencing have been performed, investigations involving metagenomic sequences are rare. We investigated the bacterial species composition in the tongue microbiome, as well as single-nucleotide variant (SNV) profiles and gene content of these species, in never and current smokers by utilizing metagenomic sequences. Among 234 never smokers and 52 current smokers, beta diversity, as assessed by weighted UniFrac measure, differed between never and current smokers (pseudo-F = 8.44, R2 = 0.028, p = 0.001). Among the 26 species that had sufficient coverage, the SNV profiles of Actinomyces graevenitzii, Megasphaera micronuciformis, Rothia mucilaginosa, Veillonella dispar, and one Veillonella sp. were significantly different between never and current smokers. Analysis of gene and pathway content revealed that genes related to the lipopolysaccharide biosynthesis pathway in Veillonella dispar were present more frequently in current smokers. We found that species-level tongue microbiome differed between never and current smokers, and 5 species from never and current smokers likely harbor different strains, as suggested by the difference in SNV frequency.

Classification of primary liver cancer with immunosuppression mechanisms and correlation with genomic alterations.

BACKGROUND: The tumor microenvironment can be classified into immunologically active "inflamed" tumors and inactive "non-inflamed" tumors based on the infiltration of cytotoxic immune cells. Previous studies on liver cancer have reported a superior prognosis for inflamed tumors compared to non-inflamed tumors. However, liver cancer is highly heterogeneous immunologically and genetically, and a finer classification of the liver cancer microenvironment may improve our understanding of its immunological diversity and response to immune therapy. METHODS: We characterized the immune gene signatures of 234 primary liver cancers, mainly virus-related, from a Japanese population using RNA-Seq of tumors and matched non-tumorous hepatitis livers. We then compared them with the somatic alterations detected using the whole-genome sequencing. FINDINGS: Liver cancers expressed lower levels of immune marker genes than non-tumorous hepatitis livers, indicating immunosuppression in the tumor microenvironment. Several immunosuppression mechanisms functioned actively and mutually exclusively, resulting in four immune subclasses of liver cancer: tumor-associated macrophage (TAM), CTNNB1, cytolytic activity (CYT), and regulatory T cell (Treg). The CYT and Treg subclasses represented inflamed tumors, while the TAM and CTNNB1 subclasses represented non-inflamed tumors. The TAM subclass, which comprised 31% of liver cancers, showed a poor survival, expressed elevated levels of extracellular matrix genes, and was associated with somatic mutations of chromatin regulator ARID2. The results of cell line experiments suggested a functional link between ARID2 and chemokine production by liver cancer cells. INTERPRETATION: Primary liver cancer was classified into four subclasses based on mutually exclusive mechanisms for immunosuppression. This classification indicate the importance of immunosuppression mechanisms, such as TAM and Treg, as therapeutic targets for liver cancer. FUNDING: The Japan Agency for Medical Research and Development (AMED).

Contribution of Rare Variants of the SLC22A12 Gene to the Missing Heritability of Serum Urate Levels.

Gout is a common arthritis caused by monosodium urate crystals. The heritability of serum urate levels is estimated to be 30-70%; however, common genetic variants account for only 7.9% of the variance in serum urate levels. This discrepancy is an example of "missing heritability." The "missing heritability" suggests that variants associated with uric acid levels are yet to be found. By using genomic sequences of the ToMMo cohort, we identified rare variants of the SLC22A12 gene that affect the urate transport activity of URAT1. URAT1 is a transporter protein encoded by the SLC22A12 gene. We grouped the participants with variants affecting urate uptake by URAT1 and analyzed the variance of serum urate levels. The results showed that the heritability explained by the SLC22A12 variants of men and women exceeds 10%, suggesting that rare variants underlie a substantial portion of the "missing heritability" of serum urate levels.

Carnosine induces intestinal cells to secrete exosomes that activate neuronal cells.

Recently, we showed that imidazole dipeptide such as carnosine contained abundantly in chicken breast meat improves brain function in a double-blind randomized controlled trial. However, the underlying molecular mechanisms remain unknown. Here, we investigated whether carnosine activates intestinal epithelial cells and induces the secretion of factors that activate brain function. We focused on exosomes derived from intestinal epithelial cells as mediators of brain-gut interaction. Results showed that exosomes derived from Caco-2 cells treated with carnosine significantly induced neurite growth in SH-SY5Y cells. To clarify the molecular basis of this finding, we performed integrated analysis of microRNAs (miRNAs) with altered expression in exosomes in response to carnosine treatment and mRNAs with altered expression in target cells in response to exosome treatment to identify related miRNAs and their target genes. The combination of miR-6769-5p and its target gene ATXN1 was found to be involved in the exosome-induced activation of neuronal cells.

Human NK cell development in hIL-7 and hIL-15 knockin NOD/SCID/IL2rgKO mice.

The immune system encompasses acquired and innate immunity that matures through interaction with microenvironmental components. Cytokines serve as environmental factors that foster functional maturation of immune cells. Although NOD/SCID/IL2rgKO (NSG) humanized mice support investigation of human immunity in vivo, a species barrier between human immune cells and the mouse microenvironment limits human acquired as well as innate immune function. To study the roles of human cytokines in human acquired and innate immune cell development, we created NSG mice expressing hIL-7 and hIL-15. Although hIL-7 alone was not sufficient for supporting human NK cell development in vivo, increased frequencies of human NK cells were confirmed in multiple organs of hIL-7 and hIL-15 double knockin (hIL-7xhIL-15 KI) NSG mice engrafted with human hematopoietic stem cells. hIL-7xhIL-15 KI NSG humanized mice provide a valuable in vivo model to investigate development and function of human NK cells.

Sensitivity analysis of agent-based simulation utilizing massively parallel computation and interactive data visualization.

An essential step in the analysis of agent-based simulation is sensitivity analysis, which namely examines the dependency of parameter values on simulation results. Although a number of approaches have been proposed for sensitivity analysis, they still have limitations in exhaustivity and interpretability. In this study, we propose a novel methodology for sensitivity analysis of agent-based simulation, MASSIVE (Massively parallel Agent-based Simulations and Subsequent Interactive Visualization-based Exploration). MASSIVE takes a unique paradigm, which is completely different from those of sensitivity analysis methods developed so far, By combining massively parallel computation and interactive data visualization, MASSIVE enables us to inspect a broad parameter space intuitively. We demonstrated the utility of MASSIVE by its application to cancer evolution simulation, which successfully identified conditions that generate heterogeneous tumors. We believe that our approach would be a de facto standard for sensitivity analysis of agent-based simulation in an era of evergrowing computational technology. All the results form our MASSIVE analysis are available at https://www.hgc.jp/~niiyan/massive.