iASeq: integrative analysis of allele-specificity of protein-DNA interactions in multiple ChIP-seq datasets

Data structure

Suppose there are D ChIP-seq datasets generated using cells from the same individual and the same cell type. Each dataset d corresponds to one TF or HM, and has J _d replicate samples (Figure 1a). Different datasets represent different TFs or HMs, or data generated by different labs. For the individual in question, assume one is interested in analyzing I heterozygote SNPs with known genotypes. We want to know whether the two alleles of each SNP behave differently in each dataset, and if possible how the AS events are correlated among datasets. For each SNP, the allele consistent with the reference genome is called the reference allele, and the other allele is called the non-reference allele.

After read mapping and data preprocessing (see Additional file 1: Supplemental Methods S.1), we count reads for each allele at each heterozygote SNP. For SNP i, dataset d and replicate sample j, let x _idj and y _idj be the read counts for the reference allele and non-reference allele respectively. Let n _idj =x _idj + y _idj be the total read count (See Figure 1a for a toy example). Protein-DNA binding can be skewed to the reference allele (SR), skewed to the non-reference allele (SN), or not allele-specific (NS). We use a binary variable b _id to indicate whether SNP i is SR (b _id =1) or not (b _id =0) in dataset d. If b _id =1, then SNP i is assumed to be SR in all replicate samples in dataset d. Similarly, we introduce another binary indicator c _id to indicate whether SNP i is SN or not in dataset d. b _id and c _id cannot be equal to one at the same time. If b _id =0 and c _id =0, then SNP i is NS in dataset d. The configuration at each SNP i can be described by two vectors B _i =(b_{i 1},⋯,b _iD ) and C _i =(c_{i 1},⋯,c _iD ) (See Figure 1d for a cartoon illustration). Based on these notations, (x _idj ,y _idj ), or equivalently (x _idj ,n _idj ), are the observed data for SNP i in sample (d,j), whereas the indicators b _id and c _id are unobserved.

Main intuition and challenge

Our primary goal is to infer for each SNP whether there is allelic imbalance in each dataset. This is equivalent to inferring b _id and c _id . A simple solution to this problem is to analyze each individual dataset separately, but this approach has low statistical power since the counts (x _idj ,n _idj ) usually are small.

If one knows how different datasets are correlated in terms of allelic imbalance, this knowledge may be used to improve the data analysis. For instance, if the allelic imbalance of two proteins A and B are closely correlated, then observing skewed read counts for protein A will provide information for inferring the allelic imbalance of protein B. Integrating the data from both A and B will increase the effective number of reads available for statistical inference, which will then lead to increased statistical power.

In reality, how different proteins are correlated is usually unknown. However, one may learn it by studying the data from many SNPs. Each SNP has three possible states in each dataset: SR, SN and NS. For D datasets, there are 3 ^D possible configurations in total. From studying many SNPs, one can know the relative frequencies (or mixing proportions) of these 3 ^D configurations. The mixing proportions will tell how different datasets are correlated. For instance, let [s₁,s₂,⋯,s _D ] be the skewness configuration of a SNP in the D datasets. If the mixing proportions for three configurations [NS,NS,⋯,NS], [SR,SR,⋯,SR] and [SN,SN,⋯,SN] are 0.9, 0.05 and 0.05, then no other configurations exist in the data and all datasets are perfectly correlated in terms of the allelic imbalance. In other words, at a particular SNP, if one dataset is SR, then all the other datasets are also SR. If one is SN, then all the others are also SN. On the other hand, if other configurations have non-zero mixing proportions, then not all datasets are perfectly correlated, and at a particular SNP, one allows the possibility that only a subset of datasets are correlated. For instance, if the mixing proportion for a configuration [SR,SR,NS,⋯,NS] is 0.03, then there will be 3% of SNPs that are skewed to the reference allele in the first two datasets but not skewed in the other datasets. Therefore, knowing the mixing proportions of all 3 ^D configurations will tell one the correlation structure in the data. This knowledge can then be used to improve statistical inference at each individual SNP by facilitating information sharing across datasets. For example, if the configuration [SR, SR, SN] has a much higher mixing proportion than [SR, SR, NS], then observing strong skewness towards the reference allele of a SNP in the first two datasets will imply that, a priori, the SNP is highly likely to be skewed to the non-reference allele in the third dataset and has much lower probability to be non-skewed for both alleles. The principle here is the same as the principle represented by the Bayesian hierarchical models in the statistical literature.

A limitation of this approach is that one has to enumerate all 3 ^D AS configurations in order to describe the correlation. As the number of datasets increases, the number of possible configurations increases exponentially. Thus this approach does not scale well with the increasing D. Later, in our analysis of GM12878 data, D=40 and 3 ^D >10¹⁹. This simple approach is clearly intractable.

To circumvent the difficulty of documenting the frequencies of all 3 ^D configurations, iASeq employs a technique that can describe the major correlation patterns in the data using a few probability vectors whose values vary from 0 to 1 rather than being dichotomous (i.e., 0 or 1). This approach significantly reduces the model complexity but keeps the flexibility to account for all 3 ^D configurations. It is easily scalable to increasing dataset number. The correlation structure in the model can then be used to improve the statistical inference of allelic imbalance at each SNP in each individual dataset.

Probability model

iASeq is based on the Bayesian hierarchical mixture model below that uses several probability vectors to describe the major correlation patterns among multiple datasets (Figure 1). The model assumes that SNPs can be grouped into K + 1 classes with different allele-specificity patterns (K≪3 ^D ), and the observed data are viewed as generated as follows:

● First, a class label a _i is randomly assigned to each SNP i according to a probability vector Π=(Π₀,Π₁,⋯,Π _K ). Here, Π _k =Pr(a _i =k) is the prior probability to assign a SNP to class k. $\sum _k{\Pi }_k=1$.

● If the class label a _i =0, then B _i =(0,⋯,0) and C _i =(0,⋯,0). In other words, all SNPs in class 0 are background SNPs, and they are NS in all datasets. If a _i =k and k≠0, then SNP i can be skewed, and its [b _id ;c _id ]s in different datasets are generated independently according to the following probabilities: Pr(b _id =1,c _id =0|a _i =k)=v _kd and Pr(b _id =0,c _id =1|a _i =k)=w _kd . We assume v _kd + w _kd <1, i.e., Pr(b _id =0,c _id =0|a _i =k)=1−v _kd −w _kd >0. The model implies that each class is associated with two vectors of probabilities V _k =(v_{k 1},⋯,v _kD ) and W _k =(w_{k 1},⋯,w _kD ). For SNPs in class k, B _i and C _i are generated according to the probabilities in V _k and W _k .

● Next, the observed read counts are generated based on the AS configurations specified by B _i s and C _i s. Consider SNP i and dataset d. If b _id =1, then (x _idj ,n _idj ) in each replicate sample (d,j) is generated according to a probability distribution Pr(x _idj ,n _idj |b _id =1,c _id =0)=Pr(n _idj |b _id =1,c _id =0)Pr(x _idj |n _idj ,b _id =1,c _id =0)≡Pr(n _idj )f_{idj 1}(x _idj ). Here we assume that the marginal distribution of n _idj does not depend on b _id and c _id , and we use f_{idj 1}(x _idj ) to denote the conditional distribution Pr(x _idj |n _idj ,b _id =1,c _id =0). Data in different replicate samples are assumed to be generated independently. Similarly, if c _id =1, then (x _idj ,n _idj )s are generated according to Pr(x _idj ,n _idj |b _id =0,c _id =1)=Pr(n _idj )f_{idj 2}(x _idj ). If b _id =0 and c _id =0, then (x _idj ,n _idj )s are generated according to Pr(x _idj ,n _idj |b _id =0,c _id =0)=Pr(n _idj )f_{idj 0}(x _idj ).

For SNP i and dataset d, we organize data from all replicates j=1,⋯,J _d into X _id =(x_{id 1},⋯,x _{idJ d} ) and N _id =(n_{id 1},⋯,n _{idJ d} ). For SNP i, X _i =(X_{i 1},⋯,X _iD ) and N _i =(N_{i 1},⋯,N _iD ) contain data from all datasets. The final observed data are X=(X₁,⋯,X _I ) and N=(N₁,⋯,N _I ) which are the ensemble of data from all SNPs.

Let A=(a₁,⋯,a _I ) be the collection of class membership indictors of all SNPs, and let B=(B₁,⋯,B _I ) and C=(C₁,⋯,C _I ) be the SR and SN indictors for all SNPs. A, B and C are the unobserved missing data one wants to infer.

Organize the probability vectors V _k and W _k from different classes into two matrices ${\mathit{V}}_{K\times D}={({\mathit{V}}_1^T,\cdots ,{\mathit{V}}_K^T)}^T$ and ${\mathit{W}}_{K\times D}={({\mathit{W}}_1^T,\cdots ,{\mathit{W}}_K^T)}^T$. V, W, and the probability vector Π that describes the class abundance are the unknown model parameters. K is assumed to be fixed. The choice of K and specification of data generating distributions Pr(n _idj ), f_{idj 0}(x _idj ), f_{idj 1}(x _idj ) and f_{idj 2}(x _idj ) will be discussed later.

Based on this model, each SNP class k (k≠0) is associated with two vectors of probabilities V _k and W _k which characterize the allelic imbalance preferences in different datasets for SNPs belonging to class k. For example, if a class has [V _k ;W _k ] = [(0.8,0.7,0.1,0.1); (0.1,0.1,0.8,0.1)], then SNPs in this class have high probability to be SR in datasets 1 and 2, and high probability to be SN in dataset 3, but they have low probability to be allele-specific in dataset 4. Since V _k and W _k are probabilities rather than 0-1 vectors, each class k can generate all 3 ^D AS configurations. Therefore, SNPs in the same class are not required to have the same AS configuration (e.g., a class can have one SNP with configuration [SR,SR,NS,NS] while at the same time another SNP with configuration [SR,NS,SR,NS]), although they usually have similar AS configurations because SNPs in the same class are all generated using the same probability vectors. Meanwhile, there are K different classes, and each class has a different [V _k ;W _k ] which specifies a different preference to generate the skewing configurations. Thus, whereas SNPs in the same class tend to have similar [B _i ;C _i ] configurations, SNPs from different classes tend to have very different configurations. Conceptually, this is similar to a model-based clustering analysis in which SNPs are grouped into K + 1 clusters based on their [B _i ;C _i ] configurations. However, an important difference here is that [B _i ;C _i ]s are unknown.

Our model assumes that [b _id ;c _id ]s of the same SNP in different datasets are a priori independent conditional on the class membership a _i . However, [b _id ;c _id ]s from different datasets are not independent marginally if one integrates out the class label a _i . For example, the marginal probability Pr([b _id ;c _id ]= [1;0])=$\sum _k\mathrm{Pr}\left(\right[b_{\mathrm{id}};c_{\mathrm{id}}]=[1;0\left]\right|a_i=k\left)\mathrm{Pr}\right(a_i=k)=$$\sum _{k=1}^K{\Pi }_k{v}_{\mathrm{kd}}$. On the other hand, the joint probability $\mathrm{Pr}\left(\right[{\mathit{B}}_i;{\mathit{C}}_i]=[(1,1,\cdots ,1);(0,0,\cdots ,0)\left]\right)=\sum _{k=1}^K{\Pi }_k\left(\prod _d{v}_{\mathrm{kd}}\right)$, which is clearly different from the product of the marginals $\prod _d\mathrm{Pr}\left(\right[b_{\mathrm{id}};c_{\mathrm{id}}]=[1;0\left]\right)=\prod _d\left(\sum _{k=1}^K{\Pi }_k{v}_{\mathrm{kd}}\right)$. This explains why our model can be used to describe the correlation among multiple datasets despite the conditional independence assumption. Intuitively, if one views the model as a clustering analysis of SNPs based on [B _i ;C _i ], then each cluster will represent a co-occurrence pattern of allele-specificity across multiple proteins. The marginal correlation among multiple datasets is described by multiple clusters, whereas within each cluster the data in different datasets are generated independently. In real data, a small K (i.e., a small number of SNP classes) usually is sufficient to describe the major correlation structure among datasets. Using Π, V and W to describe the correlation among datasets only requires O(KD) parameters, which is significantly less complex than O(3 ^D ) parameters. At the same time, the iASeq model still provides the flexibility to accommodate all 3 ^D possible [B _i ;C _i ] configurations as all of them have non-zero probability to occur.

Data generating distributions

To fully specify the model, one also needs to specify the data generating distributions Pr(x _idj ,n _idj |b _id ,c _id )=Pr(n _idj )Pr(x _idj |n _idj ,b _id ,c _id ). The primary goal of iASeq is to infer whether two alleles are different. We assume that information on allele-specificity is only contained in Pr(x _idj |n _idj ,b _id ,c _id ), and therefore the exact form of Pr(n _idj ), i.e., the marginal probability distribution of the total read count, is irrelevant for our purpose. As such, we mainly focus on modeling the conditional distribution of x _idj given n _idj , b _id and c _id , i.e., the three distributions f_{idj 0}(x), f_{idj 1}(x) and f_{idj 2}(x).

iASeq models these distributions hierarchically in two steps. First, x _idj is assumed to follow a binomial distribution x _idj |n _idj ,p _idj ∼Bin(n _idj ,p _idj ), where p _idj is the probability that a read generated at SNP i in sample (d,j) represents the reference allele. Next, we model p _idj depending on the values of b _id and c _id .

Here $C_n^k$ is the binomial coefficients “n choose k”, and B(. .) is the beta function.

Joint probabilities and model fitting

Conditional on the observed data, Pr(N) is a constant that does not contain parameters of interest, therefore it is absorbed into a proportionality constant not shown in the formula above. Using this joint posterior, an EM algorithm can be derived to search for posterior mode $(\widehat{\Pi },\widehat{\mathit{V}},\widehat{\mathit{W}})$ of $\mathrm{Pr}(\Pi ,\mathit{V},\mathit{W}|\mathit{X},\mathit{N})=\sum _{\mathit{A},\mathit{B},\mathit{C}}\mathrm{Pr}(\mathit{A},\mathit{B},\mathit{C},\Pi ,\mathit{V},\mathit{W}|\mathit{X},\mathit{N})$ in which the missing indictors A, B and C are all integrated out (see Additional file 1: Supplemental Method S.4).

For the Dirichlet prior, we use η=2 (see Additional file 1: Supplemental Method S.3 for a discussion on the choice of parameter for the Dirichlet prior). In the EM algorithm, we assume that the class number K is given. In order to choose the optimal K, we run the algorithm multiple times using different values of K. We choose the best K using the Bayesian Information Criterion (BIC) (see Additional file 1: Supplemental Method S.5).

Statistical inference of allele-specificity

Using ${\tilde{P}}_{\mathrm{id}}$, SNPs can be rank ordered for biologists to choose candidates to design follow-up studies. For each top ranked SNP, one can determine its skewing direction by comparing Pr(b _id =1,c _id =0|X _i ,N _i ,Π,V,W) and Pr(b _id =0,c _id =1|X _i ,N _i ,Π,V,W). The one with the larger value determines the direction. Finally, the posterior probabilities of top N SNPs can be converted to an estimate of false discovery rate (FDR) using $\mathrm{FDR}\left(N\right)=\sum _{i\in \mathrm{top}N\mathrm{SNPs}}(1-{\tilde{P}}_{\mathrm{id}})/N$.

Formula 7 shows that two types of information contribute to Pr(b _id ,c _id |X _i ,N _i ,Π,V,W): (1) Pr(a _i =k|X _i ,N _i ,Π,V,W), which is determined using information from all D datasets, and (2) Pr(b _id ,c _id |a _i =k,X _i ,N _i ,Π,V,W), which only uses information specific to dataset d conditional on Π, V and W. Thus for each particular dataset d, the dataset-specific information is weighted by information obtained from other datasets to determine the SNP ranking. Intuitively, if allelic imbalance in two datasets are correlated, then observing an AS event in one dataset will suggest that a relatively weak skewing event observed at the same SNP in the other dataset is very likely to be a true AS event. In contrast, if no AS event is observed in one dataset, then a relatively weak skewing event observed at the same SNP in the other dataset is likely to be a false positive. This is the underlying nature of using Pr(a _i =k|X _i ,N _i ,Π,V,W) to re-weigh information in Pr(b _id ,c _id |a _i =k,X _i ,N _i ,Π,V,W), and it provides the foundation for improving SNP ranking by borrowing information across datasets. In real applications, Π, V,W are unknown, and they are replaced by the posterior mode obtained from the EM algorithm.

GM12878 data and preprocessing

We collected 40 ENCODE [36] ChIP-seq datasets with a total of 77 samples together with a genomic DNA sample in GM12878 lymphoblastoid cells (Additional file 2: Table S1). GM12878 is a female and is one of the most extensively studied cell lines in ENCODE. Within each dataset, the number of replicate samples varied from 1 to 3. We downloaded the raw sequence reads of all 78 samples and mapped them to human genome (hg18) (see details in Additional file 1). We removed repeated sequences from the ChIP-seq datasets to avoid PCR duplicates, which may skew the determination of allelic biases. In other words, if multiple reads have exactly the same sequence, only one copy is retained. We obtained the genotype data for GM12878 from [37].

As previously described in [31], there are two different types of read mapping biases that may affect the analysis of AS events: the reference bias and the inherent bias. The reference bias often occurs when one maps sequence reads to a reference genome. If one allows the same number of mismatches in the alignment, a read from the non-reference allele is less likely to be mapped back to the reference genome compared to a read from the reference allele, since the non-reference read has one more mismatch to the reference genome. This phenomenon is known as the reference bias. This type of bias, if it exists, is automatically taken care of by the iASeq model through the parameter p_{dj 0} which models the background skewing probability and is estimated using all reads mapped to heterozygote SNPs in each sample. If there is reference mapping bias, p_{dj 0}will take a value different from 0.5 to adjust for the bias. One may remove reference bias before the analysis by masking SNPs in the reference genome during the alignment or by aligning reads to a diploid personal genome. This situation will also be automatically recognized by iASeq through the estimation of p_{dj 0} from the data (if there is no bias, p_{dj 0}=0. 5). Therefore, regardless of whether the reference bias has been removed from the data in the preprocessing or not, the iASeq model is able to automatically handle it and adjust the inference accordingly.

The intrinsic bias is a different type of bias. As shown by [31], even if the reference bias is removed (e.g., by masking SNPs in the reference genome), the inherent bias still exists. For example, suppose sequence 1 (e.g., xxxAxxx) and sequence 2(e.g., xxxTxxx) are two reads that differ only in one position (i.e., A/T). It is possible that sequence 1 is easier to be mapped back to its correct location in the genome than sequence 2 if the second sequence has many repeats in the genome. This bias reflects the inherent characteristics of the genome and cannot be removed by masking variants in the reference genome or by mapping reads to a diploid personal genome. In the above example, masking A and T in the original reads is also not a solution, since a priori one does not know which position in a read corresponds to a SNP position and therefore should be masked without first aligning the read to the genome. When a heterozygote SNP has inherent bias, one allele will have higher read counts than the other even if the two alleles have the same binding level. To avoid this bias, we used the approach described in [22, 31] to remove SNPs with the inherent bias.

We began with 1,704,166 heterozygote SNPs and filtered out 149,996 (8.8%) SNPs with inherent bias. Next, we eliminated SNPs that were not bound by any TF or associated with any HM in any dataset (see Additional file 1: Supplemental Methods S.1.1, S.1.2 and Additional file 3: Table S2 for details). After applying these filters, 94,519 heterozygote SNPs remained. These 94,519 SNPs were then analyzed by iASeq.

A simulation study

Before we apply iASeq to the real data, we first tested its performance in simulations that took into account real data characteristics. Our simulations kept the same design as the real GM12878 ChIP-seq data, with the same number of datasets and the same number of replicates within each dataset, except that the genomic DNA sample was not used here since we knew the truth in the simulations and did not need genomic DNA as a control for potential bias. To create the simulation data, we first applied iASeq to the real GM12878 data to identify 86,353 SNPs that were not skewed in any dataset using Pr(a _i =0|X _i ,N _i ,Π,V,W)>0. 5 as cutoff. To mimic the real background noise, these SNPs were resampled by a bootstrap procedure to create the background SNPs in the simulations, and we kept the read counts (x _idj ,n _idj ) of each background SNP as is in the simulated data. Next, we simulated ASB SNPs and added them to the background. Simulations were carried out under two different scenarios (Figure 2).

● Scenario 1: Two types of ASB SNPs (classes 1 and 2) were created in addition to the background SNPs (class 0). The SNP number for class 0, 1, and 2 was 85,069, 4,725 and 4,725 respectively. Thus the true Π _k for the three classes was 0.90, 0.05 and 0.05 respectively. SNPs in class 1 were SR in datasets 1 to 30 (i.e., their b _id =1 for d=1,⋯,30). SNPs in class 2 were SN in datasets 1 to 30 (i.e., c _id =1 for d=1,⋯,30). In datasets 31 to 40, no SNPs had allelic imbalance. Class 2 can be viewed as the mirror image of class 1. This symmetric design reflects the symmetry of allele-specificity, that is, the skewing to the reference allele and to the non-reference allele is approximately symmetric. The class abundance (0.90,0.05,0.05) roughly matched the abundance observed in the analysis of real GM12878 data.

● Scenario 2: Four correlation patterns (classes 1-4) were created in addition to the background class (class 0). Class 1 and class 2 were the same as in simulation 1. Classes 3 and 4 were two new patterns. SNPs in class 3 were SR in datasets 21-40, and SN in datasets 1-10. Class 4 was the mirror image of class 3. The abundance of the classes 0 to 4 was (0.90,0.025,0.025,0.025,0.025).

Given the simulated [B _i ;C _i ] configurations, we then simulated the read count data for ASB SNPs as described in detail in Additional file 1: Supplemental Methods S.6. Simulations done in this way was able to keep the major characteristics of real data while allowing us to benchmark the performance of different methods since we knew the truth.

We applied iASeq to both simulations. In both cases, iASeq was able to identify the correct number of SNP classes using BIC (Figures 2a,b,d,e). Figures 2c and 2f show that the ASB patterns reported by iASeq matched the true patterns well. In order to test whether iASeq can improve the statistical power of detecting SNPs with allelic imbalance, we compared the SNP ranking provided by iASeq with rankings provided by five other methods that analyze each dataset separately (Figure 3). In iASeq, SNPs were ranked in each dataset d according to the posterior probability ${\tilde{P}}_{\mathrm{id}}$ defined by Formula 8. Since we know the truth, we can count how many of the top N SNPs were true positives. Here the true positives were defined as SNPs that were truly allele-specific and also had the skewing direction correctly inferred. The five single-dataset based methods for ranking SNPs include a deviation statistic d, naive z statistic, naive Bayes statistic, empirical Bayes statistic and single dataset EM. These methods were applied to each individual dataset. For each dataset d, we merged data from all replicates to obtain $x_{\mathrm{id}}=\sum _{j=1}^{J_d}{x}_{\mathrm{idj}}$ and $n_{\mathrm{id}}=\sum _{j=1}^{J_d}{n}_{\mathrm{idj}}$. We then computed the statistics used for SNP ranking as described below.

1. 1.

   Deviation statistic (d): SNPs were ranked based on |x _id /n _id −p _{d 0}|. Here we estimated $p_{d0}=\frac{1}{I^{\prime }}\sum _{i:n_{\mathrm{id}}\ne 0}{p}_{\mathrm{id}}=\frac{1}{I^{\prime }}\sum _{i:n_{\mathrm{id}}\ne 0}\frac{x_{\mathrm{id}}}{n_{\mathrm{id}}}$, where I ^′ is the number of SNPs for which n _id ≠0.

    

Figure 3 compares the number of true positives, T P _d (q), in the top q SNPs reported by each method in each dataset d. In Figures 3a-c and 3e-g, T P _d (q) is plotted as a function of q in a few representative datasets. These plots show that iASeq outperformed all single-dataset based methods, and it was able to substantially improve the power for detecting allele-specificity.

In general, the observed differences between iASeq and the d, z, b and B statistics could be caused by many factors such as use of different statistical models, ranking statistics, or methods for parameter estimation. However, the comparison between iASeq and the single dataset EM represents a well-controlled comparison since these two methods used exactly the same distributional assumptions and parameter estimation methods. The only difference between them was that iASeq used information from multiple datasets whereas singleEM was based on one dataset only. This well-controlled comparison shows that jointly modeling multiple datasets is able to improve the allelic inference.

To examine whether iASeq was able to bring improvement in all datasets, we computed the Area under the Receiver Operating Characteristic (ROC) curves (AUC) for each method in each dataset using the top 2000 ranked SNPs. In each dataset, we computed the proportion of improvement in terms of AUC brought by iASeq over the best single-dataset based ranking method (i.e., $\mathrm{dAUC}=\frac{\mathrm{AU}{C}_{\mathit{\text{iAseq}}}-\mathrm{AU}{C}_{\mathit{\text{bestsingle}}}}{\mathrm{AU}{C}_{\mathit{\text{bestsingle}}}}$). dAUC>0 means iASeq is able to bring improvement. Figures 3d and 3h show the distribution of dAUC across all 40 datasets as a histogram. The results show that iASeq was able to improve the SNP ranking in almost all datasets.

In Figure 4, we converted the iASeq posterior probabilities of top N SNPs to FDR estimates and plotted the estimated FDR against the true FDR. The figure shows that iASeq was able to provide reasonable FDR estimates as well. Shown in the figure are a few representative datasets. Results in all other datasets were similar.

Analysis of real data

Our simulation study demonstrates the ability of iASeq to discover correlation patterns of allele-specificity and improve the detection of skewed SNPs. Next, we applied iASeq to analyze the 41 real datasets (78 samples) in GM12878 cells. In real data, we do not have comprehensive knowledge about the truth. Therefore, unlike simulations, we were not able to assess the FDR estimates. For this reason, we mainly focused on analyzing the correlation patterns of allele-specificity and testing whether iASeq can improve the SNP ranking.

Correlation patterns of allele-specificity

Figure 5a shows the BIC in the real data. Based on BIC, the optimal K was 2. In other words, in addition to the background class (k=0), iASeq discovered two other SNP classes, representing different allele-specificity patterns. For these two non-background classes, Π _k was estimated to be 0.0696 and 0.0691 respectively, suggesting that they cover 6.96% and 6.91% of the analyzed SNPs. Due to the background noises, not all SNPs in these two classes can be confidently detected. At the 0.90 posterior probability cutoff, iASeq reported 1868 and 2138 SNPs for classes 1 and 2 respectively (Figure 5b). Note that our simulations had similar settings as the real data analysis, and they showed that iASeq was able to discover more than two patterns if they are supported by the data. Therefore our discovery of two correlation patterns here is likely driven by the data, that is, the information in the data is only sufficient for supporting robust discovery of two patterns.

Figures 5b and 5c show the posterior mode of V _k and W _k for the two non-background classes. It turned out that these two classes corresponded to two global directions of allele-specificity, SR and SN, respectively. Since the assignment of reference or non-reference allele depends on the reference genome, the assignment per se is not of biological interest. However, recall that GM12878 is a single person, therefore at each single SNP, the nucleotide representing the reference or non-reference allele is the same across all datasets analyzed here. Given this fact, what these results essentially tell is that at each single SNP, most TFs and HMs in our analysis were highly correlated in terms of allele-specificity, and if they are skewed, they tend to be skewed toward the same direction (i.e., the same allele). For instance, for SNPs in class 1, both H3K4me3 (from the Broad Institute) and H3K27ac (Broad) had high probability to be SR, with (v _kd ,w _kd ) equal to (0.9337,0.0070) and (0.9730,0.0041) respectively (Figure 5c). The probability that one is SR and the other one is SN was small. Similarly, for SNPs in class 2, both H3K4me3 and H3K27ac were highly likely to be SN simultaneously ((v _kd ,w _kd )=(0.0061,0.9835) for H3K4me3 (Broad) and (0.0040,0.9897) for H3K27ac (Broad)). While the allelic imbalance of most TFs and HMs were highly correlated, H3K27me3, a HM involved in gene repression, was an exception. In both non-background classes, H3K27me3 had much lower skewing probabilities compared to the other proteins (Figure 5c). Within each class, the difference in the skewing probability between the two alleles was also much weaker for H3K27me3 as compared to the other proteins. For instance, in class 1, while most other proteins showed strong preference to be skewed toward the reference allele, H3K27me3 can be skewed to the reference allele at some SNPs and skewed to the non-reference allele at many other SNPs. Therefore, the allelic imbalance in H3K27me3 is not strongly correlated with the allelic imbalance of the other proteins analyzed here. For the genomic DNA which was used as control here, the skewing probabilities (v _kd ,w _kd ) in both classes were fairly low as shown in Figure 5b-c. In both classes, the probability for not being skewed in the genomic DNA (i.e., 1−v _kd −w _kd ) was bigger than 0.95. This indicates that the high probability of skewing observed in the other datasets was not an artifact.

The coordinated allelic imbalance of different proteins toward the same allele has also been observed in a recent study [38]. In that study, the authors analyzed AS of 24 TFs and found that when multiple TFs bind to the same SNP, they frequently bind to the same allele. Moreover, those authors did not observe any pair of TFs that regularly bind the same position on alternate alleles. Our observation here therefore is consistent with their finding.

Increased power for detecting allele-specificity compared with single dataset analysis

We ranked SNPs based on the posterior probabilities ${\tilde{P}}_{\mathrm{id}}$ in each dataset. The iASeq ranking was compared with the rankings provided by the five single-dataset based methods described above. Since we do not know the truth, we used two types of independent information as gold standard to benchmark the ranking results.

First, we evaluated different methods by counting how many of their top ranked SNPs were located in the non-pseudoautosomal regions of chromosome X (chrX-npa) (Figure 6). GM12878 is a female lymphoblastoid cell line. In GM12878, SNPs in chrX-npa are expected to be allele-specific due to cells rapidly become clonal in culture leading to a skewed X-inactivation [5, 38, 39]. Therefore, given a fixed number of top SNPs, the more chrX-npa SNPs one can find, the more powerful a method is. Figure 6 shows that iASeq clearly increased the power for detecting allele-specificity in each dataset compared to the single-dataset based analysis. For example, Figure 6a shows that in the H3K27ac dataset generated by the Broad Institute, iASeq was able to identify 122 chrX-npa SNPs among the top 500 SNPs. This represents 126% improvement compared to singleEM, the best single-dataset based ranking method in that dataset, which only identified 54 chrX-npa SNPs. Figures 6a-g show results in a few representative datasets. Figure 6h shows the distribution of dAUC (i.e., the proportion of improvement of AUC by iASeq over the best single-dataset based ranking method in each dataset) in all 40 datasets. These plots clearly show that iASeq outperformed all single-dataset based methods in all datasets and the average improvement in AUC was 354%.

Second, we evaluated different methods by using independent RNA-seq data. From RNA-seq, one can identify exonic ASE SNPs and use them as gold standard. We collected two RNA-seq datasets in GM12878, one from the California Institute of Technology (Caltech) and the other from the Yale/Stanford University (Yale) (Additional file 2: Table S1). From each dataset, we identified the top 400 exonic ASE SNPs using the naive Bayes statistics. Using the other methods to identify the gold standard ASE SNPs produced similar results which, for simplicity, will not be shown here. Based on these exonic ASE SNPs, we defined a SNP in our ChIP-seq analysis as truly allele-specific if there was an exonic ASE SNP in its X kb neighborhood. Here we tried both X=10kb and X=1kb and obtained similar results. Below we illustrate the results using X=10kb as an example. The corresponding results for X=1kb are shown in Additional file 4: Figures S3-S6. Among the 94,519 SNPs analyzed in the ChIP-seq data, 20,526 had one or more exonic SNPs within its 10kb neighborhood and therefore could potentially be linked to an exonic ASE SNP. Figure 7 and Additional file 5: Figure S1 compare rankings of these SNPs provided by different methods in terms of how many of the top ranked SNPs are true positives (i.e., associated with ASE). iASeq again outperformed all the other single-dataset based ranking methods. For instance, based on the Caltech gold standard, iASeq on average identified 144% more true positive SNPs among the top 500 SNPs (Figure 7a-g). According to the Yale gold standard, iASeq achieved an average of 149% improvement in terms of the true positive rate among top 500 SNPs (Additional file 5: Figure S1). The average improvement in terms of AUC (i.e., dAUC) across all 40 datasets was 148% (Figure 7h) and 165% (Additional file 5: Figure S1h) for the Caltech and Yale gold standard respectively.

To ensure that the increased statistical power was not completely attributed to X chromosome SNPs, we repeated the benchmark analysis based on RNA-seq using only SNPs in autosomal chromosomes, and we obtained similar results (Figure 8, Additional file 6: Figure S2). This shows that the increased power is not only contributed by chrX SNPs.

Interpretation of the correlation patterns

When analyzing the real data in GM12878, iASeq found two non-background AS patterns, representing two opposite directions of allelic imbalance. Since the assignment of reference and non-reference allele depends on the reference genome, whether a SNP is skewed toward reference or non-reference allele per se does not have direct biological meaning. What these two patterns essentially suggest is that the allelic imbalances of multiple proteins at a single SNP are correlated and have high preference to be skewed toward the same allele. In other words, the two patterns should be viewed as a pair and interpreted together.

In general, although one may view different allelic imbalance patterns in iASeq as different clusters of SNPs, these clusters only describe the similarities among SNPs in terms of their skewness directions, rather than the similarities in terms of their functions. The direction is defined using the reference/non-reference allele. The reference or non-reference allele for different SNPs can have different meanings (e.g., for one SNP, the maternal allele may be the reference allele, whereas for another SNP the paternal allele may be the reference allele). Therefore within each cluster, even though SNPs have similar skewness pattern, they are not necessarily functionally related to each other. One should not confuse the SNP clusters here with the clusters obtained from the traditional gene expression microarray data analysis, where co-expressed genes in a cluster often have similar functions. In iASeq, the clusters only serve as a tool to describe the correlation structure among different datasets (i.e., proteins), rather than the functional correlation among different SNPs. The correlation patterns among datasets are used by iASeq to inform one how to integrate information across datasets (i.e., which datasets are highly correlated and therefore can borrow information from each other) to improve detection of AS events for each individual SNP and dataset. In order to understand functions of the detected AS events, one needs to further correlate the iASeq results with other information (e.g., one may determine the parent-of-origin of each SNP first and then study various phenomena such as imprinting).

Our observation that different proteins prefer to be skewed in the same direction is consistent with a recent observation reported in [38] that AS of 24 different TFs are skewed toward the same allele. A number of factors could contribute to the observed correlation. First, biologically it is plausible that functionally related HMs and TFs have correlated allele-specificity. For instance, both H3K4me2 and H3K4me3 are markers for active transcription. Therefore, for a specific SNP, if the reference allele is associated with a gene with active transcription but the non-reference allele is not, then it is very likely that both H3K4me2 and H3K4me3 will be skewed toward the reference allele. For another SNP, if the non-reference allele is transcribed but the reference allele is not, then both H3K4me2 and H3K4me3 will have high probability to be skewed toward the non-reference allele. In the genome, H3K4me2 and H3K4me3 are skewed toward reference allele for some SNPs, and skewed toward non-reference allele for some other SNPs. Therefore the skewed SNPs could naturally fall into two clusters, representing two opposite AS directions. Second, as pointed out by [38], the coordinated AS could also occur as a result of the difference in the chromatin landscape between the two alleles. For instance, if the chromatin on one allele is more open and accessible, it could increase the overall binding probability of multiple different proteins, leading to correlated allelic skewing.

While our results show that most analyzed TFs/HMs tend to be skewed toward the same direction, these results do not imply that these proteins are perfectly correlated in terms of allele-specificity at each and every SNP. In iASeq, the correlation patterns V _k and W _k are probabilistic patterns rather than 0-1 vectors. Each correlation class k can generate all 3 ^D AS configurations. For instance, for a class with [V _k ;W _k ]= [(0.9,0.9,0.9,0.1); (0.1,0.1,0.1,0.1)], it is possible to have one SNP with configuration [SR,SR,NS,NS] and at the same time another SNP with configuration [SR,NS,SR,NS]. Therefore, SNPs in the same class are not required to have the same AS configuration, even though they tend to have similar AS configurations. The probabilistic patterns are used here to provide a parsimonious description of the complex correlation structure in the data, so that one can circumvent the difficulty of handling 3 ^D AS configurations whose complexity increases exponentially. As a consequence of using this parsimonious model, multiple weak correlation patterns without strong enough data support could be merged into a bigger class. For instance, consider two AS patterns [V _k ;W _k ]=[(1,1,0,0);(0,0,0,0)] (i.e., [SR,SR,NS,NS]) and [V _k ;W _k ]=[(0,0,1,1);(0,0,0,0)] (i.e., [NS,NS,SR,SR]). Suppose both patterns are equally likely to occur in the data. If each pattern is only associated with a small number of SNPs, then a parsimonious model will prefer merging them together into one single class with [V _k ;W _k ]=[(0. 5,0. 5,0. 5,0. 5);(0,0,0,0)]. For this reason, iASeq only discovers correlation patterns that have sufficient data support so that they can be distinguished from other patterns. It will not report weak patterns, which could be real but do not have enough data support to allow them to be robustly recovered. For users, this means that at the cluster level, they may not be able to see weak but real AS correlation patterns if these patterns are not associated with enough number of SNPs. On the other hand, for the purpose of inferring whether or not each SNP is allele-specific in each dataset, these parsimonious correlation patterns are sufficient for describing the correlation structure in the data and serving as a prior to guide the information sharing across datasets. The information sharing will lead the increased ASB detection power, and the eventual AS configuration at each individual SNP will be determined by the posterior probabilities of (b _id ,c _id ) (i.e., ${\tilde{P}}_{\mathrm{id}}$) rather than the cluster-level prior probabilities [V _k ;W _k ]. Therefore, in the final AS calls, the model still allows each SNP to have its own AS configuration which may not necessarily be the same as the AS configurations of other SNPs from the same cluster.

Consistent with [38], in the two non-background AS patterns discovered here, proteins skewed toward the same direction did not always correspond to known protein-protein interactions. As pointed out by [38], this could happen as a result of allelic imbalances of different proteins being caused by a common underlying factor such as allelic difference in chromatin landscape. It could also reflect unknown protein-protein interactions. For iASeq specifically, there is a third reason, that is, multiple small patterns can be merged into a bigger probabilistic class as described before. For example, because of the use of probabilistic patterns, two patterns [SR,SR,NS,NS] and [NS,NS,SR,SR] may be merged into a single SNP class (e.g., V _k ;W _k =[(0. 5,0. 5,0. 5,0. 5);(0,0,0,0)]). As a result, only looking at the pattern represented by V _k ;W _k , one cannot tell the details of protein-protein interactions, such as these interactions only exist between datasets 1 and 2, or between 3 and 4, but not between the other pairs of datasets. What one can tell from this merged pattern is that, when the allelic imbalance occurs in these four datasets, they will be skewed toward the same direction, i.e., the reference allele in this example.

In summary, while the correlation patterns in iASeq provide some insights on the correlation of allelic imbalance among different datasets, one should not over-interpret them. The primary goal of these patterns is to describe the correlation structure in the data so that information from different datasets can be shared in a principled way to increase the power of statistical inference. This also points to an important difference between this study and previous studies that reported coordinated allele-specificity among multiple proteins. The previous studies only reported the correlation as a biological finding, but did not provide a statistical method to further utilize the correlation structure to improve the statistical inference. In contrast, iASeq provides a general and rigorous statistical method that utilizes the automatically discovered correlation patterns to increase the statistical power of AS detection. As such, it represents a novel development for the analysis of allele-specificity.

Model, algorithm, and possible extensions

Unlike tools such as AlleleSeq which mainly focus on the preprocessing steps for the AS analysis (e.g., construction of diploid personal genome), iASeq is developed as a general model working downstream of the preprocessing pipelines. The input data for iASeq are the read counts in the format shown in Figure 1a. With this design, iASeq can be easily coupled with different data preprocessing protocols. For instance, some investigators may map their reads to a reference genome, while others may map their reads to a diploid personal genome. Both types of investigators can use iASeq to integrate information from multiple datasets once they obtained the allelic read counts.

In iASeq, we used an EM algorithm to find the posterior mode of parameters and carried out statistical inference accordingly. In principle, one may also use a full Bayesian approach and Markov Chain Monte Carlo (MCMC) to perform the posterior inference. However, since MCMC usually takes much longer to run for a big dataset and it is not easy for users to monitor convergence, we decided to use the posterior mode and EM-based approach in our implementation. For analyzing the GM12878 data with 94,519 SNPs, iASeq took 5 hours to run the EM algorithm to fit a single model with K=1 on a machine with 2.7 GHz CPU and 4Gb RAM. To fit a single model with K=10 on the same machine, the EM took 16 hours. Running the EM for all 10 K s between 1 and 10 on a single core took 4.6 days. However, when we run these 10 jobs in parallel on 10 cluster nodes, we were able to select the best model within 1 day. Therefore, running the algorithm on a single machine is a little time-consuming, but the computation time can be reduced by parallelization. Also, our analysis of GM12878 data indicates that the optimal K in that real data was 2. For a K not extremely large, even if running the full BIC selection on a single machine takes some time, it usually requires less than a week, which is acceptable compared to the time devoted to preparing samples and generating data.

In principle, the statistical model developed in iASeq may also be applied to analyze other types of AS events, such as ASE and ASM. In the future, we plan to improve the model by incorporating information from the spatial correlation among closely located SNPs. For example, for the ASE analysis, one may jointly model SNPs from the same gene, similar to [33].

Implications on future studies

The analysis of AS events using the high-throughput sequencing data frequently faces the problem of low statistical power due to the limited amount of information available at heterozygote SNPs. One way to increase the power is to increase the sequencing depth for one data type (e.g., MYC ChIP-seq). An alternative approach is to spend the same amount of money to generate data for multiple different but related data types (e.g., ChIP-seq for MYC, H3K4me1, H3K4me3, etc.), each with a lower coverage. One can then integrate the multiple datasets to increase the statistical power of allele-specificity analysis. The merit of the second approach is that one can collect multiple different types of information which might be useful for other purposes (e.g., in addition to studying MYC binding using MYC ChIP-seq, one may couple H3K4me1 ChIP-seq data with DNA motif information to locate active enhancers and predict binding sites of other TFs in the genome). If the second approach is used in the study design, then iASeq will offer a flexible, powerful and scalable framework for better analyzing the AS events in the data. As ChIP-seq data continue to grow rapidly, this integrative approach will allow us to use the data more efficiently to characterize the allele-specificity.