BACKGROUND
diploid eukaryotic species inherit two copies  of the same gene from both parents.
if one allele fails to work properly, the other allele can still implement a gene’s cellular function.
for some genes, however, this protective mechanism is disabled because only one allele is expressed and its failure probably leads to cellular malfunction.
these monoallelically expressed genes can be classified into one of three categories  <cit> : x-inactivated genes, which are regulated by a random process in which one of the two x chromosomes present in female mammals is silenced  <cit> ; autosomal genes subject to random monoallelic gene expression, such as the t cell receptors and natural killer cell receptors  <cit> ; and autosomal imprinted genes , which express from only one of the two alleles according to their parental origin  <cit> .
imprinted genes play important functional roles in the control of embryonic growth and development, as well as in post-natal development  <cit> .
as imprinted genes are expressed from only one of the two parental chromosomes, a de facto haploid state is caused by imprinting and leads to asymmetric functions of parental genomes and loss of diploid protection against recessive mutations  <cit> .
thus, imprinting dysregulation is linked to numerous human genetic diseases, such as developmental disorders  and cancers   <cit> .
in addition, environmental factors can influence gene expression by targeting imprinted genes  <cit> .
because imprinted genes are more susceptible to disease than non-imprinted genes  <cit> , it is of great importance to identify novel imprinted genes for human.
identifying imprinted genes experimentally has traditionally been a slow process, and the number of validated ones is much lower than the previous estimation   <cit> .
however, recently developed high-throughput screening approaches  and recently identified dna sequence characteristics  have led to the proposal of several new methods to predict imprinted genes on a global scale  <cit> .
with advances in next-generation sequencing technology  <cit> , mrna-seq is becoming a powerful tool for transcriptome profiling  <cit> .
it can generate not only the number of reads mapped to exons, which reflects the expression levels of a gene, but also the actual sequence, which may identify the allele from which the mrna is expressed.
therefore, inference can be made for predicting imprinted genes.
for example, by sequencing whole transcriptomes from mice embryos, babak et al.
measured allelic expression bias and identified six novel imprinted genes  <cit> .
however, to our knowledge, prediction of imprinted genes by deeply sequencing transcriptomes  from multiple independent tissues is still an open problem.
in this study, we proposed a bayesian model – dspig  – to predict imprinted genes based on allelic expression inferred from observed snps in mrna-seq data of independent human tissues.
with dspig, we were able to measure the imbalance of allelic expression among various tissues and calculate the posterior probability of imprinting status for each gene.
under a stringent probabilistic cut-off of the posteriors and other reasonable biological criteria, we identified  <dig> potentially imprinted genes from  <dig> diverse human tissues and  <dig> potentially imprinted genes from  <dig> cerebellar cortices, with an expected low false discovery rate .
furthermore, analysis of allelic expression of the same genes among different tissues revealed that, in some cases, even if a gene was biallelically expressed, one allele always had higher expression level than the other.
RESULTS
statistical model development
monoallelic expression generally falls into one of three categories: imprinted expression, random monoallelic expression and x-inactivation, all of which express only one of two alleles in a single cell  <cit> .
at a tissue level, however, random monoallelic expression will allow both alleles to be detected in total rna because of the “mosaicism” of the tissue  <cit>  .
because our study was based on whole transcriptomes of tissue samples, random monoallelic expression was reasonably considered as biallelic expression when averaged over the entire tissue.
x-inactivation was also excluded from this study by discarding all predictions on the x chromosome.
thus imprinting is the most likely cause of the observed monoallelic expression among transcriptomes of different tissues even though we cannot infer the parent of origin.
we used known snps from dbsnp  <cit>  to distinguish and count the two alleles of each gene.
if a gene was imprinted, we expected to observe only one of the two alleles of each snp in the exon region from the whole transcriptome.
with the allelic counts obtained from the mrna-seq data , we developed a bayesian model  to compute the posterior probability of imprinting based on each single snp.
suppose we have sequenced transcriptomes from n independent tissue samples.
for each sample, we count the alleles of all known snps, discarding those with  <dig> counts.
for each snp, let the allelic counts be: , …, , where xk and yk are the counts for two alleles x and y in the sample k .
because each gene can only have two statuses: imprinted  or non-imprinted , we consider  as a binary event vector for the imprinting status.
by denoting data = {, …, }, we have by bayes’ theorem:
pri|data=prdata|i×priprdata|i×pri+prdata|ni×prni
where pr is the posterior probability of imprinting and pr is the prior of imprinting.
based on current knowledge of prevalence of imprinted genes  <cit> , we set the prior pr=1%, thus pr =  <dig> – pr = 99%.
since samples were independent of each other and genotype has only  <dig> possible combinations , we denote the genotype as follows:
δ={1the genotype isxx2the genotype isxy3the genotype isyy
assuming p and q are the allele frequencies for allele x and y, p + q = <dig>  according to hardy-weinberg equilibrium, the prior probabilities for the three genotypes are calculated as follows:
prδ={=prx2=p2δ=1=prpr=2pqδ=2=pry2=q2δ= <dig>
since values of p and q can be retrieved from dbsnp  <cit> , p and q are treated as constants.
we used the law of total probability to calculate the likelihood pr as follows:
prdata|i=∏k=1nprxk,yk|in≥ <dig>
=∏k=1nprxk,yk|i,δ=1×p2+prxk,yk|i,δ=2×2pq+prxk,yk|i,δ=3×q <dig>
by assuming that  the transcript levels of a gene’s two alleles are equal if the gene is biallelically expressed,  two different alleles, if both expressed, have the same chance to be sequenced by mrna-seq, and  unexpressed alleles with counts > <dig> are due to sequencing errors, we had the following derivation:
prxk,yk|i,δ=1=fyk;nk,pe
prxk,yk|i,δ=3=fxk;nk,pe
prxk,yk|i,δ= <dig>
=prxk,yk|i,δ= <dig> θ=0×prθ=0+prxk,yk|i,δ= <dig> θ=1×prθ= <dig>
=fyk;nk,pe×12+fxk;nk,pe× <dig>
here, f denotes binomial distribution , nk = xk + yk is assumed fixed, and pe is the averaged sequencing error rate for each snp .
the binary variable θ is defined as follows θ={0onlyxcan be expressed due to imprinting1onlyycan be expressed due to imprinting.
since x and y have an equal chance to be inherited from either maternal or paternal genome, x and y have an equal chance to be expressed in imprinted genes.
hence, prθ={1/2θ=01/2θ= <dig>
similarly for the likelihood function pr, we have:
prdata|ni=∏k=1nprxk,yk|nin≥ <dig>
=∏k=1nprxk,yk|ni,δ=1×p2+prxk,yk|ni,δ=2×2pq+prxk,yk|ni,δ=3×q <dig>
prxk,yk|ni,δ=1=fyk;nk,pe
prxk,yk|ni,δ=3=fxk;nk,pe
based on our three assumptions  – , pr follows a binomial distribution with p= <dig> .
therefore, we have:
prxk,yk|ni,δ=2=xk+ykxk12xk+yk
computation is performed separately for each single snp.
therefore, a posterior probability of imprinting for a gene is associated with a specific snp in this gene, unless otherwise specified.
simulation-based model analysis
we generated allelic counts from simulated data by taking into account imprinting status, snp frequency, and sequencing error.
we then applied dspig to estimate the sensitivity, specificity and the fdr.
we generated two sets of allelic counts under the assumption that the locus was either imprinted or not.
the number of reads for each allele was generated assuming the presence of one dominating allele plus sequencing error for imprinted case; and presence of equal amount of two alleles plus sequencing error for non-imprinted case.
the generated allelic counts followed a distribution similar to the actual mrna-seq data .
more details of the procedure were illustrated in additional file  <dig> figure s <dig>  given an allele frequency , a sequencing error  and a prior  of imprinting, the posterior probability calculated from the allelic counts generated for imprinted genes approached  <dig> as the sample size increased, while concomitantly the posterior probability for biallelically expressed approached  <dig> .
with minor allele frequencies between  <dig>  to  <dig> , sequencing errors between  <dig>  to  <dig>  and priors between  <dig>  to  <dig> , we obtained similar results as figure 1a,b.
using  <dig>  as the cut-off for posteriors pr , sensitivity of our model-based prediction exceeded  <dig> % when sample size was > <dig>  and specificity exceeded  <dig> % when sample size was > <dig> .
under the same cut-off and the allele frequency, the fdr of predicted imprinted genes decreased as sample size increased: it dropped to 5% and 1% as sample size exceeded  <dig> and  <dig>  respectively .
using the same cut-off , we also examined how fdr varied when minor allele frequency changed from  <dig>  to  <dig>  and sample size increased from  <dig> to  <dig> .
based on these simulations, we were able to provide estimated fdrs for most of our predictions in this study .
predictions of imprinted genes
we collected two previously published mrna-seq data sets.
one set included  <dig> diverse tissue samples  from wang et al.
<dig>  and the other set had  <dig> cerebellum cortex samples  from mudge et al.
<dig>   <cit> .
wang et al.
<dig> showed that, in terms of alternative isoform expression, cerebellum tissues were clustered together and the  <dig> diverse samples were more closely correlated.
here we performed hierarchical clustering based on the imprinting-inclined snps  and obtained similar results .
as we expected, samples from the cerebellar cortex were clustered together, with caucasian and african american separated in two sub-clusters .
using caucasian allele frequency on african american samples, however, yielded a sub-cluster without separation between the two ethnicities .
this suggests that the separation observed in figure 3a was due to differences in minor allele frequencies.
as a test set,  <dig> breast cancer cell line samples were clustered together in both cases.
compared with other non-cerebellum samples, the brain sample had higher correlation with cerebellum samples in both cases, which is sensible biologically.
the result that the  <dig> diverse tissue samples were clustered together could be caused by many factors such as different experimental conditions between group i and group ii samples.
using dspig, we predicted imprinted genes for group i and group ii separately.
to call a gene imprinted, we checked the posteriors of all the snps on the same gene to make sure that  the gene had at least one snp with a posterior > <dig> , which was the same cut-off used in simulations, and  all the other snps did not show any contradictory evidence .
after applying dspig to the mrna-seq data and using the above criteria, we predicted  <dig> potentially imprinted genes for group i samples and  <dig> potentially imprinted genes for group ii samples out of a total of  <dig>  genes  that had allele frequency data in dbsnp.
the distribution of sample sizes for snps was shown in additional file  <dig> figure s <dig>  we listed the top  <dig> predictions for both groups with respect to their posterior probabilities of imprinting in table  <dig> .
surprisingly, we found only  <dig> common genes between the 57-gene list and the 94-gene list.
functional enrichment analysis performed in ingenuity® pathway analysis  showed that the two lists of genes were significantly enriched  in certain bio functions but not in any canonical pathways defined by ipa .
we also compared our predictions with the  <dig> genes that are subject to random monoallelic expression  <cit>  and found none of our predicted genes in either group overlapped with them, which further validates the quality of our predictions and strongly suggests that our predictions are not affected by random monoallelic expression.
abbreviations: rs#-snp identification number, chr-chromosome, str-strand, ss-sample size.
model analysis based on independent test sets
we collected  <dig> validated imprinted genes and  <dig> validated non-imprinted genes  and used them to check for false-positive and false-negative predictions made by dspig.
based on the mrna-seq data of group i samples,  <dig> of the  <dig> imprinted genes and  <dig> of the  <dig> non-imprinted genes had allelic counts for known snps; based on the mrna-seq data of group ii samples,  <dig> of the  <dig> imprinted genes and  <dig> of the  <dig> non-imprinted genes had allelic counts for known snps.
under the same criteria used to predict imprinted genes in table  <dig>  dspig identified  <dig> out of  <dig> imprinted genes from group i samples and  <dig> out of  <dig> imprinted genes from group ii samples, based on the mrna-seq data; under the same criteria, dspig misidentified  <dig> out of  <dig> non-imprinted genes from group i samples and  <dig> out of  <dig> non-imprinted genes from group ii samples.
moreover, among those imprinted genes, only  <dig> out of  <dig>  and  <dig> out of  <dig>  have sample sizes > <dig>  on the contrary, all  <dig> dspig-identified imprinted genes  have sample size > <dig>  therefore, for imprinted genes with sample size >  <dig>   <dig> out of  <dig> genes and  <dig> out of  <dig> genes could be identified by dspig for group i and group ii, respectively.
again, this showed that the sensitivity of dspig increased as sample size increased and dspig could probably identify more imprinted genes if the number of sequenced tissue samples further increased.
this also agreed with the simulation results.
interestingly, some of the validated imprinted genes had very small posteriors , which indicated that they had biallelic expression  in certain tissues .
posteriors were calculated based on group i or ii .
the “tissue” column indicates where genes were most likely biallelically expressed.
genes with contradictory evidence  were discarded.
igf <dig>  osbpl <dig>  smpd <dig> and gpr172a had at least  <dig> snps with posteriors less than 10- <dig>  abbreviations are the same as in table  <dig>
we used the same sets of imprinted and non-imprinted genes to determine the cut-off of the posteriors used in the prediction of imprinted genes.
because most genes  are expected to be non-imprinted, the cut-off has to yield a very high specificity  so that the overall fdr of our predictions can be low enough  for further validations.
after trying different cut-offs , we found  <dig>  to be the most appropriate cut-off in terms of the validated gene sets because  increasing the cut-off from  <dig>  only lowered the sensitivity while left the specificity unchanged , and  decreasing the cut-off from  <dig>  lowered the specificity while the sensitivity didn’t change a lot, which substantially increased the fdr of predictions .
we also showed the roc curves for dspig in additional file 7: figure s <dig>
candidates for experimental validation
we chose top  <dig> candidate genes from our predictions and listed them in table  <dig>  except one snp , all the snps in table  <dig> have high minor allele frequencies , which indicates >30% chance of observing heterozygous alleles in experiments.
in addition, these genes also met at least one of the following three criteria:  their snps  had a relatively low fdr ,  they had multiple snps with posteriors >  <dig>  , and  they were located near existing imprinted genes   <cit> .
these additional criteria further increased the possibility of imprinting.
boldface indicates proximity to known imprinted genes.
underline indicates a gene has multiple  snps with posteriors >  <dig> .
genes are sorted by their posteriors.
rpl10ap <dig> and fer are listed for both group i and group ii.
abbreviations are the same as in table  <dig>
detection of allele-preferred expression
by investigating the biallelically expressed genes identified in the mrna-seq data, we found that at least  <dig> genes expressed significantly more transcripts from one specific allele than the other among various individuals and tissues .
this indicated that the difference between expression levels of the two alleles was not caused by sequencing errors or stochastic effects in rt-pcr.
in future, as more mrna-seq data are generated, if more genes with one specific allele always under-expressed are observed, we would speculate that a sophisticated mechanism  may exist to explain this type of allelic preference in gene expression for biallelically expressed genes.
the “ratio” column shows the proportion of samples with one specific allele’s transcript level higher than the other’s.
p-values have been adjusted by “bh” correction  <cit> .
abbreviations are the same as in table  <dig>
web-based service and r package for dspig
we have provided a web-based service for dspig at http://www.shoudanliang.com/dspig/.
users need to upload either mapped mrna-seq data in the supported format or processed data files containing allelic counts for each snp .
after uploading the data, users may set the values for  the cut-off for posteriors,  the average sequencing error rate,  the prior of imprinting, and  qs .
in addition, users need to specify the human genome build and the snp build for dspig.
our server will calculate the posterior probabilities of imprinting for each gene and email the results back to the users.
in the final submission form, users may request additional analysis, such as suggesting tissues where known imprinted genes most likely have biallelic expression.
in addition, we have made an r package for dspig, which is available both on the website and in the supplementary materials .
the instruction and the sample files for this r package are in additional file  <dig>  the web service and the r package generate the same result on the predictions of imprinted genes.
for reference, we have also provided the annotated code for dspig  used in this study as additional file  <dig> and on the website.
discussion
dspig operates under the assumption that if a gene is biallelically expressed, the transcript levels of two alleles are the same.
however, because imbalance of allelic expression has been widely detected in human tissues  <cit> , our assumption may not always be correct.
indeed, within our 29-tissue samples, allelic imbalance was observed even after excluding possible expression from imprinted genes .
however, our simulation studies showed that our model can still distinguish imprinted genes from non-imprinted ones in most situations , unless one specific allele is always expressed at <13% of the other allele’s expression level across different samples .
stochastic rt-pcr amplification  <cit>  is not of particular concern in dspig because this has been taken into account as allelic imbalances.
dspig is sensitive to biallelic expression and unlikely to falsely predict imprinted genes.
a gene will get a very low posterior probability of imprinting when it obviously has biallelic expression in one tissue, even if only one allele of this gene is observed in the transcriptomes from all other tissues .
for example, a gene may have monoallelic expression because of allele-specific differences caused by either heterozygous snps or somatic mutations in the promoter region of the gene; as long as monoallelic expression caused by these conditions is not present population-wide and a large number  of independent tissue samples is used in dspig, these genes will not be falsely predicted as imprinted genes.
however, for the same reason, if transcriptomes are collected from various tissues, it becomes very hard for dspig to detect tissue-specific imprinted genes.
this partially explains why dspig predicted much less imprinted genes in group i samples than in group ii samples .
one advantage of dspig is that it is able to predict imprinted genes without sequencing the genotype.
although a homozygous allele will also lead to identification of only one allele in the transcriptome, the result will not elevate the posterior belief of imprinting in our model.
this is very important and practical for human because we no longer need to sequence the genome for genotypes.
however, a disadvantage is that dspig cannot tell the parent of origin for predicted imprinted genes, which can be verified only by other studies.
one obvious limit of dspig is that it was modeled based on single snps.
this means, if a gene has more than  <dig> snp site in its exons, it may have different single-snp-based allelic counts and thus different posterior probabilities from dspig.
therefore, one single strong posterior cannot determine the imprinting status of this gene.
instead, we have to look into all snp sites of each gene and make sure no contradictory posterior exists for our predictions .
in future, a possible improvement would be integration of all snp information of a gene and calculate a single posterior to predict imprinting status.
CONCLUSIONS
in this paper, we proposed a new method – dspig, applicable to all mammals with genomic imprinting, to predict imprinted genes based on mrna-seq data of various independent tissues.
with enough sequenced samples, dspig is capable of predicting imprinted genes on a genome-wide scale with expected low fdrs.
the power of dspig will be further enhanced after more data generated by mrna-seq technology become available in the near future.

 diploid eukaryotic species inherit two copies of the same gene from both parents. if one allele fails to work properly, the other allele can still implement a gene’s cellular function. for some genes, however, this protective mechanism is disabled because only one allele is expressed and its failure probably leads to cellular malfunction. these monoallelically expressed genes can be classified into one of three categories: x-inactivated genes, which are regulated by a random process in which one of the two x chromosomes present in female mammals is silenced; autosomal genes subject to random monoallelic gene expression, such as the t cell receptors and natural killer cell receptors; and autosomal imprinted genes, which express from only one of the two alleles according to their parental origin. imprinted genes play important functional roles in the control of embryonic growth and development, as well as in post-natal development. as imprinted genes are expressed from only one of the two parental chromosomes, a de facto haploid state is caused by imprinting and leads to asymmetric functions of parental genomes and loss of diploid protection against recessive mutations. thus, imprinting dysregulation is linked to numerous human genetic diseases, such as developmental disorders and cancers. in addition, environmental factors can influence gene expression by targeting imprinted genes. because imprinted genes are more susceptible to disease than non-imprinted genes, it is of great importance to identify novel imprinted genes for human. identifying imprinted genes experimentally has traditionally been a slow process, and the number of validated ones is much lower than the previous estimation. however, recently developed high-throughput screening approaches and recently identified dna sequence characteristics have led to the proposal of several new methods to predict imprinted genes on a global scale. with advances in next-generation sequencing technology, mrna-seq is becoming a powerful tool for transcriptome profiling. it can generate not only the number of reads mapped to exons, which reflects the expression levels of a gene, but also the actual sequence, which may identify the allele from which the mrna is expressed. therefore, inference can be made for predicting imprinted genes. for example, by sequencing whole transcriptomes from mice embryos, babak et al. measured allelic expression bias and identified six novel imprinted genes. however, to our knowledge, prediction of imprinted genes by deeply sequencing transcriptomes from multiple independent tissues is still an open problem. in this study, we proposed a bayesian model – dspig – to predict imprinted genes based on allelic expression inferred from observed snps in mrna-seq data of independent human tissues. with dspig, we were able to measure the imbalance of allelic expression among various tissues and calculate the posterior probability of imprinting status for each gene. under a stringent probabilistic cut-off of the posteriors and other reasonable biological criteria, we identified potentially imprinted genes from diverse human tissues and potentially imprinted genes from cerebellar cortices, with an expected low false discovery rate. furthermore, analysis of allelic expression of the same genes among different tissues revealed that, in some cases, even if a gene was biallelically expressed, one allele always had higher expression level than the other.  statistical model development monoallelic expression generally falls into one of three categories: imprinted expression, random monoallelic expression and x-inactivation, all of which express only one of two alleles in a single cell. at a tissue level, however, random monoallelic expression will allow both alleles to be detected in total rna because of the “mosaicism” of the tissue. because our study was based on whole transcriptomes of tissue samples, random monoallelic expression was reasonably considered as biallelic expression when averaged over the entire tissue. x-inactivation was also excluded from this study by discarding all predictions on the x chromosome. thus imprinting is the most likely cause of the observed monoallelic expression among transcriptomes of different tissues even though we cannot infer the parent of origin. we used known snps from dbsnp to distinguish and count the two alleles of each gene. if a gene was imprinted, we expected to observe only one of the two alleles of each snp in the exon region from the whole transcriptome. with the allelic counts obtained from the mrna-seq data, we developed a bayesian model to compute the posterior probability of imprinting based on each single snp. suppose we have sequenced transcriptomes from n independent tissue samples. for each sample, we count the alleles of all known snps, discarding those with counts. for each snp, let the allelic counts be:, …,, where xk and yk are the counts for two alleles x and y in the sample k. because each gene can only have two statuses: imprinted or non-imprinted, we consider as a binary event vector for the imprinting status. by denoting data = {, …, }, we have by bayes’ theorem: pri|data=prdata|i×priprdata|i×pri+prdata|ni×prni where pr is the posterior probability of imprinting and pr is the prior of imprinting. based on current knowledge of prevalence of imprinted genes, we set the prior pr=1%, thus pr = – pr = 99%. since samples were independent of each other and genotype has only possible combinations, we denote the genotype as follows: δ={1the genotype isxx2the genotype isxy3the genotype isyy assuming p and q are the allele frequencies for allele x and y, p + q = according to hardy-weinberg equilibrium, the prior probabilities for the three genotypes are calculated as follows: prδ={=prx2=p2δ=1=prpr=2pqδ=2=pry2=q2δ= since values of p and q can be retrieved from dbsnp, p and q are treated as constants. we used the law of total probability to calculate the likelihood pr as follows: prdata|i=∏k=1nprxk,yk|in≥ =∏k=1nprxk,yk|i,δ=1×p2+prxk,yk|i,δ=2×2pq+prxk,yk|i,δ=3×q by assuming that the transcript levels of a gene’s two alleles are equal if the gene is biallelically expressed, two different alleles, if both expressed, have the same chance to be sequenced by mrna-seq, and unexpressed alleles with counts > are due to sequencing errors, we had the following derivation: prxk,yk|i,δ=1=fyk;nk,pe prxk,yk|i,δ=3=fxk;nk,pe prxk,yk|i,δ= =prxk,yk|i,δ= θ=0×prθ=0+prxk,yk|i,δ= θ=1×prθ= =fyk;nk,pe×12+fxk;nk,pe× here, f denotes binomial distribution, nk = xk + yk is assumed fixed, and pe is the averaged sequencing error rate for each snp. the binary variable θ is defined as follows θ={0onlyxcan be expressed due to imprinting1onlyycan be expressed due to imprinting. since x and y have an equal chance to be inherited from either maternal or paternal genome, x and y have an equal chance to be expressed in imprinted genes. hence, prθ={1/2θ=01/2θ= similarly for the likelihood function pr, we have: prdata|ni=∏k=1nprxk,yk|nin≥ =∏k=1nprxk,yk|ni,δ=1×p2+prxk,yk|ni,δ=2×2pq+prxk,yk|ni,δ=3×q prxk,yk|ni,δ=1=fyk;nk,pe prxk,yk|ni,δ=3=fxk;nk,pe based on our three assumptions –, pr follows a binomial distribution with p=. therefore, we have: prxk,yk|ni,δ=2=xk+ykxk12xk+yk computation is performed separately for each single snp. therefore, a posterior probability of imprinting for a gene is associated with a specific snp in this gene, unless otherwise specified. simulation-based model analysis we generated allelic counts from simulated data by taking into account imprinting status, snp frequency, and sequencing error. we then applied dspig to estimate the sensitivity, specificity and the fdr. we generated two sets of allelic counts under the assumption that the locus was either imprinted or not. the number of reads for each allele was generated assuming the presence of one dominating allele plus sequencing error for imprinted case; and presence of equal amount of two alleles plus sequencing error for non-imprinted case. the generated allelic counts followed a distribution similar to the actual mrna-seq data. more details of the procedure were illustrated in additional file figure s given an allele frequency, a sequencing error and a prior of imprinting, the posterior probability calculated from the allelic counts generated for imprinted genes approached as the sample size increased, while concomitantly the posterior probability for biallelically expressed approached. with minor allele frequencies between to, sequencing errors between to and priors between to, we obtained similar results as figure 1a,b. using as the cut-off for posteriors pr, sensitivity of our model-based prediction exceeded % when sample size was > and specificity exceeded % when sample size was >. under the same cut-off and the allele frequency, the fdr of predicted imprinted genes decreased as sample size increased: it dropped to 5% and 1% as sample size exceeded and respectively. using the same cut-off, we also examined how fdr varied when minor allele frequency changed from to and sample size increased from to. based on these simulations, we were able to provide estimated fdrs for most of our predictions in this study. predictions of imprinted genes we collected two previously published mrna-seq data sets. one set included diverse tissue samples from wang et al. and the other set had cerebellum cortex samples from mudge et al.. wang et al. showed that, in terms of alternative isoform expression, cerebellum tissues were clustered together and the diverse samples were more closely correlated. here we performed hierarchical clustering based on the imprinting-inclined snps and obtained similar results. as we expected, samples from the cerebellar cortex were clustered together, with caucasian and african american separated in two sub-clusters. using caucasian allele frequency on african american samples, however, yielded a sub-cluster without separation between the two ethnicities. this suggests that the separation observed in figure 3a was due to differences in minor allele frequencies. as a test set, breast cancer cell line samples were clustered together in both cases. compared with other non-cerebellum samples, the brain sample had higher correlation with cerebellum samples in both cases, which is sensible biologically. the result that the diverse tissue samples were clustered together could be caused by many factors such as different experimental conditions between group i and group ii samples. using dspig, we predicted imprinted genes for group i and group ii separately. to call a gene imprinted, we checked the posteriors of all the snps on the same gene to make sure that the gene had at least one snp with a posterior >, which was the same cut-off used in simulations, and all the other snps did not show any contradictory evidence. after applying dspig to the mrna-seq data and using the above criteria, we predicted potentially imprinted genes for group i samples and potentially imprinted genes for group ii samples out of a total of genes that had allele frequency data in dbsnp. the distribution of sample sizes for snps was shown in additional file figure s we listed the top predictions for both groups with respect to their posterior probabilities of imprinting in table. surprisingly, we found only common genes between the 57-gene list and the 94-gene list. functional enrichment analysis performed in ingenuity® pathway analysis showed that the two lists of genes were significantly enriched in certain bio functions but not in any canonical pathways defined by ipa. we also compared our predictions with the genes that are subject to random monoallelic expression and found none of our predicted genes in either group overlapped with them, which further validates the quality of our predictions and strongly suggests that our predictions are not affected by random monoallelic expression. abbreviations: rs#-snp identification number, chr-chromosome, str-strand, ss-sample size. model analysis based on independent test sets we collected validated imprinted genes and validated non-imprinted genes and used them to check for false-positive and false-negative predictions made by dspig. based on the mrna-seq data of group i samples, of the imprinted genes and of the non-imprinted genes had allelic counts for known snps; based on the mrna-seq data of group ii samples, of the imprinted genes and of the non-imprinted genes had allelic counts for known snps. under the same criteria used to predict imprinted genes in table dspig identified out of imprinted genes from group i samples and out of imprinted genes from group ii samples, based on the mrna-seq data; under the same criteria, dspig misidentified out of non-imprinted genes from group i samples and out of non-imprinted genes from group ii samples. moreover, among those imprinted genes, only out of and out of have sample sizes > on the contrary, all dspig-identified imprinted genes have sample size > therefore, for imprinted genes with sample size > out of genes and out of genes could be identified by dspig for group i and group ii, respectively. again, this showed that the sensitivity of dspig increased as sample size increased and dspig could probably identify more imprinted genes if the number of sequenced tissue samples further increased. this also agreed with the simulation results. interestingly, some of the validated imprinted genes had very small posteriors, which indicated that they had biallelic expression in certain tissues. posteriors were calculated based on group i or ii. the “tissue” column indicates where genes were most likely biallelically expressed. genes with contradictory evidence were discarded. igf osbpl smpd and gpr172a had at least snps with posteriors less than 10- abbreviations are the same as in table we used the same sets of imprinted and non-imprinted genes to determine the cut-off of the posteriors used in the prediction of imprinted genes. because most genes are expected to be non-imprinted, the cut-off has to yield a very high specificity so that the overall fdr of our predictions can be low enough for further validations. after trying different cut-offs, we found to be the most appropriate cut-off in terms of the validated gene sets because increasing the cut-off from only lowered the sensitivity while left the specificity unchanged, and decreasing the cut-off from lowered the specificity while the sensitivity didn’t change a lot, which substantially increased the fdr of predictions. we also showed the roc curves for dspig in additional file 7: figure s candidates for experimental validation we chose top candidate genes from our predictions and listed them in table except one snp, all the snps in table have high minor allele frequencies, which indicates >30% chance of observing heterozygous alleles in experiments. in addition, these genes also met at least one of the following three criteria: their snps had a relatively low fdr, they had multiple snps with posteriors >, and they were located near existing imprinted genes. these additional criteria further increased the possibility of imprinting. boldface indicates proximity to known imprinted genes. underline indicates a gene has multiple snps with posteriors >. genes are sorted by their posteriors. rpl10ap and fer are listed for both group i and group ii. abbreviations are the same as in table detection of allele-preferred expression by investigating the biallelically expressed genes identified in the mrna-seq data, we found that at least genes expressed significantly more transcripts from one specific allele than the other among various individuals and tissues. this indicated that the difference between expression levels of the two alleles was not caused by sequencing errors or stochastic effects in rt-pcr. in future, as more mrna-seq data are generated, if more genes with one specific allele always under-expressed are observed, we would speculate that a sophisticated mechanism may exist to explain this type of allelic preference in gene expression for biallelically expressed genes. the “ratio” column shows the proportion of samples with one specific allele’s transcript level higher than the other’s. p-values have been adjusted by “bh” correction. abbreviations are the same as in table web-based service and r package for dspig we have provided a web-based service for dspig at http://www.shoudanliang.com/dspig/. users need to upload either mapped mrna-seq data in the supported format or processed data files containing allelic counts for each snp. after uploading the data, users may set the values for the cut-off for posteriors, the average sequencing error rate, the prior of imprinting, and qs. in addition, users need to specify the human genome build and the snp build for dspig. our server will calculate the posterior probabilities of imprinting for each gene and email the results back to the users. in the final submission form, users may request additional analysis, such as suggesting tissues where known imprinted genes most likely have biallelic expression. in addition, we have made an r package for dspig, which is available both on the website and in the supplementary materials. the instruction and the sample files for this r package are in additional file the web service and the r package generate the same result on the predictions of imprinted genes. for reference, we have also provided the annotated code for dspig used in this study as additional file and on the website. discussion dspig operates under the assumption that if a gene is biallelically expressed, the transcript levels of two alleles are the same. however, because imbalance of allelic expression has been widely detected in human tissues, our assumption may not always be correct. indeed, within our 29-tissue samples, allelic imbalance was observed even after excluding possible expression from imprinted genes. however, our simulation studies showed that our model can still distinguish imprinted genes from non-imprinted ones in most situations, unless one specific allele is always expressed at <13% of the other allele’s expression level across different samples. stochastic rt-pcr amplification is not of particular concern in dspig because this has been taken into account as allelic imbalances. dspig is sensitive to biallelic expression and unlikely to falsely predict imprinted genes. a gene will get a very low posterior probability of imprinting when it obviously has biallelic expression in one tissue, even if only one allele of this gene is observed in the transcriptomes from all other tissues. for example, a gene may have monoallelic expression because of allele-specific differences caused by either heterozygous snps or somatic mutations in the promoter region of the gene; as long as monoallelic expression caused by these conditions is not present population-wide and a large number of independent tissue samples is used in dspig, these genes will not be falsely predicted as imprinted genes. however, for the same reason, if transcriptomes are collected from various tissues, it becomes very hard for dspig to detect tissue-specific imprinted genes. this partially explains why dspig predicted much less imprinted genes in group i samples than in group ii samples. one advantage of dspig is that it is able to predict imprinted genes without sequencing the genotype. although a homozygous allele will also lead to identification of only one allele in the transcriptome, the result will not elevate the posterior belief of imprinting in our model. this is very important and practical for human because we no longer need to sequence the genome for genotypes. however, a disadvantage is that dspig cannot tell the parent of origin for predicted imprinted genes, which can be verified only by other studies. one obvious limit of dspig is that it was modeled based on single snps. this means, if a gene has more than snp site in its exons, it may have different single-snp-based allelic counts and thus different posterior probabilities from dspig. therefore, one single strong posterior cannot determine the imprinting status of this gene. instead, we have to look into all snp sites of each gene and make sure no contradictory posterior exists for our predictions. in future, a possible improvement would be integration of all snp information of a gene and calculate a single posterior to predict imprinting status.  in this paper, we proposed a new method – dspig, applicable to all mammals with genomic imprinting, to predict imprinted genes based on mrna-seq data of various independent tissues. with enough sequenced samples, dspig is capable of predicting imprinted genes on a genome-wide scale with expected low fdrs. the power of dspig will be further enhanced after more data generated by mrna-seq technology become available in the near future. 
here we propose a bayesian model to predict imprinted genes on the basis of allelic expression observed in mrna-seq data of independent human tissues. prediction of imprinted genestranscriptome deep sequencingmrna-seqbayesian modelanalysis of allelic expression our model was capable of identifying imprinted genes with high sensitivity and specificity and a low false discovery rate when the number of sequenced tissue samples was fairly large, according to simulations. by applying dspig to the mrna-seq data, we predicted imprinted genes in cerebellum samples and imprinted genes in diverse tissue samples with expected low false discovery rates. to date, however, fewer than imprinted genes have been identified in the human genome. for ease of use, we have set up a web service and also provided an r package for dspig at http://www.shoudanliang.com/dspig/. with the prevalence of the mrna-seq technology, dspig has become a useful tool for analysis of allelic expression and large-scale prediction of imprinted genes. with simulations, we further analyzed how imbalanced allelic expression of non-imprinted genes or different minor allele frequencies affected the predictions of dspig.
 dysregulation of imprinted genes, which are expressed in a parent-of-origin-specific manner, plays an important role in various human diseases, such as cancer and behavioral disorder. to date, however, fewer than imprinted genes have been identified in the human genome. the recent availability of high-throughput technology makes it possible to have large-scale prediction of imprinted genes. here we propose a bayesian model to predict imprinted genes on the basis of allelic expression observed in mrna-seq data of independent human tissues.  our model was capable of identifying imprinted genes with high sensitivity and specificity and a low false discovery rate when the number of sequenced tissue samples was fairly large, according to simulations. by applying dspig to the mrna-seq data, we predicted imprinted genes in cerebellum samples and imprinted genes in diverse tissue samples with expected low false discovery rates. we also assessed dspig using previously validated imprinted and non-imprinted genes. with simulations, we further analyzed how imbalanced allelic expression of non-imprinted genes or different minor allele frequencies affected the predictions of dspig. interestingly, we found that, among biallelically expressed genes, at least genes expressed significantly more transcripts from one allele than the other among different individuals and tissues.  with the prevalence of the mrna-seq technology, dspig has become a useful tool for analysis of allelic expression and large-scale prediction of imprinted genes. for ease of use, we have set up a web service and also provided an r package for dspig at http://www.shoudanliang.com/dspig/. prediction of imprinted genestranscriptome deep sequencingmrna-seqbayesian modelanalysis of allelic expression 