balanced approach where mismatches were allowed but there was high confidence in the mapping location for a fragment to be counted.) To evaluate the quality of the RNA-sequencing data, the numbers of reads that map to different annotated regions (e.g. exonic, intronic, splicing junction, intergenic, promoter, UTR) were determined with bamUtils (Breese & Liu 2013). Low quality mapped reads (including reads mapped to multiple positions) were excluded. The tool ‘featureCounts’ was used to quantify gene level read counts (Liao et al. 2014). Differential gene expression analysis was performed with edgeR (Robinson et al. 2010). According to recommendations in the edgeR user’s guide, genes with read counts <0.5 CPM in ≥ 12 samples were filtered out and a negative binomial generalized linear model with likelihood ratio test was used. A multi-dimensional scaling (MDS) plot provided a visual representation of the similarities and differences among the experimental samples (Supplementary Figure 1A,B). Distances on the plot corresponded to the Euclidian distance (root-mean-square) of leading log2 fold change (FC) for the most variable genes between each pair of samples (Ritchie et al. 2015). Comparisons were adjusted for sex by implementing a batch factor according to the edgeR manual. One sample, a HAP male, was
