A multitask clustering approach for single-cell RNA-seq analysis in Recessive Dystrophic Epidermolysis Bullosa
[['In recent years , single-cell RNA sequencing ( scRNA-seq ) has emerged as the dominant method for quantifying transcriptome-wide mRNA expression in individual cells .', 'While traditional bulk RNA-seq ignores the differences between individual cells and treats the population of cells as homogeneous , scRNA-seq identifies sub-populations of single cells and can be useful for characterizing sub-population structure , mechanisms of transcription regulation , and understanding disease progression [1] and immunology [2] .', 'A typical scRNA-seq protocol consists of several steps: isolation of single cells and RNA , reverse transcription , amplification , library generation , and sequencing .', 'In addition to the noise and bias that exist in bulk RNA-seq experiments , issues unique to scRNA-seq include those from biological sources , such as cell-cycle stage or cell size , as well as from technical/systematic sources , such as capture inefficiency , material degradation , sample contamination , amplification biases , GC content , and sequencing depth .', 'These experimental biases and limitations cause uneven coverage of the entire transcriptome and result in an abundance of zero-coverage regions [3 , 4] .', 'Typically , the cost of scRNA-Seq is much higher than bulk RNA-Seq per sample , and thus , scRNA-Seq of a large patient cohort is still prohibitively costly .', 'When a large number of single-cells from multiple samples are sequenced , more complex batch effects might be introduced .', 'Finally , some poorly sampled cell populations might only contain very few cells for the analysis .', 'To address all these challenges , proper integration of multiple scRNA-Seq datasets generated from different experiments is important .', 'When multiple single-cell populations from biological replicates or related samples such as a patient cohort are analyzed to discover the common and sample-specific cell types , technical biases and irrelevant biological variance among independent samples cannot be easily identified and removed from the signal before clustering the single cells .', 'For example , when the scRNA-seq profiles from multiple patients are pooled together for clustering , the clusters will highly overlap with the division of the single cells by the sample origins rather than similar types such as pathogenic cells vs normal cells .', 'In this paper , we introduce a multitask learning method with embedded feature selection to simultaneously capture the differentially expressed genes among cell clusters and across all cell populations to achieve better single-cell clustering .', 'The key advantage of multitask clustering is the use of multiple single-cell populations to leverage the sample size limitation in each individual dataset while allowing dataset-specific variations among the same cell types across the datasets .', 'To illustrate the objective , Fig 1 shows a simulation example of scRNA-seq data of 100 single cells from three cell populations ( n = 33 , 33 and 34 ) with 1000 expressed genes .', 'Among the 1000 genes , gene A and gene B are the hidden markers that are differentially expressed across the four cell types ( indicated by four different colors ) .', 'In the ideal scenario , there is no technical bias and the marker genes are known as shown in the ground truth in Fig 1 ( A ) .', 'Fig 1 ( B ) shows the single-cell datasets after biological variation , technical biases , and noise are introduced .', 'The data distributions are very different across the three cell populations after the rotation , re-scaling and addition of noise .', 'It is also challenging to identify the true marker genes with a limited number of samples in each population .', 'Simply pooling the single-cell data from the three populations together will confuse the clustering , even with the correct marker genes identified ( Fig 1 ( C ) ) .', 'Conversely , separated clustering on each single-cell population suffers more from the biological variation as the number of single cells is not sufficient in each individual analysis to identify the true maker genes ( Fig 1 ( D ) ) .', 'As shown in Fig 1 ( E ) , variance-driven multitask clustering of single-cell RNA-seq data ( scVDMC ) utilizes expression patterns of different single-cell populations with shared cell-type markers and corresponding similar clusters for better integration .'], ['Assume a total of D domains with each domain representing a single-cell population for clustering .', 'Let matrix X ( d ) ∈ R m × n ( d ) denote RNA-seq gene expression values from domain d ∈ {1 , 2 , … , D} , where m is the number of features ( genes ) and n ( d ) is the single-cell sample size of domain d .', 'Let U ( d ) ∈ R m × k denote the cell-type cluster centers , vector Y i , j = [ U i , j ( 1 ) , U i , j ( 2 ) , … , U i , j ( D ) ] T stack the ( i , j ) -th entry of every U ( d ) and the binary matrix V ( d ) ∈ {0 , 1}n ( d ) ×k denote the assignments of each single-cell to the clusters , where k is the number of cell types ( clusters ) .', 'With the binary vector B ∈ {0 , 1}m denoting the indicators of feature selection ( 1: selected and 0: not selected ) and DB denoting the diagonal matrix with B on the diagonal , scVDMC model outlined in Fig 2 is defined as:', 'minimize U ( d ) , V ( d ) , B 1 2 ∑ d = 1 D | | D B ( X ( d ) - U ( d ) V ( d ) T ) | | F 2 - w ∑ d = 1 D B T Var ( U ( d ) ) + α ∑ i , j B i Var ( Y ( i , j ) ) subject to ∑ B = λ , ∑ j V i , j ( d ) = 1 , ∀ i = 1 , 2 , … , n ( d ) , ∀ d = 1 , 2 , … , D , ( 1 )', 'where w and α > 0 are hyper-parameters to balance the three error terms: the reconstruction error , the cluster center separation in each cell population , and the variance of the cluster centers across the different single-cell populations .', 'λ ∈ Z + is the predefined number of features to be selected .', '| | D B ( X ( d ) - U ( d ) V ( d ) T ) | | F 2 in Eq ( 1 ) denotes the reconstruction error of the classic k-means clustering as matrix factorization with DB selecting marker genes by B , i . e . the reconstruction error is only measured on the marker genes by ignoring the irrelevant ( non-selected ) genes .', 'The second term BTVar ( U ( d ) ) is introduced to maximize the separation of the cluster centers , where Var ( U ( d ) ) is defined as a vector in which each element is the variance of the vector U i , : ( d ) ∈ R k × 1 [5] .', 'The third term Var ( Y ( i , j ) ) denotes the variance of the vector Y ( i , j ) , which is introduced to require similar gene expression centers across different single-cell populations .', 'Note that the reconstruction error encourages selection of low expression genes since the errors are usually smaller on smaller values while the second variance term encourages selection of high expression genes since the variances tend to be larger on larger values .', 'Together as the sum over all the domains , the cost function provides a balanced error on the compactness and separation of the clusters of the cell types tuned by feature selection across all the domains .', 'The unique but similar cluster centers in each domain preserves the unique expression patterns while the features are selected as common marker genes for different cell types .', 'For the three hyper-parameters in Eq ( 1 ) , λ ( the number of marker genes ) is typically a small number based on prior knowledge of the cell types , and the selection of balancing weight w and α is discussed later in this section .', 'Algorithm 1 scVDMC algorithm 1: Input: X ( d ) , α , k , w , λ , d = 1 , 2 , … , D 2: output: U ( d ) , V ( d ) , B 3: Initialize U ( d ) and V ( d ) .', '4: repeat 5: compute B with integer linear programming in Eq ( 7 ) 6: for d = 1 , 2 , … , D do 7: solve V ( d ) by Eq ( 2 ) 8: solve U ( d ) by ( 6 ) 9: end for 10: until U ( d ) , V ( d ) and B converge 11: return U ( d ) , V ( d ) and B The full scVDMC algorithm is shown in Algorithm 1 .', 'The goal is to minimize the cost function in Eq ( 1 ) to obtain the optimal U ( d ) , V ( d ) and B . We employ an alternating update strategy to solve the optimization problem .', 'First , we fix the feature selection B , all the cluster centers U ( i ) , i = 1 , 2 , … , D and all other V ( i ) , i ≠ d , to obtain a certain V ( d ) .', 'minimize V ( d ) 1 2 | | D B ( X ( d ) - U ( d ) V ( d ) T ) | | F 2 subject to ∑ j V i , j ( d ) = 1 , ∀ i = 1 , 2 , … , n ( d ) .', '( 2 )', 'This is equivalent to assigning samples to the nearest centers U ( d ) by the Euclidean distance in the features selected by B , where each column of DB X ( d ) is a sample and each column of DB U ( d ) is a center .', 'Then the distance of a sample to every center is calculated and the nearest center is chosen to assign 1 to the corresponding V ( d ) .', 'The time complexity for assigning each sample to one of the k clusters over the λ marker genes will be O ( n × k × λ ) , where n is the total number of samples in all the domains .', 'Next , we fix the feature selection B , all clustering assignments V ( i ) , i = 1 , 2 , … , D , and all other U ( i ) , i ≠ d , to solve a certain U ( d ) , rewritten as:', 'minimize U ( d ) 1 2 ∑ i = 1 m B i | | ( X i , : ( d ) - U i , : ( d ) V ( d ) T ) | | 2 2 - w ∑ i = 1 m B i Var ( U i , : ( d ) ) + α ∑ i , j B i Var ( Y ( i , j ) ) , ( 3 )', 'where Var ( U i , : ( d ) ) is the variance of vector U i , : ( d ) defined as', 'Var ( U i , : ( d ) ) = 1 k ( U i , : ( d ) - U i , : ( d ) 1 k 1 k T k ) ( U i , : ( d ) - U i , : ( d ) 1 k 1 k T k ) T = 1 k U i , : ( d ) ( I k - 1 k 1 k T k ) ( I k - 1 k 1 k T k ) T U i , : ( d ) T = 1 k U i , : ( d ) ( I k - 1 k 1 k T k ) U i , : ( d ) T , ( 4 )', 'where Ik denotes the identity matrix of size k and 1k is a length k column vector of all ones .', 'Similarly , we have', 'Var ( Y ( i , j ) ) = 1 d Y ( i , j ) T ( I d - 1 d 1 d T d ) Y ( i , j ) .', '( 5 )', 'As shown in S1 Appendix , the analytical solution of Eq ( 3 ) when Bi = 1 is', 'Ui , : ( d ) T= ( V ( d ) TV ( d ) −2wkΨ+2αkΦd , ddIk ) −1 ( V ( d ) TXi , : ( d ) T−2αkd∑l≠dΦdlUi , : ( l ) T ) .', '( 6 )', 'The time complexity is O ( k3 ) for the matrix inversion and O ( n × k2 ) for computing V ( d ) T V ( d ) .', 'Since the matrix inversion is common to all the genes and only needs to be computed once , the total time complexity is only O ( n × k × λ ) .', 'Finally , to update binary vector B , we fix all U ( d ) and V ( d ) to optimize', 'minimize B ∑ i = 1 m B i ( 1 2 ∑ d = 1 D | | ( X i , : ( d ) - U i , : ( d ) V ( d ) T ) | | 2 2 - w ∑ d = 1 D Var ( U i , : ( d ) ) + α ∑ j = 1 k Var ( Y ( i , j ) ) ) subject to ∑ B = λ , ( 7 )', 'which is a standard constrained linear binary integer programming problem that can be easily solved by sorting the coefficients of B and taking the top λ entries .', 'The time complexity is O ( m × n × k ) for computing the construction error terms , O ( D × m × k ) for computing the variances and O ( m log m ) for sorting the coefficients .', 'The overall time complexity is O ( m × n × k ) assuming n × k > log m .', 'Thus , the total time complexity of each iteration in Algorithm 1 will be O ( ( m + λ ) × n × k ) , which is comparable to k-means when λ <<m .', 'There are four hyper-parameters to tune for the scVDMC algorithm , α and w: weights of the two variance terms , k: the number of clusters and λ: the number of marker genes .', 'Below we describe our strategies for tuning α , w and k assuming that λ can be approximately informed by prior knowledge of the cell types .', 'Tuning α: The role of α is to weight the cost term on the cross-domain variance of the cluster centers .', 'The larger the α the more similar the cluster centers are across the domains .', 'Ideally , α should be relatively small to allow smaller reconstruction error but yet meet the consistency requirement across the domains .', 'The strategy is to start with a small α and measure the total difference between the cluster centers of the corresponding cluster across the domains , and then increase α to reduce the difference until the total difference does not change much .', 'This selection can also be achieved by visualization of the cluster centers with Principle Component Analysis ( PCA ) or other dimension reduction methods .', 'After clustering , we can project the data in each domain into the first two PCs .', 'The distance between the cluster centers of the same cluster in each domain can be compared for choosing an appropriate α .', 'Several examples are shown later in the experiments .', 'Deriving the upper bound of w: Eq ( 3 ) is a sum of a few quadratic terms of variable U i , : ( d ) .', 'The global minimum of U i , : ( d ) can be solved in closed-form if the Hessian below is positive semi-definite ,', 'H = V ( d ) T V ( d ) - 2 w k Ψ + 2 α k Φ d , d d I k .', '( 8 ) In the following , we show that an upper bound on w will guarantee that H is positive semi-definite .', 'By Gershgorin circle theorem ( For any eigenvalue δ of matrix H , |δ − Hii| ≤ ∑j≠i |Hij| for ∀i ⇔ Hii − ∑j ≠ i |Hij| ≤ δ ≤ Hii + ∑j ≠ i |Hij| . ) , the sufficient condition of H ≽ 0 is Hii − ∑j ≠ i |Hij| ≥ 0 for ∀i .', 'This is equivalent to stating that H is diagonally dominant and only has non-negative diagonal entries .', 'H can be rewritten as follows ,', 'H i i = c i + 2 w ( 1 - k ) k 2 + 2 α k ( d - 1 ) d 2 , ∀ i = 1 , . . . , k H i j = 2 w k 2 , ∀ i ≠ j ,', 'where ci is the ith diagonal entry of matrix V ( d ) T V ( d ) , i . e . , the size of cluster i .', 'Then we have', 'c i + 2 w ( 1 - k ) k 2 + 2 α k ( d - 1 ) d 2 ≥ 2 w ( k - 1 ) k 2', 'and thus ,', 'w ≤ k 2 c m i n 4 ( k - 1 ) + α k 3 ( d - 1 ) 2 d 2 ( k - 1 )', 'where cmin is the minimum of ci , ∀i = 1 , … , k .', 'Since cmin ≥ 1 ( no empty cluster ) , we obtain a loose upper bound of w = k 2 4 ( k - 1 ) + α k 3 ( d - 1 ) 2 d 2 ( k - 1 ) .', 'In all the experiments , we set w to be smaller than the upper bound for feasible implementation .', 'Determining the number of clusters k: The number of clusters k is selected by an “elbow” plot of the within-clusters sum of squares Ts computed as follows:', 'T s = ∑ d = 1 D | | D B ( X i , : ( d ) - U i , : ( d ) V ( d ) T ) | | 2 2 .', '( 9 ) Ts represents the amount of variance to minimize for better clustering .', 'Larger k will lead to smaller Ts .', 'By plotting Ts under different options of k , we can select the best k at the so-called “elbow” of the curve .', 'S6 and S7 Figs show the “elbow” plot on two datasets in the experiments .', 'In addition , when an empty cluster is created , the calculation of cluster center variance will be invalid .', 'To address the possible issue , we use a simple splitting procedure to handle empty clusters .', 'Specifically , if there is an empty cluster in V ( d ) ( i . e . the whole column is 0 ) we randomly split the largest cluster into two clusters .', 'This procedure is repeated until there are exactly k clusters .', 'This strategy is similar to commonly used k-mean rerun when a cluster center is collapsed on a single data point or no data point .', 'To identify sub-populations producing homing signals that could attract bone marrow-derived cells to injured skin , we captured single dermal fibroblasts from six patients with severe generalized RDEB and their HLA-matched healthy siblings using the Fluidigm C1 system .', 'The demographics information of the patients and donors are shown in S1 Table .', 'Cell culture: Dermal fibroblasts from patients with severe generalized RDEB and their human leukocyte antigen ( HLA ) matched healthy siblings were obtained from skin biopsies and cultured in DMEM high glucose ( Thermo Fisher Scientific ) containing 10% fetal bovine serum ( MilliporeSigma ) , 1% Pen/Strep ( Thermo Fisher Scientific ) , 1% L-glutamine ( Thermo Fisher Scientific ) , and 1% MEM NEAA ( Thermo Fisher Scientific ) .', 'For sub-culture , the medium was removed and cells were washed with 1X PBS ( Thermo Fisher Scientific ) and detached using Trypsin/EDTA ( Thermo Fisher Scientific ) .', 'Experiments were performed with fibroblasts at passages 4-9 .', 'Single-cell capture and RNA-seq: Fibroblasts were collected by trypsinization and resuspended in 5 μL of fibroblast medium for loading into the capture chip .', 'The medium- ( 10-17 μm diameter ) and large-size ( 17-25 μm diameter ) chips were used to capture cells with the C1 system ( Fluidigm ) .', 'Cells were loaded at a concentration of 2 . 5 x 105 per μL and stained with the Live/Dead Viability/Cytotoxicity kit ( Thermo Fisher Scientific ) .', 'Cells were imaged with phase-contrast and fluorescence microscopy to assess cell number and viability at each capture point .', 'Capture sites with single , live cells were selected while capture sites with multiple , no , or an unclear number of cells were excluded from further analysis .', 'Images for each single-cell used in this study are available upon request .', 'In total , 295 patient cells and 248 sibling cells were selected .', 'On the device , cDNA was created from the selected cells using the SMARTer Ultra Low RNA kit designed for the C1 system ( Clontech ) .', 'mRNA libraries were constructed using the Nextera XT kit ( Illumina ) according to the manufacturer’s protocol .', 'The libraries were sequenced on an Illumina MiSeqv3 with 75bp paired-end reads to a depth of 19-22 million reads per lane .', 'Target sequencing depth for each library was 200K reads .', 'The RNA-seq data have been deposited in NCBI’s Gene Expression Omnibus and are accessible through GEO series accession number GSE108849 ( https://www . ncbi . nlm . nih . gov/geo/query/acc . cgi ? acc=GSE108849 ) .', 'Processing of RNA-seq data: Paired-end 75bp reads were mapped to the UCSC human transcriptome ( hg19 ) using Bowtie2 ( version 2 . 2 . 4 ) and Tophat ( version 2 . 0 . 9 ) .', 'Gene expression levels were calculated using Cuffquant ( Cufflinks version 2 . 2 . 1 with parameters -u -max-bundle-frags 10000000 ) and Cuffnorm ( Cufflinks version 2 . 2 . 1 ) .', 'FPKM values as estimated by Cufflinks were added a value of 1 ( to avoid zeros ) and log2 transformed .', 'We removed nine single-cell samples with low read counts ( < 50K ) and sub-sampled two single-cell samples sequenced as population controls with high read counts ( > 1 . 5M ) ( random sub-sampling , 10% of total reads ) .', '11 single-cell samples were excluded as outliers .', 'We excluded lowly expressed genes ( average log2 ( FPKM ) < 1 . 5 ) from further analysis .', 'The remaining 543 single-cell samples met the requirement of expressing at least 2 , 000 of these remaining 5 , 196 genes .', 'For each individual , the number of single-cells used in the analysis and the average number of reads for those single-cells is summarized in Table 1 .', 'The total number of the reads and the number of aligned reads in each cell are also shown in S3 Fig . Flow cytometry: Fibroblasts were collected by trypsinization ( as above ) and resuspended in fibroblast medium .', 'A BD Cytofix/Cytoperm™ kit ( BD Biosciences ) was used to prepare the cells for intracellular staining .', 'Cells were fixed for 15 min with 150 μl Fixation/Permeabilization solution before being resuspsended in 300 μl 1X BD Perm/Wash Buffer and incubated at 4°C for 20 min .', 'Primary antibodies ( S2 Table ) were diluted in 100 μl 1X BD Perm/Wash Buffer and cells were resuspended in this for 20-30 min at 4°C , followed by one wash with 500 μl 1X BD Perm/Wash Buffer .', 'Secondary antibodies ( S3 Table ) were diluted in 300 μl 1X BD Perm/Wash Buffer and cells were resuspended in this for 20-30 min at 4°C , followed by one wash with 500 μl 1X BD Perm/Wash Buffer , and resuspension in 300 μl 1X BD Perm/Wash Buffer .', 'Flow cytometry experiments were carried out on a BD LSRII system equipped with FACsDiva 8 . 0 software ( BD Biosciences ) and analyzed using FlowJo ( Tree Star Inc . ) .', 'Most existing methods focus only on sub-population clustering and differential gene expression detection among the learned cell clusters with one ( pooled ) cell population .', 'Some of these methods were directly adopted from traditional bulk RNA-seq analysis and/or classical dimension reduction algorithms such as Principal Component Analysis [6–8] , hierarchical clustering [9] , t-SNE [10–12] , Independent Component Analysis [13] and Multi-dimensional Scaling [14] .', 'Other methods focus on special properties of scRNA-seq data , such as high variance and uneven expressions .', 'For example , SNN-Cliq [15] uses a ranking measurement to get reliable results on high dimensional data; [16] proposed a special dimension reduction method to handle the large amount of zeros in scRNA-seq; [17] proposed a Latent Dirichlet Allocation model with latent gene groups to measure cell-to-cell distance; CellTree method [17] clusters single cells by a detected tree structure outlining the hierarchical relationship between single-cell samples to introduce biological prior knowledge; Seurat [18] was proposed to infer cellular localization by integrating single-cell RNA-seq data with in situ RNA patterns; and more recently a consensus clustering approach SC3 [19] was proposed to improve the robustness of clustering through combining multiple clustering solutions by consensus .', 'Mixed multiple batch strategies [9 , 20] have been proposed to reduce the technical variance , which does not directly improve clustering .', 'To the best of our knowledge , multitask clustering with an embedded feature selection has not been previously applied to scRNA-seq data analysis .', 'All patients gave consent for samples to be taken per the Declaration of Helsinki .', 'This research was approved by the University of Minnesota’s Institutional Review Board: IRB 1301M26601: MT2013-02R ( Establishment of a Cell and Tissue Repository for Human Cell Reprogramming and Derivation of iPS Cell Lines to Investigate Mechanisms and Treatment of Human Disease ) .'], ['scVDMC was compared with six baseline methods: ( 1 ) k-means clustering on each domain separately , ( 2 ) pooling all domains and applying k-means clustering , ( 3 ) SNN-Cliq [15] , ( 4 ) CellTree [17] , ( 5 ) Seurat [18] and ( 6 ) SC3 [19] .', 'Pooled k-means ( 2 ) was used to obtain the initialization for scVDMC .', 'To apply the SNN-Cliq method [15] , we used the provided MATLAB code to transform the data into the SNN graph , then used the Python code to produce the clustering result by ranking measurement .', 'There are three hyper-parameters: k ( size of the nearest neighbor list ) , r ( parameter for quasi-clique finding , range ( 0 , 1] ) , and m ( parameter for cluster merging range ( 0 , 1] ) .', 'We tested multiple combinations of the three hyper-parameters using k = 3 , 5 , 7 , r = 0 . 1 , 0 . 2 , … , 0 . 9 and m = 0 . 1 , 0 . 2 , … , 0 . 9 .', 'We also required the program to annotate all the data instead of leaving singletons unlabeled ( −n ) .', 'Since SNN-Cliq identifies the number of clusters automatically , we only reported the results with the correct number of clusters .', 'In all experiments with SNN-Cliq , we further removed genes with low expression and log-transformed the data , as recommended in [15] .', 'To apply the CellTree method [17] , we used the provided R package to first fit a Latent Dirichlet Allocation ( LDA ) model with the default method ( joint MAP estimation ) to choose the number of topics followed by learning a pair-wise distance for all cells .', 'Then we ran hierarchical clustering with four different methods for computing cluster distance ( ‘ward’ , ‘complete’ , ‘single’ , ‘average’ ) and selected the best clustering results .', 'To apply Seurat [18] , Seurat v2 . 0 R package was downloaded from SATIJA LAB .', 'The scRNA-seq data were converted into the required format ( gene index | cell index | gene expression ) as the input .', 'The parameter “Resolution” tunes the granularity of the downstream clustering , with increased values resulting a larger number of clusters .', 'We tested a range [0 . 5 , 1 . 5] to get the exact number of clusters for comparison with other methods .', 'The reported result of Seurat is computed with the resolution parameter that gives the exact number of clusters and the lowest error .', 'To apply the SC3 [19] we downloaded SC3 v1 . 7 . 2 R package from Bioconductor .', 'All parameters in SC3 are set to default .', 'In the experiments with more than 5000 instances for clustering , the SVM mode will be trigged to run a second stage supervised learning to improve the scalability .', 'To further test separated cluster , pooled clustering and SC3 combined with feature selection , we chose the genes with larger variance as the marker genes .', 'Since the other three baselines use a different strategy for clustering and do not provide marker-gene selection , we only focused on the clustering result for these three baselines .', 'The true cluster labels are obtained as the validated clusters with high confidence in the mESC data [21] and Lung data [22] , and the known PBMC populations from donor A sorted with FACS analysis [23] .', 'We downloaded the single-cell expression data for 250 mESCs [21] from the European Bioinformatics Institute’s ( EBI ) ESpresso database .', 'These 250 mESCs cultured in serum conditions were captured using the Fluidigm C1 on three different days from three different passages ( biological replicates , n = 81 , 90 , and 79 ) .', 'After removing genes expressed uniformly within a single replicate , 12 , 114 genes remained .', 'To tune α for scVDMC , we examined the positions of the cluster centers across the domains and show the visualization by PCA in S4 ( A ) & S4 ( B ) Fig . Based on the visualization , α = 0 and 1 are chosen since the relative positioning of cluster centers are similar in all the three domains .', 'For the SNN-Cliq method , we further removed genes with log-transformed average expression less than 20 .', 'Fig 3 ( A ) shows the clustering results .', 'Compared with the six baselines , scVDMC shows a consistently lower error with different choices of λs .', 'Within a reasonable range of λ , such as from 20 to 200 , scVDMC shows significant improvement compared with the baseline methods .', 'When λ is too small , such as 10 genes selected , there are not enough markers to capture the difference among the cell types such that the error is larger .', 'When λ is too big , scVDMC will consider almost all the genes and the variance selection will not play a role .', 'As such , scVDMC will eventually degrade into separated k-means and the error will also increase .', 'As shown in S1 ( A ) Fig , it is worth noting that the results are less sensitive to the choice of the parameter w , for which the upper bound for w is 9 8 in this case .', 'It is also interesting that the CellTree method performed better than both pooled and separated k-means , while SNN-Cliq and SC3 performed better than separated k-means but worse than pooled k-means .', 'Under various tuning of the parameters , Seurat still performed poorly on this dataset .', 'Both separated k-means and pooled k-means performed much worse with the feature selection by variance , indicating that simple feature selection strategies will not identify correct markers in this dataset .', 'Running scVDMC with α = 1 performed the best when 20 marker genes are selected but the overall performance is very similarly as running with α = 0 , indicating that the control of the cross-domain variance could play a role in improving the results .', 'However , since the cluster centers are already not very different when running with α = 0 , the improvement will only be marginal .', 'Fig 3 ( B ) shows the detailed clustering errors by scVDMC , pooled k-means and separated k-means .', 'Compared with the pooled k-means and separated k-means , scVDMC captures relatively high variance in the leading principle components and achieves improved clustering in every domain ( fewer mixed-color dots ) .', 'In S2 ( A ) Fig , we also show the convergence of scVDMC by the number of iterations .', 'Analysis of the mESC transcriptome data using scVDMC yielded comparable results on marker gene selection in the original paper [21] as well as pooled and separated k-means .', 'Both analyses were able to detect and highly rank the known markers of differentiation Krt8 , Krt18 , Anxa1 , Anxa3 , and Acta1 .', 'Further , scVDMC detected several additional genes that pooled k-means , separated k-means and the original paper did not .', 'These included Cav1 , which is required for normal lung development [24] and Dsp , variants of which are associated with idiopathic pulmonary fibrosis [25] .', 'We downloaded the single-cell expression data for 80 embryonic mouse lung epithelial cells [22] .', 'These 80 single-cell samples were taken from three different mice ( biological replicates , n = 20 , 34 , and 23 ) and contained five cell types: ciliated , Clara , AT1 , and AT2 cells , as well as a bi-potential progenitor ( BP ) .', 'Since only one replicate contained ciliated cells , we removed these from the analysis , leaving 77 single-cell samples .', 'After removing genes expressed uniformly within a single replicate , 7 , 357 genes remained .', 'To tune α for scVDMC , we examined the positions of the cluster centers across the domains and show the visualization by PCA in S4 ( C ) & S4 ( D ) Fig . α = 1 is chosen as the optimal parameter to achieve similar relative positioning of cluster centers in all the three domains .', 'For the SNN-Cliq method , we further removed genes with log-transformed average expression less than 2 .', 'With the limited number of single-cell samples in this dataset , scVDMC still much improved clustering over the baselines in the range of λ ∈ [30 , 80] shown in Fig 3 ( C ) .', 'In Fig 3 ( D ) , PCA plots of the top 50 genes show a trend similar to the ESC dataset , where scVDMC’s top genes capture more variance and show less clustering error .', 'Both SNN-Cliq and CellTree performed better than pooled k-means and separated k-means , with SNN-Cliq leading CellTree by a very small margin .', 'Similarly , Seurat also performed poorly while SC3 performed well on the dataset with only 5 mistakes .', 'It is also interesting to observe that running scVDMC with α = 1 performed significantly better than running with α = 0 , indicating that the control of the cross-domain variance played an important role in improving the results .', 'Since the cluster centers are very different when running with α = 0 , the improvement is significant .', 'Another interesting observation is that the clustering performance is more sensitive to the number of marker genes to select by scVDMC .', 'In particular , selection of 20-80 genes with scVDMC ( α = 1 ) will give the optimal clustering results while selection of more than 90 genes will give much higher error .', 'This is due to the small clusters in this dataset ( e . g . purple cluster in domain 2 and yellow cluster in domain 1 ) , which could be sensitive to the number of selected genes in low-read-coverage samples .', 'Thus , the error will be more sensitive to the gene selection in this small dataset .', 'On this dataset , both separated k-means and pooled k-means performed better with the feature selection by variance but never achieved zero clustering error as scVDMC does .', 'As shown in S1 ( B ) and S2 ( B ) Figs , scVDMC behaved similarly by the choices of the w parameters and the convergence .', 'Analysis of the mouse lung epithelial transcriptome data using scVDMC yielded comparable results in the original paper [22] as well as pooled and separated k-means .', 'Both analyses were able to detect and highly rank the known marker genes of the different cell types: Clara ( Scgb1a1 ) , AT1 ( Pdpn , Ager ) , and AT2 ( Sftpc , Sftpb ) .', 'Further , scVDMC detected several additional genes that pooled k-means , separated k-means and the original paper did not .', 'These included two components of the Notch signaling pathway ( Notch1 and Nrarp ) previously shown to be critical for the development of lung alveolar spaces , with AT2 cells being major sites of Notch activation [26] .', 'We downloaded the peripheral blood mononuclear cells ( PBMC ) data generated by [23] from the 10xGenomics website .', 'In the original data , there are 10 bead-enriched subpopulations of PBMC from a fresh donor ( Donor A ) with 93802 cells in total .', 'In addition , there are also PBMC from two other frozen donors ( Donor B and C ) with 7783 and 9519 cells , respectively .', 'A massive droplet-based method was applied to count the mRNAs in the tens of thousands of cells in parallel .', 'To better evaluate the multitask learning setting , we sampled from each of the 10 s
['Introduction', 'Materials and methods', 'Results', 'Discussion']
['Single-cell RNA sequencing ( scRNA-seq ) has been widely applied to discover new cell types by detecting sub-populations in a heterogeneous group of cells .', 'Since scRNA-seq experiments have lower read coverage/tag counts and introduce more technical biases compared to bulk RNA-seq experiments , the limited number of sampled cells combined with the experimental biases and other dataset specific variations presents a challenge to cross-dataset analysis and discovery of relevant biological variations across multiple cell populations .', 'In this paper , we introduce a method of variance-driven multitask clustering of single-cell RNA-seq data ( scVDMC ) that utilizes multiple single-cell populations from biological replicates or different samples .', 'scVDMC clusters single cells in multiple scRNA-seq experiments of similar cell types and markers but varying expression patterns such that the scRNA-seq data are better integrated than typical pooled analyses which only increase the sample size .', 'By controlling the variance among the cell clusters within each dataset and across all the datasets , scVDMC detects cell sub-populations in each individual experiment with shared cell-type markers but varying cluster centers among all the experiments .', 'Applied to two real scRNA-seq datasets with several replicates and one large-scale droplet-based dataset on three patient samples , scVDMC more accurately detected cell populations and known cell markers than pooled clustering and other recently proposed scRNA-seq clustering methods .', 'In the case study applied to in-house Recessive Dystrophic Epidermolysis Bullosa ( RDEB ) scRNA-seq data , scVDMC revealed several new cell types and unknown markers validated by flow cytometry .', 'MATLAB/Octave code available at https://github . com/kuanglab/scVDMC .']
['scRNA-seq enables detailed profiling of heterogeneous cell populations and can be used to reveal lineage relationships or discover new cell types .', 'In the literature , there has been little effort directed towards developing computational methods for cross-population transcriptome analysis of multiple single-cell populations .', 'The cross-cell-population clustering problem is different from the traditional clustering problem because single-cell populations can be collected from different patients , different samples of a tissue , or different experimental replicates .', 'The accompanying biological and technical variation tends to dominate the signals for clustering the pooled single cells from the multiple populations .', 'In this work , we have developed a multitask clustering method to address the cross-population clustering problem .', 'The method simultaneously clusters each individual cell population and controls variance among the cell-type cluster centers within each cell population and across the cell populations .', 'We demonstrate that our multitask clustering method significantly improves clustering accuracy and marker discovery in three public scRNA-seq datasets and also apply the method to an in-house Recessive Dystrophic Epidermolysis Bullosa ( RDEB ) dataset .', 'Our results make it evident that multitask clustering is a promising new approach for cross-population analysis of scRNA-seq data .']
['flow cytometry', 'sequencing techniques', 'medicine and health sciences', 'cell cycle and cell division', 'applied mathematics', 'cell processes', 'marker genes', 'fibroblasts', 'simulation and modeling', 'algorithms', 'developmental biology', 'organism development', 'mathematics', 'connective tissue cells', 'molecular biology techniques', 'lung development', 'rna sequencing', 'research and analysis methods', 'animal cells', 'gene expression', 'connective tissue', 'biological tissue', 'organogenesis', 'molecular biology', 'spectrophotometry', 'cytophotometry', 'cell biology', 'anatomy', 'genetics', 'biology and life sciences', 'cellular types', 'physical sciences', 'spectrum analysis techniques']