dataesr/scientific-paragraphs-categorization huggingface dataset:
We then conducted differential expression analysis on the clusters and identified key marker genes for each cell type (Supplemental Table  1 ). To compare the consistency between cluster assignment in the snATAC-seq data and the scRNAseq data, we next derived a gene activity score for the top 3,000 highly variable genes in each was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which this version posted June 4, 2020. ; https://doi. org/10. 1101/2020. 05. 24. 113910 doi: bioRxiv preprint 7/50 snATAC-seq cluster and computed the Pearson's correlation coefficient between each snATAC cluster and scRNA cluster (Methods). This analysis indicated good concordance between the two datasets (Figures  1g, S1p ).
While the correlation between gene expression and inferred gene activity score was high, we noted some differences in cell proportions, which was mostly related to the sample preparation-induced cell drop-out (Figures  S1g, i ). Consistent with previous observations that single cell preparations better capture immune cells than single nuclear preparations  21  , we noted that the immune cell repertoire was limited in the snATAC-seq dataset; on the other hand, stromal cells were better captured by the nuclear preparation. Finally, to allow the interactive use of this dataset by the community, we not only made the raw data available but also the processed dataset via our searchable website (susztaklab. com/developing_adult_kidney/snATAC/ for snATAC-seq data, susztaklab. com/developing_adult_kidney/scRNA/ for scRNA-seq data, and susztaklab.
We reasoned that a subset of the cell type-specific cis-regulatory elements should regulate cell type-specific gene expression in cis. Inspired by Zhu et al. 23  , we aligned DAPs and differentially expressed genes from our snATAC-seq and scRNA-seq datasets, and inferred the putative regulatory peaks by their proximity (Methods). Such cis-regulatory elements predictions were confirmed by comparing with cis-regulatory elements inferred previously  24  , as we recapitulated roughly 20% of elements from their analysis. In addition, our analysis was able to identify several known enhancers such as for Six2 and Slc6a18  24, 25  (Figure  S2g ).
11/50 The regulatory trajectory of nephron progenitor differentiation All cells in the body differentiate from the same genetic template. Cell type-specific chromatin opening and closing events associated with TF binding changes set up the cell type-specific regulatory landscape resulting in cell type specification and development. We found that closing of open chromatin regions was the predominant event during the nephron progenitor differentiation (Figure  S2a ). We then evaluated the cellular differentiation trajectory in the snATAC-seq and scRNA-seq datasets (Methods). We identified multiple nephron progenitor subgroups (Figures  3a-b ), which will need to be carefully mapped to prior gene expression-and anatomical location-driven nephron progenitor sub-classification. Consistently, across both data modalities, we identified that the podocyte precursors differentiated early from the nephron progenitor pool (Figures  3a-b ).
The tubule cell trajectory was more complex with a shared intermediate stage and later differentiation into proximal tubules and distal tubules/loop of Henle (Figures  3a-b, S3a-c ). We also integrated snATAC-seq and scRNA-seq data to obtain a single trajectory (Methods). The cell types in this dataset were correctly mapped and the trajectory resembled the path observed in individual analyses of the scRNA and snATAC datasets (Figures  S3e-g ). The robustness of developmental trajectories was further supported by obtaining similar results when performing RNA velocity analysis using Velocyto 30 (Figure  S3d ) and by comparing with previous human and mouse kidney developmental studies  9, 13, 14. Building on both the SCENIC-generated gene regulatory network and the robust differentiation trajectories of the snATAC-seq and scRNA-seq datasets, we next aimed to understand chromatin dynamics, identify TFs and driver pathways for cell type specification and differentiation.
Discussion In summary, here we present the first cellular resolution open chromatin map for the developing and adult mouse kidney. Using this dataset, we identified key cell type-specific regulatory networks for kidney cells, defined the cellular differentiation trajectory, characterized regulatory dynamics and identified key driving TFs for nephron development, especially for the terminal differentiation of epithelial cells. Furthermore, our results shed light on the cell types and target genes for genetic variants associated with kidney disease development. By performing massively parallel single cell profiling of chromatin state, we were able to define the key regulatory logic for each kidney cell type by investigating cis-regulatory elements and TFtarget gene interaction. We found that most cell type-specific open chromatin regions are within distal regulatory elements and intronic regions.
In addition, one needs to be aware of the limitations when interpreting different computational analyses, for example, the motif enrichment analyses such as implemented by HOMER, SCENIC, and chromVAR, are not able to distinguish between TFs with similar binding sites. Future highthroughput studies that analyze open chromatin and gene expression information from the same was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which this version posted June 4, 2020. ; https://doi. org/10. 1101/2020. 05. 24. 113910 doi: bioRxiv preprint 22/50 cells will be exceedingly helpful to correlate open chromatin and gene expression information along the differentiation trajectory  24, 60, 61. In summary, our dataset provides critical novel insight into the cell type-specific gene regulatory network, cell differentiation program, and disease development.
Dimension reduction, batch effect correction and clustering Clustering was conducted using snapATAC  16  , a single-cell ATAC-seq algorithm scalable to large dataset. Previous benchmarking evaluation has shown that snapATAC was one of the bestperforming methods for snATAC-seq clustering  65. Diffusion map was applied as a dimension reduction method using function runDiffusionMaps. To remove batch effect, we used Harmony 17 , in which the low dimensional embeddings obtained from the diffusion map were used as input. Harmony iteratively pulled batch-specific centroid to cluster centroid until convergence to remove the variability across batches. After batch correction, a graph was constructed using k Nearest Neighbor (kNN) algorithm with k=15, which was then used as input for Louvain clustering.
Peak calling and visualization Peak calling was conducted for each cell type separately using MACS2  18. We aggregated all fragments obtained from the same cell types to build a pseudo-bulk ATAC data and conducted peak calling with parameters "--nomodel --keep-dup all --shift 100 --ext 200 --qval 1e-2 -B --SPMR --call-summits". By specifying "--SPMR", MACS2 generated "fragment pileup per million reads" pileup files, which were converted to bigwig format for visualization using UCSC bedGraphToBigWig tool. We also visualized public chromatin ChIP-seq data and RNA-seq data obtained from ENCODE Encyclopedia with the following identifiers: ENCFF338WZP, ENCFF872MVE, ENCFF455HPY, ENCFF049LRQ, ENCFF179NTO, ENCFF071PID, ENCFF746MFH, ENCFF563LOO, ENCFF184AYF, ENCFF107NQP, ENCFF465THI, ENCFF769XWI, ENCFF591DAX.
Genomic elements stratification Mouse mm10 genome annotation files were download from UCSC Table  Browser  (https://genome. ucsc. edu/cgi-bin/hgTables) using GENCODE VM23. TSS upstream 5 kb regions were included as promoter regions, but the results were similar when using 2 kb upstream regions as promoters. We then studied the number of overlapped regions between open chromatin regions identified from the snATAC-seq and bulk ATAC-seq dataset and genome annotations. Since one open chromatin region could overlap with multiple genomic elements, we defined an order of genomic elements as exon > 5'-UTR > 3'-UTR > intron > promoter > distal elements. To be more specific, if one peak overlapped with both exon and 5'-UTR, the algorithm would count it as an exon-region peak.
Bulk ATAC sequencing analysis Bulk ATAC-seq raw fastq files were processed using the end-to-end tool ENCODE ATAC-seq pipeline (Software and Algorithms). This tool provided a standard workflow for ATAC-seq data was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which this version posted June 4, 2020. ; https://doi. org/10. 1101/2020. 05. 24. 113910 doi: bioRxiv preprint 33/50 quality control, adaptor removal, alignment, and peak calling. To obtain high quality ATAC-seq peaks, peak calling results from two biological replicates were compared and only those peaks that were present in both replicates were kept, which were further used to compare with snATAC-seq peaks.
Correlation of bulk and single nuclei ATAC sequencing data snATAC-seq reads were aggregated to a pseudo-bulk data for the comparison purpose. To prevent the effect of sex chromosome and mitochondria chromosome, reads from chromosome X, Y and M were excluded from our analysis. We used multiBigwigSummary tool from deeptools  69  to study the correlation between different samples. Specifically, the whole genome was binned into equally sized (10 kb) windows, and the reads in each bin were aggregated, generating a bin-read count vector for each of the sample. The correlation of these vectors was computed as a measure of pairwise similarity between samples. To compare the number of peaks in these two datasets, we used as input the narrowpeak files from the snATAC-seq and bulk ATAC-seq analysis. We filtered out bulk ATAC-seq peaks with q value > 0.
Comparison between single nuclei ATAC sequencing data and single cell RNA sequencing data In order to compare the cluster assignment between snATAC-seq data and scRNA-seq data, we obtained the average gene expression values and peak accessibility in each cluster for P0 and adult samples separately. We next transformed snATAC-seq data by summing up the reads within gene body and 2 kb upstream regions to build gene activity score matrix, as suggested in Seurat  20. Then, we normalized the data and computed the mean expression and mean gene activity scores in each cell type, and calculated z scores of each gene. Pearson's correlation coefficient was then calculated among top 3,000 highly variable genes between snATAC-seq data and scRNA-seq data. was not certified by peer review) is the author/funder.
Alignment and quality control Raw fastq files were aligned to the mm10 (Ensembl GRCm38. 93) reference genome and quantified using CellRanger v3. 1. 0. Seurat R package v3. 0 was used for data quality control, preprocessing and dimensional reduction analysis. After gene-cell data matrix generation of both P0 and adult datasets, matrices were merged and poor-quality cells with <200 or >3,000 expressed genes and mitochondrial gene percentages >50 were excluded, leaving 25,138 P0 and 18,498 adult cells for further analytical processing, respectively (Figures  S1j-k ).
Cells were re-clustered by Monocle3 using a resolution of 0. 0005 with k-nearest neighbor (kNN) k=29. Highly variable genes along pseudotime were identified using differential GeneTest function and cells were ordered along pseudotime trajectory. NP cluster was defined as earliest principal node. In order to find genes differentially expressed along pseudotime, trajectories for podocytes, loop of Henle, and proximal tubule clusters were analyzed separately with the fit_models function of Monocle3. Genes with a q value <0. 05 in the differential GeneTest analysis were kept. In an alternate approach, graph_test function of Monocle3 was used and trajectory-variable genes were collected into modules at a resolution of 0. 01.
RNA velocity To calculate RNA velocity, Python-based Velocyto command-line tool as well as Velocyto. R package were used as instructed  30. We used Velocyto to calculate the single-cell trajectory/directionality using spliced and unspliced reads. From loom files produced by the command-line tool, we subset the exact same cells that were previously selected randomly for Monocle trajectory analysis. This subset was loaded into R using the SeuratWrappers v0. 1. 0 package. RNA velocity was estimated using gene-relative model with k-nearest neighbor (kNN) cell pooling (k = 25). The parameter n was set at 200, when visualizing RNA velocity on the UMAP embedding.
Gene regulatory network inference In order to identify TFs and characterize cell states, we employed cis-regulatory analysis using the R package SCENIC v1. 1. 2. 2  28  , which infers the gene regulatory network based on co-expression and DNA motif analysis. The network activity is then analyzed in each cell to identify recurrent cellular states. In short, TFs were identified using GENIE3 and compiled into modules (regulons), which were subsequently subjected to cis-regulatory motif analysis using RcisTarget with two gene-motif rankings: 10 kb around the TSS and 500 bp upstream. Regulon activity in every cell was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which this version posted June 4, 2020.
Ligand-receptor interactions To assess cellular crosstalk between different cell types, we used the CellPhoneDB repository to infer cell-cell communication networks from single cell transcriptome data  40. We used the Python package CellPhoneDB v2. 1. 2 together with the database v2. 0. 0 to predict cell type-specific ligandreceptor complexes as per the authors' instructions. Only receptors and ligands expressed in more than 5% of the cells in the specific cluster were considered. 1,000 iterations were used for pairwise comparison between cell types and considered for further statistical analysis.
D) Genome browser view of read density in each snATAC-seq cluster at cell type marker gene transcription start sites. We used Uncx for nephron progenitors, Nphs2 for podocytes, Akr1c21 for proximal tubules, Slc12a1 for loop of Henle, Slc12a3 for distal convoluted tubule, Trpv5 for connecting tubule, Aqp2 for collecting duct principal cells, Atp6v1g3 for intercalated cells, Egfl7 for endothelial cells, C1qb for immune cells and Col3a1 for stroma. Additional marker gene examples are shown in Figure S1h. (E) Comparison of peaks identified from snATAC-seq data and bulk ATAC-seq data. Peaks that are identified in both datasets are colored blue, and peaks that are dataset-specific are grey. (F) Violin plots showing cell type-specific gene expression in scRNA-seq data. With the exception of proximal tubule, the same marker genes as in snATAC-seq data were used (Slc5a2 and Slc22a30 for proximal tubule S1 and S3, respectively).
(G) Correlation between snATAC-seq gene activity scores and gene expression values in P0 data. The correlation of the adult dataset is shown in FigureS1p. was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. The copyright holder for this preprint (which this version posted June 4, 2020. ; https://doi. org/10. 1101/2020. 05. 24.
; https://doi. org/10. 1101/2020. 05. 24. Distribution of number of unique molecular identifiers (UMIs, x axis) and promoter ratio (y axis) in 5 samples shown by dot plot. (E) Violin plots representing the number of accessible peaks across different clusters in the snATAC-seq dataset indicating similar distributions. (F) UMAP representation of the snATAC-seq dataset colored by batches. (G) Stacked bar graphs representing absolute numbers and percentages of identified cell types across snATAC-seq batches. (H) Genome browser view of cell type-specific peaks at the TSS of marker genes for 13 cell types in the snATAC-seq dataset. (I) From left to right: Stacked bar graphs showing the percentage of different cell types in the P0 and adult scRNA-seq datasets, tables showing the number of cells in each cell type (nCells) and corresponding percentage.
NP, nephron progenitor; Podo, podocyte; PT, proximal tubule; S1, segment 1; S3, segment 3; LOH, loop of Henle; DCT, distal convoluted tubule cells; PC, collecting duct principal cells; IC, collecting duct intercalated cells; Endo, endothelial cells; Macro, macrophages; Neutro, neutrophils. (J) UMAP representation of scRNA-seq data colored by the mitochondrial gene ratio (Mt %). (K) Violin plots showing number of informative genes per single cell and unique molecular identifiers (UMIs) per single cell. Blue denotes adult kidney, orange denotes P0 kidney. (L, M) Principal component (PC) representation of combined adult and P0 scRNA-seq dataset (left panel) and violin plots of corresponding embeddings values (right panel) before (L) and after (M) batch correction using Harmony. (N) Dot plot of cell type-specific marker genes.
Dot size denotes percentage of cells expressing the marker. Color scale represents average expression, orange denotes P0, blue denotes adult kidney. (O) Feature plots of representative marker genes projected on UMAP dimension. (P) Correlation between snATAC-seq gene activity scores and gene expression values in adult data, which is complementary to Figure 1g. (Q) We provide the processed chromatin accessibility dataset via a searchable, interactive website (susztaklab. com/developing_adult_kidney/igv/). Ace2 was used as an example, and we show proximal tubule-specific enrichment of peaks at transcription start sites of the Ace2 (Angiotensinconverting enzyme 2) gene (red boxes).
data files and data matrix are being uploaded onto GEO and an accession number will be provided when it becomes available. The annotated and analyzed data can be viewed at susztaklab. com/developing_adult_kidney/snATAC/, susztaklab. com/developing_adult_kidney/scRNA/, and susztaklab. com/developing_adult_kidney/igv/. Further information and requests for resources and was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission.
Cell Ranger ATAC v. 1. 1. 0 10X Genomics https://support. 10xgenomics. com/sin gle-cellatac/software/downloads/latest Cell Ranger v. 3. 1. 0 10X Genomics https://support. 10xgenomics. com/sin gle-cell-geneexpression/software/downloads/latest CellPhoneDB v. 2. 1. 2 open source https://www. cellphonedb. org ChromVAR v. 3. 1. 0 open source http://bioconductor. org/packages/rele ase/bioc/html/chromVAR. html Cicero v. 1. 5. 5 open source https://github. com/cole-trapnelllab/cicero-release deeptools v. 2. 0 open source https://deeptools. readthedocs. io/en/d evelop/ ENCODE ATAC-seq pipeline open source https://github. com/kundajelab/atac_d nase_pipelines HOMER v. 4. 10. 4 open source http://homer. ucsd. edu/homer/motif/ IGV open source http://software. broadinstitute. org/soft ware/igv/ MACS2 v.
2. 2. 6 open source https://github. com/taoliu/MACS Monocle3 v. 0. 1. 3 open source http://cole-trapnelllab. github. io/monocle-release/ Harmony open source https://github. com/immunogenomics /harmony SCENIC v. 1. 1. 2. 2 open source https://aertslab. org/#scenic Seurat R package v. 3. 0 open source https://satijalab. org/seurat/ snapATAC open source https://github. com/r3fang/SnapATA C UCSC bedgraphtobigwig open source http://hgdownload. soe. ucsc. edu/admi n/exe/linux. x86_64. v385/ UCSC liftOver open source https://genome. ucsc. edu/cgibin/hgLiftOver VelocytoR open source https://github. com/velocytoteam/velocyto. R VisCello open source https://github. com/qinzhu/VisCello.
( D) Genome browser view of read density in each snATAC-seq cluster at cell type marker gene transcription start sites. We used Uncx for nephron progenitors, Nphs2 for podocytes, Akr1c21 for proximal tubules, Slc12a1 for loop of Henle, Slc12a3 for distal convoluted tubule, Trpv5 for connecting tubule, Aqp2 for collecting duct principal cells, Atp6v1g3 for intercalated cells, Egfl7 for endothelial cells, C1qb for immune cells and Col3a1 for stroma. Additional marker gene examples are shown in Figure S1h. (E) Comparison of peaks identified from snATAC-seq data and bulk ATAC-seq data. Peaks that are identified in both datasets are colored blue, and peaks that are dataset-specific are grey. (F) Violin plots showing cell type-specific gene expression in scRNA-seq data. With the exception of proximal tubule, the same marker genes as in snATAC-seq data were used (Slc5a2 and Slc22a30 for proximal tubule S1 and S3, respectively).
67) GSM1051156 Mouse kidney H3K27ac and ENCODE project ENCFF338WZP, ENCFF872MVE, or RESOURCE H3K4me1 CHIP-seq, WGBS, and RNA-seq SOURCE ENCFF455HPY, ENCFF049LRQ, IDENTIFIER ENCFF179NTO, ENCFF071PID, ENCFF746MFH, ENCFF563LOO, Chemicals, Antibodies, Peptides, and Recombinant Proteins Guinea pig anti-Foxl1 own production ENCFF184AYF, ENCFF107NQP, (ref. 73) Mouse anti-E-cadherin BD Transducton ENCFF465THI, ENCFF769XWI, Cat#610182 Cy2-conjugated donkey secondary antibody Jackson ENCFF591DAX Cat#715-225-150 Six2 ChIP-seq data in nephron (ref. 66) GUDMAP database (RID:Q-Y4CY) ImmunoResearch progenitor cells Laboratories Software and Algorithms Cy3-conjugated donkey bedtools v. 2. 29. 2 Jackson open source Cat#715-165-150 https://bedtools. readthedocs. io/en/lat secondary antibody ImmunoResearch est/ Laboratories 37/50 38/50 was.
Lead Contact and Data Availability Raw data files and data matrix are being uploaded onto GEO and an accession number will be provided when it becomes available. The annotated and analyzed data can be viewed at susztaklab. com/developing_adult_kidney/snATAC/, susztaklab. com/developing_adult_kidney/scRNA/, and susztaklab. com/developing_adult_kidney/igv/. Further information and requests for resources and was not certified by peer review) is the author/funder. All rights reserved. No reuse allowed without permission. bioRxiv preprint doi: ; https://doi. org/10. 1101/2020. 05. 24. 113910 this version posted June 4, 2020. The copyright holder for this preprint (which 39/50.
Hi-C exploits contact frequencies between pairs of loci to bridge and order contigs during genome assembly, resulting in chromosome-level assemblies. Because few robust programs are available for this type of data, we developed instaGRAAL, a complete overhaul of the GRAAL program, which has adapted the latter to allow efficient assembly of large genomes. instaGRAAL features a number of improvements over GRAAL, including a modular correction approach that optionally integrates independent data. We validate the program using data for two brown algae, and human, to generate near-complete assemblies with minimal human intervention.
The latter procedure derives from techniques aiming at recovering snapshots of the higherorder organization of a genome  [8, 9]. When applied to genomics, Hi-C-based methods are sometimes referred to as proximity ligation approaches, as they quantify and exploit physical contacts between pairs of DNA segments in a genome to assess their collinearity along a chromosome, and the distance between the segments  [10]. Early studies using control datasets demonstrated that Hi-C can be used to scaffold and/or correct a wide range of eukaryotic DNA regions  [11] [12] [13] [14] , i. e. stretches of bp, whether they be small-scale contigs or full chromosomes. The Hi-C scaffolder GRAAL (Genome Re-Assembly Assessing Likelihood from 3D) is a probabilistic program that uses a Markov Chain Monte Carlo (MCMC) approach.