Joint estimation of genetic and parent-of-origin effects using RNA-seq data from human
RNA sequencing allows one to study allelic imbalance of gene expression, which may be due to genetic factors or genomic imprinting. It is desirable to model both genetic and parent-of-origin effects simultaneously to avoid confounding and to improve the power to detect either effect. In a study of experimental cross, separation of genetic and parent-of-origin effects can be achieved by studying reciprocal cross of two inbred strains. In contrast, this task is much more challenging for an outbred population such as human population. To address this challenge, we propose a new framework to combine experimental strategies and novel statistical methods. Specifically, we propose to collect genotype data from family trios as well as RNA-seq data from the children of family trios. We have developed a new statistical method to estimate both genetic and parent-of-origin effects from such data sets. We demonstrated this approach by studying 30 trios of HapMap samples. Our results support some of previous finding of imprinted genes and also recover new candidate imprinted genes.
Joint estimation of genetic and parent-of-origin effects using RNA-seq data from human–References 2016 \artmonthOctober
NA-seq; imprinting; allele-specific expression; expression Quantitative Trait Loci (eQTL)
By combining RNA-seq data and phased genotype of a diploid genome, one can estimate allele-specific expression (ASE). Allelic imbalance of gene expression (i.e., the ASE of the two alleles are different) may be controlled by genetic factors or genomic imprinting. Genetic loci that affect gene expression variation are often referred to as gene expression quantitative trait loci (eQTLs). Sophisticated statistical methods have been developed to identify eQTLs using both ASE and total expression (Sun, 2012; Sun and Hu, 2013; Hu
et al., 2015). Genomic imprinting is a phenomenon in which the expression of a gene depends whether an allele has maternal or paternal origin. Some genes are (completely) imprinted in the sense that they are exclusively expressed from one of the two parental alleles. However, recent studies have shown that many genes are partially imprinted, i.e., both alleles are expressed while the allele from one parent has higher expression the other allele (Crowley
et al., 2015). Although only a small percentage of human genes are (completely or partially) imprinted, studying of imprinting is of great importance for at least two reasons. First, the imprinted genes play important roles in many common diseases such as obesity or psychiatric disorders (Peters, 2014). Second, genomic imprinting is also a great model to study epigenetic regulation of gene expression (Barlow, 2011).
It is important to jointly estimate genetic and parent-of-origin effects on allelic imbalance of gene expression in order to avoid confounding of these two factors and to improve the power to detect either effect. This goal can be achieved by studying RNA-seq data from reciprocal cross of two inbred strains, e.g., F1 mouse population (Zou et al., 2014; Crowley
et al., 2015). Specifically, denote two inbred strains by and , a reciprocal cross refers to two crosses and . Following standard naming convention, a cross is written as maternal strain paternal strain, and thus in a cross , is the maternal strain and is the paternal strain. To make it more clear, we write the two crosses as and , where the superscript indicates maternal or paternal origin. In one cross , if the expression of allele is larger than allele, it is not clear whether it is due to genetic effect (i.e., allele to higher expression) or parent-of-origin effect (i.e., the maternal allele has higher expression). However, given a reciprocal cross, such confounding can be resolved.
In an outbred population such as the human population, separation of genetic and parent-of-origin effect becomes much more challenging. Previous studies have employed genotyping arrays to assess allelic imbalance using both RNA and genomic DNA (gDNA) samples and identify potential imprinted genes if the degrees of allelic imbalance is significantly different between RNA and gDNA samples. Those potential imprinted genes are then screened using family pedigree (Morcos et al., 2011) or experimental validation (Barbaux et al., 2012). A more direct approach is to collect RNA-seq data from the children of family trios (a family trio includes both parents and one child) and collect genotype data from all three family members of each trio. The genotype can be phased and imputed to a larger number of SNPs. Then ASE can be quantified based on phased genotype data. In this work, we developed a new statistical method to separate genetic and parent-of-origin effect while jointly modeling ASE and total read count (TReC) across a group of family trios.
We assume genetic effects can be captured by a small number of cis-acting eQTLs. Here we define cis-eQTLs as those eQTLs that influence allelic-imbalance of gene expression (Sun and Hu, 2013). This assumption can be justified by the following rationales. First, by definition, only cis-acting eQTLs are relevant for ASE study. Second, the number of cis-acting eQTLs is limited because most cis-acting eQTLs are local eQTLs that are close to the gene of interest, and due to linkage dis-equilibrium (LD) structure around a gene, the number of independent cis-eQTLs of a gene is relatively small. In the following, we assume there is only one cis-eQTL to simplify the discussion. Our method can be easily extended to the cases with multiple cis-eQTLs.
2.1 Allele specific read counts (ASReCs)
To simplify the notation, we assume the haplotypes connecting a candidate eQTL and the gene of interest are known. Let and be the two haplotypes of a gene of interest in the -th individual (for ), and let the and be the allele-specific read counts of these haplotypes and , respectively. Let . Denote the two alleles of the candidate eQTL as and , and denote its genotype in the -th individual as . We assign different meanings to genotypes and such that means haplotypes and harbor the and alleles, respectively, and means haplotypes and harbor the and , respectively. We model by a beta-binomial distribution (denoted by ):
where is over-dispersion parameter, and and are defined as follows:
2.2 Total read counts (TReCs)
The TReC of the gene of interest in the -th individual, denoted by , is modeled by a negative binomial distribution with mean and over-dispersion parameter , denoted by . We can write the mean structure for negative binomial distribution as:
where , , is the regression coefficient for the -th covariate (e.g., age, gender, batch effects). is defined as:
Note that the additive genetic and parent-of-origin effects are parametrized by and , which are the same as the and in equation (1) for allele-specific read counts.
2.3 Joint likelihood
The joint likelihood of total read count (TReC) and allele-specific read count (ASReC) is
where . In this model we assume common genetic and parent-of-origin effect for TReC and ASReC. We test the genetic and parent-of-origin effects by testing whether or equals to 0, respectively, with the likelihood ratio test.
2.4 Optimization Algorithm
For a given initial values for non-linear terms (, , , ), our optimization algorithm includes the following steps.
Step 1: Optimize linear terms given the initial values of non-linear terms:
where is a diagonal matrix.
Step 2: Iteratively estimate and . Note that the following two steps are redundant, but improves the robustness of the algorithm.
Optimize and together using BFGS method.
Optimize and separately using Brent algorithm.
Step 3: optimize over-dispersion parameter and separately using Brent algorithm. For stability we limit range of over-dispersion to be between and .
Step 4: If likelihood change is larger than a small number , go to step 1, otherwise finish the estimation process.
We simulated ASE and TReC from the model described by joint likelihood in equation (3). We chose the smallest sample size to be 32, which is similar to the sample size we have for real data analysis. To demonstrate the asymptotic properties, we also simulated data of larger sample sizes of 64, 128, and 256. We set the over-dispersion parameters for beta-binomial () and negative binomial () distributions to be 1/4, and 4/3, respectively, which are fairly typical in real data. We further scale the expected counts so that mean total read count was about 250. We set the proportion of counts that are allele-specific to be 10% which is reasonably close to observed value in our real data.
Our simulation results show that the estimates of and from our method are unbiased. In contrast, if a model is fit with only genetic or parent-of-origin effect while in fact both effects are present, the parameter estimates have significant bias (Figure 1).
We also compared the type I error and the power of our joint model versus naive fit using function to fit Negative-Binomial model using total read counts only and function to fit Beta-Binomial model using allele specific counts only. As shown in Table 1, the simple models don’t control type I error as well as the joint model and that the joint model has higher power to detect both genetic and parent-of-origin effects.
|Parameters||Genetic||Parent of Origin|
The statistical power of the joint model is illustrates in Figure 2. Even at sample size of 32, our method has around 80% of power to detect either genetic or parent-of-origin effect at two-fold change, which corresponds to effect size . The power to detect parent-of-origin effect is higher than the power to detect genetic effect. This is because the ASE of all samples can be used to quantify parent-of-origin effect (i.e., comparing paternal vs. maternal allele). In contrast, ASE can be used to quantify genetic effect only if the genotype of the candidate eQTL is heterozygous, which only happen in a subset of the samples.
The computational time of our method scales well with the sample size. The time per gene remains under a second for sample sizes under 200 individuals: see Figure 6 in Supplementary Materials .
4.1 Data collection
We collected RNA-seq data from 30 HapMap CEU (Caucasian) samples (15 males + 15 females). All of these samples are children of family trios, where the genotypes/haplotypes of three family members have been reported in the HapMap project, and their lymphoblastoid cell lines are available through Coriell (http://www.coriell.org/). For most of these samples, the RNA reads were 150 bp paired-end reads, with an additional smaller run of 75 bp paired-end reads. The median of the total number of reads for these 30 samples is about 20 million. These reads were mapped to hg38 human reference genome using Tophat2. Hapmap project genotyped about 3.9 million SNPs for these 30 trios. We phased and imputed the genotypes of these 30 trios using shapeit2 (Delaneau et al., 2014) and impute2 (Howie et al., 2012), against 1000 Genome reference panel containing 2,504 individuals with 82 million SNPs. Finally, based on phased and imputed genotype, we extracted allele-specific reads (i.e., those RNA-seq reads that overlap with at least one heterozygous SNP), and counted the number of allele-specific reads for each haplotype of a gene.
4.2 Identification of candidate cis-eQTLs
The parents of these 30 HapMap family trios are part of the samples included in the 1000 Genomes project (1000 Genomes Project Consortium, 2012). To improve the power and precision for eQTL mapping, we first mapped cis-acting eQTLs using the Caucasian samples of the 1000 Genome Project. Then given the eQTLs, we jointly estimated genetic and parent-of-origin effects in the 30 children of HapMap trios. We downloaded fastq files for 462 samples from the Geuvadis consortium (Lappalainen et al., 2013). We mapped all the reads to hg38 reference, and then performed the same process of phasing and imputation of the genotypes of these samples. Then we calculated total read counts and allele-specific read counts per gene and per sample. After filtering out non-European samples as well as samples with low quality in terms of systematic bias of allele-specific expression (see section A.9 in Supplementary Materials for details), we performed eQTL mapping for 227 samples using the TReCASE method (Hu et al., 2015). For each gene, we identified one candidate cis-acting eQTL by selecting the one with the smallest p-value. Note that some genes do not have any cis-acting eQTL because all the nearby SNPs show certain evidence of trans-acting eQTL, i.e., the eQTL genotype is associated with total expression, but not allele-specific expression. The eQTL effect sizes estimated from these 227 samples are consistent with the estimates from the 30 children of HapMap trios (see section B.1 in Supplementary Materials for details).
4.3 Identification of imprinted genes
We used our method to estimate genetic (eQTL) and imprinting effects for 12,386 genes in the 30 children of the family trios. These genes were selected because they had enough expression in the 30 samples and they had at least one candidate cis-eQTL based on eQTL mapping results from the 227 samples of 1000 Genome Project. For the negative binomial model of total read count (TReC), we fit our model with additional covariates to capture batch effect (the RNA-seq data were collected through 3 batches with 10 samples per batch) and read-depth. No additional covariate is needed for the analysis of allele-specific read counts (ASReC) because our model for ASE compared the expression of one allele versus another allele, and thus the effects of such covariates are canceled. We found 16 genes with significant imprinting effects at q-value cutoff 0.05 (Table 2).
Out of these 16 genes, 10 were known to be imprinted from previous studies and 6 were novel candidates for imprinting. For 14 of these 16 genes, the paternal allele had higher expression than the maternal allele. For all of those 10 known imprinted genes, our estimates of allelic imbalance agree with what were reported before. At a more liberal cutoff of q-value 0.25, we identified 15 additional candidates of imprinted genes. By manually examining the expression pattern of these 15 genes, we suggest that 9 of them missed the q-value cutoff 0.05 due to power - either because parent-of-origin effect was smaller (3 genes) or the number of allele-specific reads was small (6 genes). The other 6 are likely false discoveries, because there is no apparent pattern of allelic imbalance for one parental allele.
We illustrate the read count data of a clearly imprinted gene, ZNF497, which has higher expression on maternal allele (Figure 3).
The imprinting effect can be observed from both TReC and ASReC. We denote the genotype such as the first allele is maternal allele, i.e., genotype AB means A and B are from maternal and paternal allele, respectively. For TReC, the two groups with genotype AA and AB have similar expression because they share the same maternal allele and maternal allele has higher expression. Similarly, the two groups BA and BB have similar expression. The effect of imprinting is more apparent in allele-specific reads where we observe the proportion of reads from paternal allele is far below 50%.
We selected known imprinted genes based on the list reported by Morison et al. (2005) or those genes recorded in the Geneimprint database (Jirtle, 2016). There is a total of 90 genes classified as imprinted. Out of those, 32 were expressed in our dataset and had a candidate eQTL. Of those 32 genes 10 were found to be significant (q-value 0.05) by our method. For several genes (such as RB1, KCNQ1, PEG3, and PLAGL1), we observed signal of imprinting, but the signal was too weak to produce significant q-value. In general, though, we observed that even for insignificant results, those with relatively smaller q-value tend to have estimated imprinting direction matched with the reported imprinting direction.
4.4 Locations of genes with parent-of-origin effect
The imprinted genes in mouse genome tend to form clusters. It is interesting to check whether this is true for human genome. Based on the findings from our study, the imprinted genes on chromosomes 15 and 16 form two clusters (Figure 4 and Table 3). However, the total number of imprinted genes is small and thus it is difficult to formally assess the statistical significance of such spatial clusters. On the other hand, while most imprinted genes have higher expression on the paternal allele, there is some variation at chromosome level. We can test whether the proportion of paternally imprinted (or equivalently maternally imprinted) genes is uniform across all chromosomes. In order to have enough genes for testing, we define potential imprinted genes at very liberal q-value cutoffs 0.25 or 0.5, and only consider the chromosomes with at least 2 imprinted genes. At q-value cutoff and , the p-values for Fisher’s exact test are and , respectively, suggesting that genes within one chromosome have more consistent pattern of parent-of-origin expression than expected by chance. Next, for each chromosome, we assess whether the proportion of paternally imprinted genes of this chromosome is different from all the other chromosomes. We performed a test only if there are at least 5 imprinted genes in a chromosome (Table 3). These formal tests confirm the clusters at chromosome 16 and 15, at imprinting q-value cutoff 0.25 and 0.5, respectively; although at the q-value cutoff 0.5, the p-value for chromosome 15 is not significant after multiple testing correction.
We developed a systematic approach to jointly estimate genetic and parent-of-origin effect in human. Our results recovered about one third of known imprinted genes, and if we excluded genes with low expression and included weaker, but non-contradicting cases, there was a good overall consistency. None of the genes classified as “predicted imprinting” instead of “imprinting” in Geneimprint database were detected as high confidence in our short list. One possible reason is that most of these genes were selected based on a screening paper (Luedi et al., 2007) and were false positives in that work. We also noted that for these “predicted imprinting” genes, the proportion of genes with predominantly paternal expression roughly 50%, while this proportion is about 68% (61/90) for known imprinting genes. Several imprinted genes that we found are non-coding RNAs, such as RP11-701H24.4 and RP1-309F20.3, which warrants further studies to elucidate the functional consequence if imprinted non-coding RNAs.
The raw RNA-seq data were submitted to NCI SRA database. The R package of our method will be submitted to CRAN and posted at http://research.fhcrc.org/sun/en/software.html.
Web Appendix, referenced in Section 2-4, is attached below.
6 RNA-seq data of 30 children of HapMap family trios
6.1 Data collection and mapping
RNA-seq data were collected from 30 children of HapMap family trios (15 males + 15 females) in five runs: 10 individuals sequenced in November 2014 with 150 bp paired-end reads (batch 1), 10 sequenced twice in December 2014 with 150 bp paired-end reads (batch 2), and 10 sequenced twice in January 2015 with most of the reads sequenced with 150bp paired-end reads and another smaller run with 75 bp paired-end reads (batch 3). RNA-seq reads were mapped to hg38 reference assembly using TopHat v2 Trapnell
et al. (2009), filtered with the following criteria: 3 read mismatches, read gap length 3, read edit distance 3, and read realign edit distance equals to 0.
6.2 Quality control
We further filtered RNA-seq reads with average base quality 30, mapping quality 20, and keeping only those uniquely mapped reads, by the prepareBAM function from R package asSeq. In a typical sample, most of the RNA-seq reads pass these filters. Figure 5 illustrates mean base quality for the first batch (November 2014).
6.3 Matching genotype coordinates to hg19
During phasing and imputation we worked with the reference panel mapped to hg19. Therefore the available genotype data was mapped to hg19 coordinates using liftover tool Rhead et al. (2009). After phasing and imputation was done, the resulting SNPs were mapped to hg38 coordinates, because the RNAseq data were mapped to hg38.
For this phasing step and the next imputation step, we used the 1000 Genome reference panel Howie
et al. (2011) (as of summer 2015) containing 2,504 individuals with 82 million SNPs. For the 90 individuals from 30 family trios, phasing was done using shapeit v.2. Family information of the 30 trios was used to improve the phasing accuracy. We started with 3.9 million SNPs of these 30 trios obtained from HapMap project, 1.3 million of which were filtered out due to following reasons: 30% of them were absent in the reference panel, 70% had mismatching information. Effective size of the population was set to the suggested value --effective-size 20000 and random seed was set to .
The results of phasing can be directly used as input for impute2, with imputation done using the latest reference available at impute2 web page Howie
et al. (2009). Imputation was done by splitting the genome into blocks of 5 Mb (no more than 7 Mb according to the instruction). We also used the same population size option as the one used in phasing step (-Ne 20,000), other options used: -align_by_maf_g and -seed 12345.
6.6 Creating the lists of heterozygous SNPs for each individual
We obtained the list of heterozygous SNPs for each sample from the results of imputation. The main output file (other files have additional extensions to the name as _haps, _allele_probs etc) has 3 columns with genotype probabilities per SNP. This file represents probability of observing respectively. We selected heterozygous locations (i.e. with high probability of ) to output phased genotypes of these locations.
6.7 Finding allele specific reads
Once we got a list of heterozygous SNPs, we extracted allele-specific reads for two haplotypes using R function extractASReads from R package asSeq Sun (2012). The extracted reads were saved in two additional bam files. We assigned the parent of origin of the two haplotypes based on how well they match the paternal and maternal haplotypes.
6.8 Gene level summarization
Finally, we summarized main bam file as well as the parent-specific bam files using R function summarizeOverlaps from R package GenomicAlignments Lawrence et al. (2013) with option “mode = Union” and “yieldSize = 1,000,000”. Across all the 30 samples, the percentage of reads being allele-specific are consistent (Figure 6). Approximately of the reads can be classified as allele specific. Comparing with percent of allele-specific counts from unrelated population we conclude that having a family structure allowed to increase the percentage of allele-specific reads by about . Additionally we can see that the proportion of reads with conflicting haplotype information is quite low in all samples (Table 4).
6.9 Processing of 462 lymphoblastoid cell lines from the 1000 Genomes samples
A similar process as in previous subsections was applied on RNA-seq data of 462 lymphoblastoid cell lines from the 1000 Genomes Lappalainen
et al. (2013): unrelated human lymphoblastoid cell line samples from the CEU, FIN, GBR, TSI and YRI populations (see E-GEUV-1 data link below) were sequenced with Illumina HiSeq2000 platform, with paired-end 75-bp reads. We mapped these reads to the hg38 reference with the same options as for original 30 individuals. To produce allele specific reads we started with 2,123 individuals of African or European descent with 2.2 million SNPs per individual produced by the Omni Array data originally from Delaneau
et al. (2012). Since African sub-population is quite different from Western European sub-population, we removed this sub-population and also removed samples, for which we observed many genes having extreme proportions of reads attributed to one haplotype - this is likely due to the fact of having a notable proportion of SNPs with problematic genotypes. Figure 7 shows typical samples we classify bad/good samples. A principal component analysis of total expression for the remaining 227 European samples suggested that western European individuals (CEU), which are the most comparable to our 30 samples, are also comparable to other sub-European populations. Thus we used these 227 samples for eQTL mapping.
7 eQTL analysis for 227 Caucasian samples from 1000 genome project
7.1 Comparing 227 Caucasian samples versus our 30 HapMap samples
We mapped cis-eQTLs using TReCASE model for the 227 Caucasian samples as well as our 30 samples (30 children of family trios). Figure 8 shows the location distribution of these cis-eQTLs relative to the transcription starting site of the corresponding genes. The distributions derived from the 227 Caucasian samples and our 30 samples are very similar.
In addition, the estimates of genetic effects using 227 individuals vs. 30 children at the same SNP (considering only SNPs declared significant at were also consistent (Figure 9).
8 Additional simulations
We define one model fit as the following. First fit a TReC only model with to obtain reasonable initial estimate of , and then fitting full TReCASE model estimating and , as well as two short models for additive genetic effect or parent-of-origin effect only. The computational time needed for one model fit (i.e., for one gene) is linearly dependent on sample size and increases for parameters being farther away from zero (Figure 10). For the sample size of around 200 the computation time is about 1 second per gene.
9 Software and Data
RNA-seq data from 1000G samples (E-GEUV-1 data) Lappalainen
et al. (2013) is available from
Shapeit v2 Delaneau et al. (2014): https://mathgen.stats.ox.ac.uk/genetics_software/shapeit/shapeit.html
IMPUTE v2 Howie
et al. (2009): http://www.stats.ox.ac.uk/~marchini/software/gwas/gwas.html
1000 Genome Reference Panel Howie
et al. (2011): https://mathgen.stats.ox.ac.uk/impute/impute_v2.html#reference
R/GenomicAlignments Lawrence et al. (2013): https://bioconductor.org/packages/release/bioc/html/GenomicAlignments.html
R/asSeq Sun (2012): http://research.fhcrc.org/sun/en/software/asSeq.html
This research is supported in part by NIH grant GM105785.
- 1000 Genomes Project Consortium (2012). An integrated map of genetic variation from 1,092 human genomes. Nature 491, 56–65.
- Barbaux, S., Gascoin-Lachambre, G., Buffat, C., Monnier, P., Mondon, F., Tonanny, M.-B., Pinard, A., Auer, J., Bessières, B., Barlier, A., et al. (2012). A genome-wide approach reveals novel imprinted genes expressed in the human placenta. Epigenetics 7, 1079–1090.
- Barlow, D. P. (2011). Genomic imprinting: a mammalian epigenetic discovery model. Annual review of genetics 45, 379–403.
- Crowley, J. J., Zhabotynsky, V., Sun, W., Huang, S., Pakatci, I. K., Kim, Y., Wang, J. R., Morgan, A. P., Calaway, J. D., Aylor, D. L., et al. (2015). Analyses of allele-specific gene expression in highly divergent mouse crosses identifies pervasive allelic imbalance. Nature genetics 47, 353–360.
- Delaneau, O., Marchini, J., Consortium, . G. P., et al. (2014). Integrating sequence and array data to create an improved 1000 genomes project haplotype reference panel. Nature communications 5,.
- Delaneau, O., Marchini, J., and Zagury, J.-F. (2012). A linear complexity phasing method for thousands of genomes. Nature methods 9, 179–181.
- Howie, B., Fuchsberger, C., Stephens, M., Marchini, J., and Abecasis, G. R. (2012). Fast and accurate genotype imputation in genome-wide association studies through pre-phasing. Nature genetics 44, 955–959.
- Howie, B., Marchini, J., and Stephens, M. (2011). Genotype imputation with thousands of genomes. G3: Genes, Genomes, Genetics 1, 457–470.
- Howie, B. N., Donnelly, P., and Marchini, J. (2009). A flexible and accurate genotype imputation method for the next generation of genome-wide association studies. PLoS Genet 5, e1000529.
- Hu, Y.-J., Sun, W., Tzeng, J.-Y., and Perou, C. M. (2015). Proper use of allele-specific expression improves statistical power for cis-eQTL mapping with RNA-seq data. Journal of the American Statistical Association 110, 962–974.
- Jirtle, R. L. (2016). Gene imprint, imprinted gene database (internet)Available at http://www.geneimprint.com/site/genes-by-species.
- Lappalainen, T., Sammeth, M., Friedländer, M. R., ACât Hoen, P., Monlong, J., Rivas, M. A., Gonzàlez-Porta, M., Kurbatova, N., Griebel, T., Ferreira, P. G., et al. (2013). Transcriptome and genome sequencing uncovers functional variation in humans. Nature 501, 506–511.
- Lawrence, M., Huber, W., Pages, H., Aboyoun, P., Carlson, M., Gentleman, R., Morgan, M. T., and Carey, V. J. (2013). Software for computing and annotating genomic ranges. PLoS Comput Biol 9, e1003118.
- Luedi, P. P., Dietrich, F. S., Weidman, J. R., Bosko, J. M., Jirtle, R. L., and Hartemink, A. J. (2007). Computational and experimental identification of novel human imprinted genes. Genome research 17, 1723–1730.
- Morcos, L., Ge, B., Koka, V., Lam, K. C., Pokholok, D. K., Gunderson, K. L., Montpetit, A., Verlaan, D. J., and Pastinen, T. (2011). Genome-wide assessment of imprinted expression in human cells. Genome biology 12, 1.
- Morison, I. M., Ramsay, J. P., and Spencer, H. G. (2005). A census of mammalian imprinting. TRENDS in Genetics 21, 457–465. http://www.otago.ac.nz/IGC.
- Peters, J. (2014). The role of genomic imprinting in biology and disease: an expanding view. Nature Reviews Genetics 15, 517–530.
- Rhead, B., Karolchik, D., Kuhn, R. M., Hinrichs, A. S., Zweig, A. S., Fujita, P. A., Diekhans, M., Smith, K. E., Rosenbloom, K. R., Raney, B. J., et al. (2009). The ucsc genome browser database: update 2010. Nucleic acids research page gkp939.
- Sun, W. (2012). A statistical framework for eQTL mapping using RNA-seq data. Biometrics 68, 1–11.
- Sun, W. and Hu, Y. (2013). eqtl mapping using rna-seq data. Statistics in biosciences 5, 198–219.
- Trapnell, C., Pachter, L., and Salzberg, S. L. (2009). Tophat: discovering splice junctions with RNA-seq. Bioinformatics 25, 1105–1111.
- Zou, F., Sun, W., Crowley, J. J., Zhabotynsky, V., Sullivan, P. F., and de Villena, F. P.-M. (2014). A novel statistical approach for jointly analyzing RNA-seq data from f1 reciprocal crosses and inbred lines. Genetics 197, 389–399.