BayCount: A Bayesian Decomposition Method for Inferring Tumor Heterogeneity using RNA-Seq Counts
Fangzheng Xie Mingyuan Zhou Yanxun Xu



[
    [
    December 30, 2023
=====================





Tumor is heterogeneous – a tumor sample usually consists of a set of subclones with distinct transcriptional profiles and potentially different degrees of aggressiveness and responses to drugs. Understanding tumor heterogeneity is therefore critical to precise cancer prognosis and treatment. 
In this paper, we introduce BayCount, a Bayesian decomposition method to infer tumor heterogeneity with highly over-dispersed  RNA sequencing count data. 




Using negative binomial factor analysis, 
BayCount takes into account both the between-sample  and gene-specific random effects 
on raw counts of sequencing reads mapped to each gene.



For posterior inference, we develop an efficient compound Poisson based 
blocked Gibbs sampler.  
Through extensive simulation studies and analysis of The Cancer Genome Atlas lung cancer and kidney cancer RNA sequencing count data, we  show that 

BayCount is able to accurately estimate the number of subclones,  the proportions of these subclones in each tumor sample, and the gene expression profiles in each subclone. Our method represents the first effort in characterizing tumor heterogeneity using RNA sequencing 
 count data that simultaneously removes the need of normalizing the counts, 
achieves statistical robustness, and obtains biologically/clinically meaningful insights. 

KEY WORDS:  Cancer genomics,  compound Poisson,  Markov chain Monte Carlo, negative binomial, over-dispersion 




§ INTRODUCTION



Tumor heterogeneity (TH) is a phenomenon that describes distinct molecular profiles of different cells in one or more tumor samples. 
TH arises during the formation of a tumor as a fraction of cells acquire and accumulate different somatic events (e.g., mutations in different cancer genes), resulting in heterogeneity within the same biological tissue sample and between different ones,

spatially and temporally <cit.>. As a result, tumor cell populations are composed of different subclones (subpopulations) of cells, characterized by distinct genomes, transcriptional profiles <cit.>, as well as other molecular profiles, such as copy number alterations.  
Understanding TH is critical to precise cancer prognosis and treatment. Heterogenetic tumors may exhibit different degrees of aggressiveness and responses  to drugs among different samples due to genetic or gene expression differences. The level of heterogeneity itself can be used as a biomarker to predict treatment response or prognosis since more heterogeneous tumors are more likely to contain treatment-resistant subclones <cit.>. 
This will ultimately facilitate the rational design of combination treatments, with each distinct compound targeting a specific tumor subclone based on its transcriptional profile.

Large-scale sequencing techniques provide  valuable information for understanding tumor complexity and open a door for the desired statistical inference on TH. Previous studies have focused on reconstructing the subclonal composition by quantifying  the  structural subclonal copy number variations  <cit.>,  somatic mutations   
<cit.>, or both <cit.>. In this paper, we aim to learn tumor transcriptional heterogeneity using RNA sequencing (RNA-Seq) data. 

In the analysis of  gene expression data, matrix decomposition models have been extensively studied in the context of microarray and normalized RNA-Seq data <cit.>. Generally, given gene expression data matrix X=(x_ij)_G× S, where the (i,j)th element records the expression value of the ith gene in the jth sample, they decompose X by modeling  x_ij with ∑_k=1^Kϕ_ikθ_kj,
where ϕ_ik encodes the expression level of the ith gene in the kth subclone, θ_kj represents the mixing 
weight of the kth subclone in the jth sample, and K is the number of subclones. The decomposition can be solved by either optimization algorithms <cit.> or statistical inference by assuming a normal distribution on x_ij. While it is reasonable to assume normality for microarray gene expression data, it is often inappropriate to adopt such an assumption for directly modeling RNA-Seq data, which involve nonnegative integer observations. 
If 
a model based on normal distribution is used,  one often needs to first normalize RNA-Seq data 
 before performing any downstream analysis.  See <cit.> for a review on normalization methods. 
 
 Although normalization often destroys the nonnegative and discrete nature of the RNA-Seq data, it remains the predominant way for data preprocessing 
 
 
 due to not only the computational convenience in modeling normalized data, but also the lack of appropriate count data models. Distinct from previously proposed methods, in this paper, we  propose an attractive class of count data models in 
 decomposing RNA-Seq count matrices. 











There are, nevertheless, statistical challenges with RNA-Seq count data. First, the distributions of the RNA-Seq count data 

are typically over-dispersed and  sparse. Second, the scales of the read counts in sequencing data across samples can be enormously different  due to the mechanism of the sequencing experiment such as the variations in technical lane capacities. The larger the library sizes (i.e., sequencing depth) are, the larger the read counts tend to be. In addition, the differences in gene lengths or GC-content <cit.> can bias 
gene differential expression analysis, particularly for lowly expressed genes <cit.>. A number of count data models have been developed for RNA-Seq data <cit.>. For example, 
<cit.>  proposed a Poisson factor model on microRNA to reduce the dimension of count data and identify low-dimensional features, followed by a clustering procedure over tumor samples. <cit.> developed a method using a mixture of negative binomial and Poisson distributions to model single cell RNA-Seq data for gene differential expression  analysis. None of these methods, however, address the problem of TH. 




To this end, we propose BayCount, a Bayesian 

matrix decomposition model 
built upon the negative binomial model  <cit.>, to infer tumor transcriptional heterogeneity using RNA-Seq count data. BayCount accounts for 
both the  between-sample and gene-specific random effects and infers the number of latent subclones, the proportions of these subclones in each sample, and subclonal expression simultaneously.  





The remainder of the paper is organized as follows. 
 In Section <ref>, we  introduce BayCount, a hierarchical Bayesian model for RNA-Seq  count data, and develop an  efficient compound Poisson based blocked Gibbs sampler. We investigate the performance of posterior inference and robustness of the BayCount model through extensive simulation studies in Section 3, and apply our proposed BayCount model to analyze two 
 
 real-world RNA-Seq datasets from The Cancer Genome Atlas (TCGA) <cit.>
 in Section 4. We conclude the paper 
 in Section <ref>.






§ HIERARCHICAL BAYESIAN MODEL AND INFERENCE

In this section we  present the proposed hierarchical model for RNA-Seq count data, develop the corresponding posterior inference, and discuss how to 
determine the number of subclones. 



 §.§ BayCount Model

We assume that S 
 tumor samples are available from the same or different patients. 
Consider a G× S count matrix Y=(y_ij)_G× S, where each row represents a gene, each column represents a tumor sample, and the element y_ij records the read count of the ith gene from the jth tumor sample. The Poisson distribution Pois(λ) with mean λ>0 is commonly used for modeling count data.

Poisson factor analysis (PFA) <cit.>  factorizes the count matrix Y as 

y_ij∼(∑_k=1^Kϕ_ikθ_kj), 

where Φ=(ϕ_ik)_G× K∈ℝ_+^G× K is 
the factor loading matrix 
and Θ=(θ_kj)_K× S∈ℝ_+^K× S is the factor score matrix. 
Here K is an integer indicating the number of latent factors, and each column of Φ is subject to the constraint that ∑_i=1^Gϕ_ik=1 and ϕ_ij≥ 0. However, the restrictive equidispersion property of the Poisson distribution that the variance and mean are the same limits the application of PFA in modeling sequencing data, which are often highly over-dispersed.
For this reason, one may consider 

negative binomial factor analysis (NBFA)  of <cit.> that factorizes Y as

y_ij∼NB(∑_k=1^Kϕ_ikθ_kj,p_j), where p_j∈(0, 1). 

We denote y∼NB(r, p) as a negative binomial distribution with shape parameter r>0 and success probability p∈(0,1), whose 

mean and variance 

are rp/(1-p) and rp/(1-p)^2, respectively, with the variance-to-mean ratio as 1/(1-p).













Denote the jth column of  as _j=(y_1j, y_2j, …, y_Gj)^T,  the count profile of the jth tumor sample. 
To account for both the between-sample and gene-specific random effects when modeling RNA-Seq count data, we  propose 

    y_ij|λ,α_i,ζ_j, p_j,Φ,Θ   ∼   NB(λα_i+∑_k=1^Kϕ_ikθ_kjζ_j, p_j),

where α_i accounts for the gene-specific random effect of the ith gene, λ and p_j control the scales of the gene-specific effect and between-sample effect of the jth sample, respectively,  and ∑_k=1^Kϕ_ikθ_kjζ_j represents the average effect of the K subclones on the expression of the ith gene in the jth sample.

To see this, recall that the mean of y_ij
based on (<ref>) is

    𝔼[y_ij]=(λα_i+∑_k=1^Kϕ_ikθ_kjζ_j) p_j/1-p_j.

Since p_j is sample-specific, the term 
p_j/(1-p_j)
describes the effect of sample j on read counts due to technical or biological reasons 

(e.g., different library sizes, biopsy sites, etc). 



We assume the relative expression of the ith gene in the kth subclone is described by ϕ_ik, where ϕ_ik≥ 0.  
Since the sample-specific effect has already been captured by p_j, for modeling convenience, we normalize the gene expression so that the expression levels sum to one for each subclone. Namely, ∑_i=1^Gϕ_ik=1 for all k=1,⋯,K. Furthermore, we assume that θ_kj represents the proportion of the kth subclone in the jth sample, where θ_kj≥ 0 and ∑_k=1^Kθ_kj=1. We can interpret θ_kjζ_j as the population frequencies of the kth subclone in the jth sample, where parameter ζ_j controls the scale. 


Together, the summation ∑_k=1^Kϕ_ikθ_kjζ_j represents the aggregated expression level of the ith gene across all K subclones for the jth sample. 

To further account for 

the gene-specific random effects that are 


 independent of the samples and subclones, we introduce an additional term λα_i to describe the random effect of the ith gene on the read counts such as GC-content and gene length. 
 We assume ∑_i=1^Gα_i=1 so that α_i represents the relative gene-specific random effect of the ith gene 
with respect to all the genes and λ controls the overall scale of the gene-specific random effects. 


Following <cit.>, the model in (<ref>) has an augmented representation as

    y_ij   =   x_ij+z_ij, 
    
    x_ij   =   ∑_k=1^Kx_ijk, 
    
     z_ij|λ,α_i,p_j   ∼   NB(λα_i,p_j), 
    
    x_ijk|_k,_j,ζ_j, p_j   ∼   NB(ϕ_ikθ_kjζ_j,p_j).

From (<ref>), 


the raw count y_ij of the ith gene in the jth sample can be interpreted as coming from multiple sources: x_ijk represents the count of the ith gene contributed by the kth subclone in the jth sample, where k=1, …, K, while z_ij is the count contributed by the gene-specific random effect of the ith gene in the jth sample.




Denote y_ j = ∑_i=1^Gy_ij. Since ∑_i=1^G ϕ_ik=1 and ∑_k=1^Kθ_kj = 1 by construction,  under (<ref>), by the additive property of independent negative binomial random variables with the same 
success probability, we have

    y_ j|λ,α_i,ζ_j, p_j,Φ,Θ∼NB(λ+ζ_j, p_j),





and, in particular, the mean as 𝔼[y_ j] = (λ + ζ_j)p_j/(1-p_j) and the variance as Var(y_ j)= 𝔼[y_ j] + 𝔼^2[y_ j]/(λ+ζ_j). 






It is clear that p_j, the between-sample random effect of the jth sample,   governs the variance-to-mean ratio of y_ j, whereas λ+ζ_j, the sum of the scale λ of 
the gene-specific random effects and the scale ζ_j for the jth sample,  controls the quadratic relationship between Var(y_ j) and  𝔼[y_ j].










We complete the model by setting the following priors that will be shown to be amenable to posterior inference: 


    _k   ∼Dirichlet(η,⋯,η),        ∼Dirichlet(δ,⋯,δ),
    _j| r_1,⋯,r_K   ∼Dirichlet(r_1,⋯,r_K),         p_j∼Beta(a_0,b_0),
    ζ_j| r_1,⋯,r_K,c_j   ∼Gamma(∑_k=1^Kr_k,c_j^-1),        λ∼Gamma(u_0,v_0^-1),

where _k=(ϕ_1k,⋯,ϕ_Gk)^T, _j=(θ_1j,⋯,θ_Kj)^T, =(α_1,⋯,α_G)^T, Gamma(a, b) denotes a gamma distribution with mean ab and variance ab^2, and Dirichlet(η_1,⋯,η_d) denotes a d-dimensional Dirichlet distribution with parameter vector (η_1,⋯,η_d). We further impose the hyperpriors, expressed as r_k|γ_0, c_0∼Gamma(γ_0/K,c_0^-1), c_j∼Gamma(e_0,f_0^-1), γ_0∼Gamma(g_0,h_0^-1), and c_0∼Gamma(e_0,f_0^-1)
to construct a more flexible model.














Shown in Figure <ref> 
is  the 
graphical representation of our BayCount model. 





 §.§ Gibbs Sampling via Data Augmentation

 

For the proposed BayCount model, while the conditional posteriors of p_j, c_j and c_0 
are straightforward to derive due to conjugacy, 




a variety of data augmentation techniques are used to derive the closed-form Gibbs sampling update equations
  for all the other model parameters. 


 Rather than going into the details here, let us first assume that we have already sampled the latent counts x_ijk given the observations y_ij and model parameters, which, according to Theorem 1 of  <cit.>, can be realized by sampling from the Dirichlet-multinomial distribution; given x_ijk, we show how to derive the Gibbs sampling update equations for Φ and Θ via data augmentation; and we will describe in the Supplementary Material a compound Poisson based blocked Gibbs sampler that completely removes the need of sampling x_ijk.






  §.§.§ Sampling Φ and Θ

We introduce an auxiliary variable ℓ_ijk that follows a Chinese restaurant table (CRT) distribution, denoted by ℓ_ijk| x_ijk,ϕ_ikθ_kjζ_j∼CRT(x_ijk, ϕ_ikθ_kjζ_j), with probability mass function

    p(ℓ_ijk| x_ijk,ϕ_ikθ_kjζ_j) = Γ(ϕ_ikθ_kjζ_j)/Γ(x_ijk+ϕ_ikθ_kjζ_j)|s(x_ijk, ℓ_ijk)| (ϕ_ikθ_kjζ_j)^ℓ_ijk,

supported on {0,1,2,⋯,x_ijk},
where 
s(x_ijk, ℓ_ijk) are Stirling numbers of the first kind <cit.>. Sampling ℓ∼CRT(x,r) 
can be realized by 
taking the summation of m independent Bernoulli random variables: ℓ=∑_t=1^x b_t, where b_t∼Bernoulli(r/(r+t-1)) independently. Following <cit.>, the joint distribution of ℓ_ij and x_ij described by 

    ℓ_ijk| x_ijk ,ϕ_ik,θ_kj,ζ_j
       ∼   CRT(x_ijk,ϕ_ikθ_kjζ_j),
    
    x_ijk|ϕ_ik,θ_kj,ζ_j,p_j
       ∼   NB(ϕ_ikθ_kjζ_j,p_j),

can be equivalently characterized under the compound Poisson representation

    x_ijk|ℓ_ijk ,p_j
       ∼   SumLog(ℓ_ijk,p_j),
    ℓ_ijk|ϕ_ik,θ_kj,ζ_j, p_j
       ∼   Pois(-ϕ_ikθ_kjζ_jlog(1-p_j)),

where x∼SumLog(ℓ,p) denotes the sum-logarithmic distribution generated as 
x = ∑_t=1^ℓu_t, where (u_t)_t=1^ℓ are independent, and identically distributed (i.i.d.) according to the logarithmic distribution <cit.> with probability mass function p(u)=-p^u/[ulog(1-p)], supported on {1,2,⋯}. 

Under this augmentation, the likelihood of ϕ_ik, θ_kj and ζ_j becomes 

    ℒ(ϕ_ik,θ_kj,ζ_j)∝(ℓ_ijk|-ϕ_ikθ_kjζ_jlog(1-p_j)),

where (·|λ) denotes the probability mass function of the Poisson distribution with mean λ .  It follows immediately 
 that the full conditional posterior distributions for _k and _j are

    (_k|-)   ∼   Dirichlet(η+
    ∑_j=1^Sℓ_1jk,⋯,η+∑_j=1^Sℓ_Gjk),
    
    (_j|-)   ∼   Dirichlet(
    r_1+∑_i=1^Gℓ_ij1,⋯,r_K+∑_i=1^Gℓ_ijK).




Using data augmentation, we can similarly  derive the full conditional posterior distributions for ζ_j, , r_k and γ_0, as described in detail in the Supplementary Material.

 



 §.§ Determining the Number of Subclones K

We have so far assumed a priori that K is fixed. Determining the number of factors in factor analysis 
is, in general, challenging. <cit.> suggested adaptively truncating  K during Gibbs sampling iterations. 
This adaptive truncation procedure, which is designed to fit the data well, may tend to choose a large number of factors, some of which may be highly correlated to each other and hence 
appear to be redundant. 
To facilitate the interpretation of the model output,
we seek a model selection procedure that estimates K in a more conservative manner. 
To select a moderate K that is large enough to fit the data reasonably well, but at the same time is small enough for the sake of interpretation, 


we generalize the deviance reduction-based approach in <cit.> and calculate the estimated log-likelihood of the model under different numbers of subclones using post-burn-in MCMC samples. These samples are obtained by running the compound Poisson based blocked Gibbs sampler for different K's. 
The estimate of K can be  identified by an apparent decrease in the slopes of segments that connect the log-likelihood values of two consecutive K values. Formally, we denote the log-likelihood function logℒ(K) as a function of K, and define the second-order finite difference Δ^2logℒ(K) of the log-likelihood function by

Δ^2logℒ(K):=2logℒ(K) -logℒ(K-1) - logℒ(K+1),

 for K=K_min+1,⋯,K_max-1, where K_min and K_max are the lower and upper limits of K, respectively.  
 Then an estimate of K is given by 

    K̂ = argmax_KΔ^2logℒ(K).
 








§ SIMULATION STUDY


In this section, we evaluate the proposed BayCount model through simulation studies.  Two different scenarios are considered. 


	
  * Scenario I:  We simulate the data according to the BayCount model itself in (<ref>).  In particular, we generate the subclone-specific 
	gene expression data matrix Φ = (ϕ_ik)_G× K^o∈ℝ_+^G× K^o by i.i.d. draws of _k∼Dirichlet(0.05,⋯,0.05), 
	
	the proportion matrix Θ = (θ_kj)_K^o× S by i.i.d. draws of _j∼Dirichlet(0.5,⋯,0.5), and ζ_j by i.i.d. draws of ζ_j∼Gamma(0.5K^o,1), where i=1, ⋯, G, j=1, ⋯, S, and k=1, ⋯, K^o. 
Here  G is the number of genes, S is the number of samples, and K^o is the simulated number of subclones.  
	We set λ = 1, draw  from Dirichlet(0.5,⋯,0.5), and generate p_j from a uniform distribution  such that the variance-to-mean ratio p_j/(1-p_j) of y_ j ranges from 100 to 10^6, encouraging the simulated data to be over-dispersed.  
	
  * Scenario II: To evaluate the robustness of the BayCount model, under scenario II we consider simulating the data from a model that is different from the BayCount. 
	We generate the subclone-specific 
	
	gene expression data matrix W = (w_ik)_G× K^o∈ℝ_+^G× K^o by i.i.d. draws of w_ik∼Gamma(0.05,10), 
	and the proportion matrix Θ = (θ_kj)_K^o× S by i.i.d. draws of _j∼Dirichlet(0.5,⋯,0.5). 
	We set λ=1,  draw  from Dirichlet(0.5,⋯,0.5), and generate p_j from a uniform distribution such that the variance-to-mean ratio p_j/(1-p_j) of y_ j ranges from 100 to 10^6. 
	
	The count matrix Y=(y_ij)_G× S is generated from 
	y_ij∼NB(λα_i+∑_k=1^K^ow_ikθ_kj,p_j). Note 
	that in scenario II the scales of W=(w_ik)_G× K^o are not subject to the constraint ∑_i=1^Gw_ik=1. 





We will show that  BayCount  can accurately recover both the subclone-specific gene expression patterns and subclonal proportions. The hyperparameters are set to be η = 0.1, a_0=b_0=0.01, e_0=f_0=1, g_0=h_0=1, and u_0=v_0=100. We consider K∈{2,3,⋯,10}. The compound Poisson based blocked Gibbs sampler is implemented with an initial burn-in of B=1000 iterations and a total of n=2000 iterations. The posterior means and  95% credible intervals for all parameters are computed using the 1000 post-burn-in MCMC samples.




 §.§ Synthetic data with K^o=3



We first simulate two 
datasets with G=100, S=20, and K^o=3 under both scenario I and scenario II. Under scenario I, the data generation scheme is the same as the BayCount model. Figure S1 in the Supplementary Material plots  Δ^2logℒ(K) versus K, indicating K̂=3, which is the same as the simulation truth. The estimated subclone-specific gene expression matrix Φ̂ and subclonal proportions Θ̂ are computed as the posterior means of the post-burn-in MCMC samples. 
Figure S2 and S3 compare the simulated true Φ and Θ with the estimated Φ̂ and Θ̂, respectively. We can see that both the subclone-specific gene expression patterns and the subclonal proportions are successfully recovered. 
 


The analysis under scenario II is of greater interest, since the focus is to evaluate the robustness of BayCount. BayCount yields an estimate of K̂=3, as shown in Figure S4. 









We then focus on the posterior inference based on K̂=3. 

Figure <ref> compares the estimated subclonal proportions Θ̂ with the simulated true subclonal proportions across samples,

 along with the posterior 95% credible intervals. The results show that the estimated Θ̂ approximates 
  the simulated true Θ well.

  





We then report the posterior inference on the subclone-specific gene expression Φ. Under the BayCount model, ∑_i=1^Gϕ_ik=1,  
hence the estimated Φ̂ by BayCount and the unnormalized gene expression profile matrix W used in generating the simulated data  are not directly comparable. 
To see whether the gene expression pattern is recovered, we first normalize W by its column sums as Ŵ = WΛ^-1, where 
Λ=diag(∑_i=1^Gw_i1,⋯,∑_i=1^Gw_iK), so that ŵ_ik represents the relative expression level of the ith gene in the kth subclone, and then compare Φ̂ with Ŵ. 

For visualization, the genes with small standard deviations (less than 0.01)

are filtered out due to their indistinguishable expressions across different subclones.  
Figure <ref> compares the heatmap of Φ̂, with the heatmap of the simulated true (normalized) subclone-specific gene expression Ŵ on selected differentially expressed genes. 
It is clear 
 that the pattern of subclone-specific gene expression estimated by BayCount 
 closely matches the simulation truth.  


















 §.§ Synthetic data with K^o=5

Similarly as in Section <ref>, we simulate  two datasets with G=1000, S=40, and K^o=5 under scenarios I and  II, respectively. Under scenario I, BayCount yields an estimate of K̂=5 (Figure S5), and from 

Figures S6 and S7, both the subclone-specific gene expression pattern and the subclonal proportions are successfully captured. 


Under scenario II, BayCount yields an estimate of K̂=5 (Figure S8). 


For the subclonal proportions Θ=(θ_kj)_K× S, Figure  <ref> shows that the estimated Θ̂ successfully recovers the simulated true proportions. Notice that the credible bands are narrower than those in Figure <ref>, implying relatively smaller variability in estimating subclonal proportions for larger dataset. Figure S9

presents the autocorrelation plots of the posterior samples of some randomly selected proportions  by the compound Poisson based blocked Gibbs sampler, indicating  that the Markov chains mix well. 



 
Figure <ref> compares the simulated true (normalized) subclone-specific gene expression Ŵ
with the estimated Θ̂ under the inference of BayCount. 
For this dataset we pre-screen Ŵ with a threshold 0.008 on the across-subclone standard deviation for all genes for visualization. 
The high concordance between the heatmaps of the estimated and true expression patterns of the  differentially expressed genes indicates that the subclone-specific gene expression patterns have been successfully recovered as well.  





















In summary, the BayCount model can accurately identify the number of subclones, estimate the subclonal proportions in each sample, and recover the subclone-specific gene expression pattern of the differentially expressed genes. 






§ REAL-WORLD DATA ANALYSIS



We implement and evaluate the proposed BayCount model on the RNA-Seq data from The Cancer Genome Atlas (TCGA) <cit.> to study tumor heterogeneity (TH) in both lung squamous cell carcinoma (LUSC) and kidney renal clear cell carcinoma (KIRC). We first run the proposed Gibbs sampler for each fixed K∈{2,3,⋯,10}, compute both the posterior mean and  95% credible interval of the log-likelihood for each fixed K, and estimate K by maximizing Δ^2logℒ(K) over K. Next, based on the estimated K̂ and the posterior samples generated by the proposed Gibbs sampler, we estimate the proportions of the identified subclones in each tumor sample and 

the subclone-specific gene expression,  which in turn can be used for a variety of downstream analyses.


































 §.§ TCGA LUSC Data Analysis

We apply the proposed BayCount model to the TCGA RNA-Seq data in lung squamous cell carcinoma (LUSC), which is a common type of lung cancer that  causes nearly one million

deaths worldwide every year. The raw RNA-Seq data for 200 LUSC tumor samples were processed using the Subread algorithm in the Rsubread package <cit.> to obtain gene level read counts <cit.>. We select 382 previously reported important lung cancer genes <cit.> for analysis, such as KRAS, STK11, BRAF, and RIT1. 


BayCount yields an estimate of five subclones (Figure S10) and their proportions in each tumor sample are shown in Figure <ref>. To identify the dominant subclone for each sample, we compare the estimated Θ̂ of the five subclones in each tumor sample, and use them to cluster the patients. Formally, for each patient j=1,⋯,S, we compute the dominant subclone k_j=_k=1,⋯,Kθ̂_kj, and then cluster patients according to {j: k_j=k}, k=1, …, K̂. 
 


That is to say, the patients with the same dominant subclone belong to the same cluster. We next check if the identified subclones have any clinical utility, e.g., stratification of patients in terms of overall survival.  Figure <ref>a shows the Kaplan-Meier plots of the overall survival of the patients in the five clusters identified by their dominant subclones. Indeed, patients stratified by these five BayCount-identified groups exhibit very distinct survival patterns (log-rank test p = 0.0194).
 
 




Figure <ref>b shows the expression levels of the top 30 differentially expressed genes  (ranked by the standard deviations of the  subclone-specific gene expression levels ϕ_ik's in an increasing order) in these five subclones. Distinct expression patterns are observed among  different subclones. For example, the FTL level is elevated in subclone 1; the expression levels of several genes encoding keratins (KRT5, KRT6A, etc.) are elevated in subclone 3; and the COL1A1 and COL1A2 expression levels are elevated in subclone 4. 
Interestingly, the patients with these dominant subclones also show the expected survival patterns. The subclone-1 dominated patients have better overall survival. Previous studies show 
 that the expression of FTL is decreased in lung tumors compared to normal tissues <cit.>, and one plausible explanation is that subclone 1 may descend from less malignant cells and therefore resemble (or 
  consist  
  of) normal cells. Keratins and collagen I (encoded by COL1A1 and COL1A2) are known to play key roles in epithelial-to-mesenchymal transition (EMT), which subsequently initiates metastasis and promotes tumor progression <cit.>. 
  
 This agrees with our observation of worse prognosis in patients who have either subclone 3 (with elevated Keratin-coding genes) or subclone 4 (with elevated collagen I coding genes) as their dominant subclone.
 




 §.§ Kidney Cancer (KIRC) Data Analysis


Similarly, we obtain gene level read counts <cit.> for 200 TCGA kidney renal clear cell carcinoma (KIRC)  tumor RNA-seq samples and analyze them with BayCount. 
Among a total of 23,368 genes, 966 significantly mutated genes <cit.> in KIRC patients are selected, including VHL, PTEN, MTOR, etc. 

BayCount yields an estimate of five subclones in KIRC (Figure S11). Figure <ref> shows the Kaplan-Meier curves of the overall survival of the patients grouped by their dominant subclones (panel a) and the heatmap of the top 30 differentially expressed genes (panel b). Since we have a large number of genes to begin with, whereas ∑_i=1^Gϕ_ik=1 for all k=1,⋯,K, the subclone-specific gene expression estimates Φ̂ will be small. For better visualization, we plot Φ̂ in the logarithmic scale. 



The subclonal proportions across 200 KIRC tumor samples are shown in  Figure S12. 
As shown in Figure <ref>, the patients with these dominant subclones again show distinct survival patterns. One of the poor survival groups (dominated by subclone 5) is characterized by elevated expression of TGFBI, which is known to be associated with poor prognosis <cit.> and 
matches our observation here. 



One distinction of our method from conventional subgroup analysis methods is that we focus on characterizing the underlying subclones (i.e., biologically meaningful subpopulations), by not only their individual molecular profiles but also their proportions. Instead of grouping the patients by their dominant subclones, we examine the proportion itself in terms of clinical utility. Interestingly, as shown in Figure <ref>a, the proportion of subclone 2 increases with tumor stage: i.e., as subclone 2 expands and eventually outgrows other subclones, the tumor becomes more aggressive. In contrast, the proportion of subclone 3 decreases with tumor stage (Figure <ref>b). Subclone 3 might be characterized by the less malignant (or normal-like) cells and takes more proportion in the beginning of the tumor life cycle. 
As tumor progresses to more advanced stages, subclone 3 could be 
suppressed by more aggressive subclones (e.g., subclone 2) and takes a decreasing proportion. Unsurprisingly,

the survival patterns 
 agree with our speculations about subclones 2 and 3, with 
the patients dominated by subclone 2 (the more aggressive subclone) and subcolone 3 (the less aggressive subclone) showing the worst and best survivals, respectively. 


More excitingly, we find that the proportions of these two subclones can complement clinical variables in further stratifying patients. For patients at early stage where the event rate is low and clinical information is relatively limited, the proportions of subclones 2 and 3 serve as a potent factor in further stratifying  patients (Figure S13) when dichotomizing at a natural cutoff. Combining our observations above, subclone proportions may provide additional insights into the progression course of tumors, assistance in biological interpretation, and potentially more accurate clinical prognosis. 




















































§ CONCLUSION
 

The emerging high-throughput sequencing technology provides us with massive information for understanding tumors' complex microenvironment and allows us to develop novel statistical models for inferring tumor heterogeneity. Instead of 
normalizing RNA-Seq data that may bias 
downstream analysis, we propose BayCount 
 to directly analyze the raw RNA-Seq count data. Overcoming the natural challenges of analyzing raw RNA-seq count data, BayCount is able to factorize 

them while adjusting for both the between-sample and gene-specific 
 random effects.
Simulation studies show that BayCount can accurately recover the subclonal inference used to generate the simulated data. 

We apply BayCount to the TCGA LUSC and KIRC datasets, followed by correlating the subclonal inferences with their clinical utilities for comparison. In particular, by grouping patients according to their 
dominant subclones, we observe distinct and biologically sensible overall survival patterns for both LUSC and KIRC patients. Moreover, the proportions of the subclones may complement clinical variables in further stratifying patients.
In addition to prognosis value, tumor heterogeneity may be used as a biomarker to predict treatment response. 
For example, tumor samples with large proportions of cells bearing higher expressions on clinically actionable genes should be treated differently from those that have no or a small proportion of such cells. In addition, metastatic or recurrent tumors may possess very different compositions of subclones and should be treated differently. 




BayCount provides a general framework for inference on latent structures arising naturally in many other biomedical applications involving count data. For example, analyzing single-cell data is a potential further application of BayCount due to their sparsity and over-dispersion nature. 
<cit.> describe Drop-Seq, a technology for profiling more than 40,000 single cells at one time. The unique characteristic of dropped-out events <cit.> in single cell sequencing limits the applicability of normalization methods in bulk RNA-Seq data. Also, 
such huge amount number of single-cells and high levels of sparsity pose difficulties for dimensionality reduction methods such as principal component analysis. 

Inferring distinct cell populations in single-cell RNA count data will be an interesting extension of BayCount. 








§ ACKNOWLEDGEMENT

Yanxun Xu's research is partly supported by Johns Hopkins inHealth and Booz Allen Hamilton. 


apalike




