Conditional Random Field for inferring Neandertal local ancestry

Consider a haploid genome in a test population that carries Neandertalancestry, for example Europeans. Given the allelic states of a sequence of SNPsalong this haplotype, we would like to infer the ancestral state of the alleleat each SNP, specifically, whether it has entered modern humans throughNeandertal gene flow. In addition to the test haplotype, the data we analyzeconsist of a panel of haplotypes from the sub-Saharan African Yoruba (YRI) whowe assume harbor no Neandertal ancestry¹. To determine the allelic state of the Neandertals, weuse a high-coverage Neandertal genome². We determine the ancestral and derived allele at eachSNP using a 6-primate consensus sequence¹³. To estimate the genetic distance between adjacentSNPs, we use the Oxford combined linkage disequilibrium (LD) map²⁹. We specify the distributionof the unobserved Neandertal ancestry states at each SNP given the observedgenetic data as a Conditional Random Field (CRF)⁹. Intuitively, we specify featurefunctions that relate the observed data and the unobserved ancestral state ateach SNP (“emission functions”) as well as feature functionsthat relate the unobserved ancestral states at adjacent SNPs(“transition functions”). Thus, the model is a linear-chain CRF.The feature functions and their associated parameters fully specify thedistribution of the unobserved ancestral states given the observed data. Giventhe parameters and the observed data, we are able to infer the marginalprobability of Neandertal ancestry at each SNP of the haploid genome. We computethe marginal probabilities efficiently using the forward-backward algorithm^9,30. SI 1 presents the mathematicaldetails.

Feature functions

The emission functions couple the unobserved ancestral state at a SNP tothe observed features. We use two classes of emission functions.

The first class of emission function captures information from the jointpatterns observed at a single SNP across Europeans, Africans and Neandertals.These features are indicator functions that assume the value “1”when a specific pattern is observed at a SNP and “0” otherwise.We use feature functions that pick out two classes of allelic patterns. One ofthese features is 1 if at a given SNP, the test haplotype carries the derivedallele, all the YRI haplotypes carry the ancestral allele, and either of the twoNeandertal alleles is derived. SNPs with this joint configuration have anincreased likelihood of Neandertal ancestry. In the CRF, an increased likelihoodassociated with this feature is reflected in the fact that the parameter ispositive with a magnitude determined by the informativeness of the feature. Thesecond feature is 1 if at a given SNP, the test haplotype carries a derivedallele that is polymorphic in the panel of Africans but absent in theNeandertal. SNPs with this joint configuration have a decreased likelihood ofNeandertal ancestry.

The second class of emission functions uses multiple SNPs to capture thesignal of Neandertal ancestry. Specifically, we compare the divergence of thetest haplotype to the Neandertal sequence to the minimum divergence of the testhaplotype to all African haplotypes over non-overlapping 100 kilobase (kb)windows (the size scale we expect for Neandertal haplotypes today based on thetime of Neandertal gene flow into modern humans¹⁰). In a region of the genome where thetest haplotype carries Neandertal ancestry, we expect the test haplotype to becloser to the Neandertal sequence than to most modern human sequences (albeitwith a large variance), and we expect the pattern to be reversed outside theseregions. While computing distance to the Neandertal sequence, we build aNeandertal haplotype by choosing the derived allele at heterozygous sites sothat this distance is effectively the minimum distance of the potentiallyintrogressed test haplotype to one of the two Neandertal haplotypes.

The transition feature function modulates the correlation of theancestral states at adjacent SNPs. We define this feature function as anapproximation, at small genetic distance, to the log of the transitionprobabilities of a standard Markov process of admixture between two populations.This approximation makes parameter estimation in the CRF efficient.

Parameter Estimation

To estimate the parameters of the CRF, we need haplotypes labeled withNeandertal ancestries; that is, training data. Since we do not in fact know thetrue Neandertal state in any individual, we estimated the CRF parameters on datasimulated under a demographic model. We estimate parameters by maximizing theL2-regularized conditional log likelihood using a limited-memory version ofLBFGS³¹. We fixed the valueof the parameter associated with the L2-penalty at 10 although a broad range ofvalues appear to work in practice. We assumed a simple demographic modelrelating Africans, Europeans and Neandertals with Neandertal-modern humanadmixture occurring 1,900 generations ago¹⁰ and a fraction of Neandertal ancestry of 3%¹. The model parameterswere broadly constrained by the observed allele frequency differentiationbetween the West African YRI and European American CEU populations and by theobserved excess sharing of alleles between Europeans and Africans relative toNeandertals. The simulations incorporated hotspots of recombination¹¹ as well as the reduced powerto detect low-frequency alleles from low-coverage sequencing data¹².

Validation of the CRF

We assessed the accuracy of the CRF to predict Neandertal ancestry usingsimulated data. Given the marginal probabilities estimated by the CRF, weestimated the precision (fraction of predictions that are truly Neandertal) andthe recall (fraction of true Neandertal alleles that are predicted) as we varythe threshold on the marginal probability for an allele to be declaredNeandertal. We also evaluated the accuracy when the haplotype phase needs to beinferred and when the genetic map has errors. Since the CRF parameter estimationassumes a specific demographic model, we were concerned about the possibilitythat the inferences might be sensitive to the model assumed. We thereforeperturbed each demographic parameter in turn and applied the CRF to datasimulated under these perturbed models, fixing the parameters of the CRF to theestimates obtained under the original model. For each of these perturbed models,we evaluated the false discovery rate (1-precision) when we restrict to sites atwhich the CRF assigns a marginal probability of at least 0.90. SI 2 presents thedetails.

Preparation of 1000 Genomes Data

We applied the CRF to the computationally phased haplotypes in each ofthe 13 populations in the 1000 Genomes project¹² (1KG), excluding the west African Yoruba(YRI). The CRF requires reference genomes from Africans and Neandertals. For theAfrican population, we used 176 haplotypes from 88 YRI individuals. For theNeandertal genome, we used the genotypes called from the recently generatedhigh-coverage Neandertal sequence². We restricted our analysis to sites passing the filtersdescribed in ref.² and forwhich the genotype quality score GQ ≥ 30. These filters discard sitesthat are identified as repeats by the Tandem Repeat Finder³², that have Phred-scaled mapping qualityscores of MQ < 30, or that map to regions where the alignment is ambiguous orwhich fall within the upper or lower 2.5^th percentile of thesample-specific coverage distribution (applied within the regions of uniquemappability binned according to the GC-content of the reference genome). For themappability filter, we used the liberal filter that requires that at least50% of all 35-mers that overlap a position do not map to any otherposition in the genome allowing up to one mismatch. We further restricted ouranalysis to sites that are biallelic across the Neandertal and the 1000 Genomesproject samples. For each haplotype analyzed, we also restricted to the set ofpolymorphic sites in the population containing the haplotype. After filtering,we were able to analyze 26,493,206 SNPs on the autosomes and 817,447 SNPs onchromosome X. We obtained genetic distances from the Oxford combined LD maplifted over to hg19 coordinates²⁹. For the X chromosome, we obtained an appropriatesex-averaged map by scaling the X chromosome LD-based map by 2/3.

Statistics for measuring Neandertal ancestry

We computed several statistics to summarize the Neandertal ancestryinferred by the CRF. We estimated the proportion of an individual diploid genomethat is confidently inferred to be Neandertal as the fraction of sites for whichthe marginal probability is ≥ 90%. To assess variation in theproportion of Neandertal ancestry along the genome, we computed the fraction ofalleles across individuals with marginal probability greater than a specifiedthreshold. This statistic is likely to be affected by variation in power alongthe genome. Hence, we also consider an estimate of the ancestry proportionobtained by averaging the marginal probability across individuals. Depending onthe analyses, these statistics are estimated at a single SNP or innon-overlapping windows of a specified size.

To assess whether the predictions made by the CRF are sensible, weinferred Neandertal ancestry using the low-coverage genome from the VindijaNeandertals¹. For thisanalysis, we restricted to sites at which there is at least one read withmapping quality score between 60 and 90 and base quality of at least 40. As asecond validation analysis, we applied the CRF to the sub-Saharan African Luhya(LWK), using the parameters optimized for non-Africans. We empirically assessedthe accuracy of the CRF on the 1000 Genomes project data by assuming that LWKhas no Neandertal ancestry, that the false discovery rate in each non-Africanpopulation is equal to the false discovery rate in LWK, and using thegenome-wide proportion of Neandertal ancestry estimated in ref.² (SI 3). We computed thetheoretical standard deviation in the proportion of Neandertal ancestry³³ assuming a pulsemodel of admixture with 2% Neandertal ancestry followed by 2,000generations of random mating, and 2.03 Gigabases as the number of bases of thehigh-coverage Neandertal genome that pass filters.

Tiling path of Neandertal haplotypes

We identified Neandertal haplotypes as runs of consecutive alleles alonga test haplotype assigned a marginal probability of > 90%. Wefiltered haplotypes smaller than 0.02 cM. At each SNP that is covered by atleast one such haplotype, we estimated the allelic state as the consensus alleleacross the spanning haplotypes. See SI 4.

Functional analysis of introgressed alleles

We defined two subsets of Consensus Coding Sequence (CCDS) genes³⁴. We define a genewith “low Neandertal ancestry” as one in which all allelesacross all individuals have a marginal probability ≤ 10%. Wealso require that the genes included in this analysis include at least ≥100 SNPs within a 100 kb window centered at its midpoint (this excluded geneswith low power). We define a gene with “high Neandertalancestry” as one that is in the top 5% of CCDS genes ranked bythe average marginal probability across individual haplotypes.

Functional enrichment analysis

We tested for enrichment of Gene Ontology (GO)³⁵ categories in genes with low or highNeandertal ancestry, using the hypergeometric test implemented in the FUNCpackage³⁶. We reportmultiple-testing corrected P-values estimated from 1000 permutations for the GOenrichment analysis (Family-Wise Error Rate – FWER). Given the observedcorrelation between Neandertal ancestry and B-statistic¹⁸, a concern is that the functionalcategories may not be randomly distributed with respect to B-statistic. Tocontrol for this, we assigned a B-statistic to each gene (estimated as anaverage of the B-statistic over the length of the gene) as well as a uniformrandom number. This resulted in 17,249 autosomal genes. Genes were binned into20-equal sized bins based on the gene-specific B-statistic. Within each bin,genes were sorted by Neandertal ancestry and then by the random number. Genesranked within the top 5% within each bin were used for the analysis. SeeSI 6 for details.

Identifying alleles born in Neandertals and cross-correlation withassociation study data

To infer whether an allele segregating in a present-day human populationwas introduced by Neandertal gene flow, we defined a probable Neandertal alleleas one with marginal probability of ≥ 90% and a non-Neandertalallele as having a marginal probability of ≤ 10%. A SNP at whichall of the confident non-Neandertal alleles as well as all alleles in YRI areancestral and all of the confident Neandertal alleles are derived is inferred tobe of Neandertal origin. This allows for some false negatives in the predictionof the CRF. This procedure yields 97,365 Neandertal-derived SNPs when applied tothe predictions in Europeans and East Asians. We downloaded the variants listedin the NHGRI GWAS catalog¹⁶,retaining entries for which the reported association is a SNP with an assignedrs-number, and where the nominal P-value is<5×10⁻⁸. This resulted in 5,022 associations,which we then intersected with the Neandertal-derived list. See SI 7 fordetails.

Identification of genomic regions deficient in Neandertal ancestry

We measured the fraction of alleles across individuals and SNPs thathave less than a specified proportion of Neandertal ancestry measured withinnon-overlapping 10 Megabase (Mb) windows. We chose a threshold of marginalprobability of >25% as this threshold was found to lead to highrecall in our empirical assessment (SI 3). We reported windows for which thisstatistic is < 0.1%.

To understand the causes for variation in Neandertal ancestry, we testedfor correlation to a B-statistic. Each SNP was annotated with the B-statisticlifted over to hg19 coordinates. We assessed correlation between B and estimatesof Neandertal ancestry proportion at a nucleotide level as well as at differentsize scales, and separately on the autosomes and chromosome X (SI 8). We alsoused wavelet decomposition³⁷to analyze the correlation of the inferred Neandertal ancestry between Europeansand East Asians at multiple size scales (SI 10).Figure 2 reports the relation between the mean marginal probabilityof Neandertal ancestry across individuals and quintiles of B-statistic at eachSNP. To assess significance, we estimated Spearman’s correlationρ_Spearman and standard errors using a block jackknife³⁸ with 10 Mb blocks(SI 8).

To understand the contribution of demography to variation in Neandertalancestry along the genome, we measured the coefficient of variation, at a 10 Mbsize scale, of the proportion of ancestry estimated as defined above. We thenapplied the CRF to data simulated under diverse demographic models and comparedthe coefficient of variation to the observed value (SI 8).

Power to infer Neandertal ancestry as a function of demography and genomicfeatures

To assess the power of the CRF to infer Neandertal ancestry, wesimulated data under diverse demographic models³⁹. In one simulation series, we variedeffective population size to approximate the effect of background selection andmeasured recall at a precision of 90%. In a second series, we assessedpower on chromosome X versus the autosomes by matching the effective populationsize, recombination rate and mutation rate to estimated values for chromosome X.See SI 2.

Unbiased estimate of the proportion of Neandertal ancestry as a function ofB-statistic

To estimate the proportion of Neandertal ancestry in an unbiased way, wedivided the genome into quintiles of B, and estimated the proportion ofNeandertal ancestry using a statistic first published in ref.¹⁹. This statistic measures howmuch closer a non-African individual is to Denisova than an African individual,divided by the same quantity replacing the non-African individual withNeandertal. We report the estimated proportion of Neandertal ancestry in eachquintile as a fraction of the genome-wide mean and obtain standard errors usinga block jackknife with 100 blocks.

We analyzed data from 27 deeply sequenced genomes: 25 present-day humansand the high-coverage Neandertal and Denisova¹⁴ genomes. For each, we required thatsites pass the more stringent set of the two filters described in ref.², have a genotypequality of GQ ≥ 45, and have an ancestral allele that can be determinedbased on comparison to chimpanzee and at least one of gorilla or orangutan. Wecomputed a Z-score for the difference in the ancestry across the bin of highestB-statistic versus the rest and used a Bonferroni correction for ten hypotheses(5 hypotheses based on which set of bins we merge and a 2-sided test in each).In our main analysis, we analyzed both transitions and transversions and pooledgenomes for all non-African samples. We also analyzed other subsets of the data:transversions only in non-Africans, in Europeans and in eastern non-Africans.See SI 9 for details.

Tissue-specific expression

We defined tissue-specific expression levels using the Illumina BodyMap2.0 RNA-seq data²⁵, whichcontains expression data from 16 human tissues. We identified genes that areexpressed in a tissue-specific manner using the DESeq package⁴⁰ and used a P-value cutoff of0.05. We tested enrichment of tissue-specific genes in regions of high or lowNeandertal ancestry. A concern when testing for enrichment is that clustering ofsimilarly expressed genes coupled with the large size of regions of lowNeandertal ancestry might lead to spurious signals of enrichment. Hence, wedevised a permutation test that randomly rotated the annotations of genes(treating each chromosome as a circle) while maintaining the correlation withingenes and within Neandertal ancestry as well as between Neandertal ancestry andgenes. We tested enrichment on the whole genome, on the autosomes alone, and onchromosome X alone. We generated 1,000 random rotations for each test except forthe X chromosome for which we generated all possible rotations. We computed thefraction of permutations for which the P-value of Fischer’s exact testis at least as low as the observed P-value (See SI 6).

• 1.Green RE, et al. A draft sequence of the Neandertal genome. Science. 2010;328:710–722. doi: 10.1126/science.1188021. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 2.Prüfer K, et al. The complete genome sequence of a Neandertal from the AltaiMountains. Nature. 2014;505:43–49. doi: 10.1038/nature12886. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 3.Abi-Rached L, et al. The shaping of modern human immune systems by multiregionaladmixture with archaic humans. Science. 2011;334:89–94. doi: 10.1126/science.1209202. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 4.Mendez FL, Watkins JC, Hammer MF. A haplotype at STAT2 Introgressed from neanderthals and serves asa candidate of positive selection in Papua New Guinea. Am J Hum Genet. 2012;91:265–274. doi: 10.1016/j.ajhg.2012.06.015. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 5.Mendez FL, Watkins JC, Hammer MF. Neandertal origin of genetic variation at the cluster of OASimmunity genes. Mol Biol Evol. 2013;30:798–801. doi: 10.1093/molbev/mst004. [DOI] [PubMed] [Google Scholar]
• 6.Yotova V, et al. An X-linked haplotype of Neandertal origin is present among allnon-African populations. Mol Biol Evol. 2011;28:1957–1962. doi: 10.1093/molbev/msr024. [DOI] [PubMed] [Google Scholar]
• 7.Wall JD, et al. Higher levels of neanderthal ancestry in East aSians than inEuropeans. Genetics. 2013;194:199–209. doi: 10.1534/genetics.112.148213. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 8.Lachance J, et al. Evolutionary history and adaptation from high-coveragewhole-genome sequences of diverse African hunter-gatherers. Cell. 2012;150:457–469. doi: 10.1016/j.cell.2012.07.009. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 9.Lafferty J, McCallum A, Pereira FCN. Conditional Random Fields: Probabilistic Models for Segmentingand Labeling Sequence Data. ICML ‘01 Proceedings of the Eighteenth InternationalConference on Machine Learning; 2001. pp. 282–289. [Google Scholar]
• 10.Sankararaman S, Patterson N, Li H, Paabo S, Reich D. The date of interbreeding between Neandertals and modernhumans. PLoS Genet. 2012;8:e1002947. doi: 10.1371/journal.pgen.1002947. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 11.Hellenthal G, Stephens M. msHOT: modifying Hudson’s ms simulator to incorporatecrossover and gene conversion hotspots. Bioinformatics. 2007;23:520–521. doi: 10.1093/bioinformatics/btl622. [DOI] [PubMed] [Google Scholar]
• 12.Abecasis GR, et al. An integrated map of genetic variation from 1,092 humangenomes. Nature. 2012;491:56–65. doi: 10.1038/nature11632. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 13.Paten B, Herrero J, Beal K, Fitzgerald S, Birney E. Enredo and Pecan: genome-wide mammalian consistency-basedmultiple alignment with paralogs. Genome Res. 2008;18:1814–1828. doi: 10.1101/gr.076554.108. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 14.Meyer M, et al. A high-coverage genome sequence from an archaic Denisovanindividual. Science. 2012;338:222–226. doi: 10.1126/science.1224344. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 15.23andMe White Paper 23-05. Neanderthal Ancestry Estimator. <https://23andme.https.internapcdn.net/res/pdf/hXitekfSJe1lcIy7-Q72XA_23-05_Neanderthal_Ancestry.pdf>.
• 16.Hindorff LA, et al. Potential etiologic and functional implications of genome-wideassociation loci for human diseases and traits. Proc Natl Acad Sci U S A. 2009;106:9362–9367. doi: 10.1073/pnas.0903103106. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 17.The SIGMA Type 2 Diabetes Consortium. Sequence variants in SLC16A11 are a common risk factor for type 2diabetes in Mexico. Nature. 2014;506(7486):97–101. doi: 10.1038/nature12828. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 18.McVicker G, Gordon D, Davis C, Green P. Widespread genomic signatures of natural selection in hominidevolution. PLoS Genet. 2009;5:e1000471. doi: 10.1371/journal.pgen.1000471. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 19.Reich D, et al. Genetic history of an archaic hominin group from Denisova Cave inSiberia. Nature. 2010;468:1053–1060. doi: 10.1038/nature09710. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 20.Tucker PK, Sage RD, Wilson AC, Eichler EM. Abrupt cline for sex chromosomes in a hybrid zone between twospecies of mice. Evolution Int J Org Evolution. 1992;46:1146–1163. doi: 10.1111/j.1558-5646.1992.tb00625.x. [DOI] [PubMed] [Google Scholar]
• 21.Good JM, Dean MD, Nachman MW. A complex genetic basis to X-linked hybrid male sterility betweentwo species of house mice. Genetics. 2008;179:2213–2228. doi: 10.1534/genetics.107.085340. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 22.Presgraves DC. Sex chromosomes and speciation in Drosophila. Trends Genet. 2008;24:336–343. doi: 10.1016/j.tig.2008.04.007. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 23.Coyne JAO, HA . In: Speciation and its Consequences. Otte D, Endler JA, editors. 1989. pp. 180–207. [Google Scholar]
• 24.Wu CI, Davis AW. Evolution of postmating reproductive isolation: the compositenature of Haldane’s rule and its genetic basis. American Naturalist. 1993;142:187–212. doi: 10.1086/285534. [DOI] [PubMed] [Google Scholar]
• 25.Derrien T, et al. The GENCODE v7 catalog of human long noncoding RNAs: analysis oftheir gene structure, evolution, and expression. Genome Res. 2012;22:1775–1789. doi: 10.1101/gr.132159.111. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 26.Keinan A, Mullikin JC, Patterson N, Reich D. Measurement of the human allele frequency spectrum demonstratesgreater genetic drift in East Asians than in Europeans. Nat Genet. 2007;39:1251–1255. doi: 10.1038/ng2116. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 27.Bhatia G, et al. Genome-wide scan of 29,141 African Americans finds no evidence ofselection since admixture. arXiv:1312.2675 [q-bio.PE] [DOI] [PMC free article] [PubMed] [Google Scholar]
• 28.Orr HA, Turelli M. The evolution of postzygotic isolation: accumulatingDobzhansky-Muller incompatibilities. Evolution. 2001;55:1085–1094. doi: 10.1111/j.0014-3820.2001.tb00628.x. [DOI] [PubMed] [Google Scholar]
• 29.Myers S, Bottolo L, Freeman C, McVean G, Donnelly P. A fine-scale map of recombination rates and hotspots across thehuman genome. Science. 2005;310:321–324. doi: 10.1126/science.1117196. [DOI] [PubMed] [Google Scholar]
• 30.Sutton C, McCallum A. In: Introduction to Statistical Relational Learning. Getoor Lise, Taskar Ben., editors. Ch. 4. MIT Press; 2007. pp. 93–128. [Google Scholar]
• 31.Byrd RH, Nocedal J, Schnabel RB. Representations of quasi-Newton matrices and their use in limitedmemory methods. Mathematical Programming. 1994;63.1-3:129–156. [Google Scholar]
• 32.Tandem Repeat Finder. <url={http://hgdownload.soe.ucsc.edu/goldenPath/hg19/database/simpleRepeat.txt.gz}>.
• 33.Gravel S. Population genetics models of local ancestry. Genetics. 2012;191:607–619. doi: 10.1534/genetics.112.139808. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 34.Pruitt KD, et al. The consensus coding sequence (CCDS) project: Identifying acommon protein-coding gene set for the human and mousegenomes. Genome Res. 2009;19:1316–1323. doi: 10.1101/gr.080531.108. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 35.Ashburner M, et al. Gene ontology: tool for the unification of biology. The GeneOntology Consortium. Nat Genet. 2000;25:25–29. doi: 10.1038/75556. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 36.Prufer K, et al. FUNC: a package for detecting significant associations betweengene sets and ontological annotations. BMC Bioinformatics. 2007;8:41. doi: 10.1186/1471-2105-8-41. [DOI] [PMC free article] [PubMed] [Google Scholar]
• 37.Percival DB, Walden AT. Wavelet Methods for Time Series Analysis. Cambridge University Press; 2005. [Google Scholar]
• 38.Kunsch HR. The jackknife and the bootstrap for general stationaryobservations. The Annals of Statistics. 1989;17:1217–1241. [Google Scholar]
• 39.Hudson RR. Generating samples under a Wright-Fisher neutral model of geneticvariation. Bioinformatics. 2002;18:337–338. doi: 10.1093/bioinformatics/18.2.337. [DOI] [PubMed] [Google Scholar]
• 40.Anders S, Huber W. Differential expression analysis for sequence countdata. Genome Biol. 2010;11:R106. doi: 10.1186/gb-2010-11-10-r106. [DOI] [PMC free article] [PubMed] [Google Scholar]

Supplementary Materials