Daniel J. Gaffney

A survey of best practices for RNA-seq data analysis

RNA-sequencing (RNA-seq) has a wide variety of applications, but no single analysis pipeline can be used in all cases. We review all of the major steps in RNA-seq data analysis, including experimental design, quality control, read alignment, quantification of gene and transcript levels, visualization, differential gene expression, alternative splicing, functional analysis, gene fusion detection and eQTL mapping. We highlight the challenges associated with each step. We discuss the analysis of small RNAs and the integration of RNA-seq with other functional genomics techniques. Finally, we discuss the outlook for novel technologies that are changing the state of the art in transcriptomics.

DNase I sensitivity QTLs are a major determinant of human expression variation

The mapping of expression quantitative trait loci (eQTLs) has emerged as an important tool for linking genetic variation to changes in gene regulation. However, it remains difficult to identify the causal variants underlying eQTLs, and little is known about the regulatory mechanisms by which they act. Here we show that genetic variants that modify chromatin accessibility and transcription factor binding are a major mechanism through which genetic variation leads to gene expression differences among humans. We used DNase I sequencing to measure chromatin accessibility in 70 Yoruba lymphoblastoid cell lines, for which genome-wide genotypes and estimates of gene expression levels are also available. We obtained a total of 2.7 billion uniquely mapped DNase I-sequencing (DNase-seq) reads, which allowed us to produce genome-wide maps of chromatin accessibility for each individual. We identified 8,902 locations at which the DNase-seq read depth correlated significantly with genotype at a nearby single nucleotide polymorphism or insertion/deletion (false discovery rate = 10%). We call such variants ‘DNase I sensitivity quantitative trait loci’ (dsQTLs). We found that dsQTLs are strongly enriched within inferred transcription factor binding sites and are frequently associated with allele-specific changes in transcription factor binding. A substantial fraction (16%) of dsQTLs are also associated with variation in the expression levels of nearby genes (that is, these loci are also classified as eQTLs). Conversely, we estimate that as many as 55% of eQTL single nucleotide polymorphisms are also dsQTLs. Our observations indicate that dsQTLs are highly abundant in the human genome and are likely to be important contributors to phenotypic variation.

Accurate inference of transcription factor binding from DNA sequence and chromatin accessibility data

Accurate functional annotation of regulatory elements is essential for understanding global gene regulation. Here, we report a genome-wide map of 827,000 transcription factor binding sites in human lymphoblastoid cell lines, which is comprised of sites corresponding to 239 position weight matrices of known transcription factor binding motifs, and 49 novel sequence motifs. To generate this map, we developed a probabilistic framework that integrates cell- or tissue-specific experimental data such as histone modifications and DNase I cleavage patterns with genomic information such as gene annotation and evolutionary conservation. Comparison to empirical ChIP-seq data suggests that our method is highly accurate yet has the advantage of targeting many factors in a single assay. We anticipate that this approach will be a valuable tool for genome-wide studies of gene regulation in a wide variety of cell types or tissues under diverse conditions.

DNA Sequence-Dependent Compartmentalization and Silencing of Chromatin at the Nuclear Lamina

A large fraction of the mammalian genome is organized into inactive chromosomal domains along the nuclear lamina. The mechanism by which these lamina associated domains (LADs) are established remains to be elucidated. Using genomic repositioning assays, we show that LADs, spanning the developmentally regulated IgH and Cyp3a loci contain discrete DNA regions that associate chromatin with the nuclear lamina and repress gene activity in fibroblasts. Lamina interaction is established during mitosis and likely involves the localized recruitment of Lamin B during late anaphase. Fine-scale mapping of LADs reveals numerous lamina-associating sequences (LASs), which are enriched for a GAGA motif. This repeated motif directs lamina association and is bound by the transcriptional repressor cKrox, in a complex with HDAC3 and Lap2β. Knockdown of cKrox or HDAC3 results in dissociation of LASs/LADs from the nuclear lamina. These results reveal a mechanism that couples nuclear compartmentalization of chromatin domains with the control of gene activity.

Dense fine-mapping study identifies new susceptibility loci for primary biliary cirrhosis

We genotyped 2,861 cases of primary biliary cirrhosis (PBC) from the UK PBC Consortium and 8,514 UK population controls across 196,524 variants within 186 known autoimmune risk loci. We identified 3 loci newly associated with PBC (at P < 5 × 10−8), increasing the number of known susceptibility loci to 25. The most associated variant at 19p12 is a low-frequency nonsynonymous SNP in TYK2, further implicating JAK-STAT and cytokine signaling in disease pathogenesis. An additional five loci contained nonsynonymous variants in high linkage disequilibrium (LD; r2 > 0.8) with the most associated variant at the locus. We found multiple independent common, low-frequency and rare variant association signals at five loci. Of the 26 independent non–human leukocyte antigen (HLA) signals tagged on the Immunochip, 15 have SNPs in B-lymphoblastoid open chromatin regions in high LD (r2 > 0.8) with the most associated variant. This study shows how data from dense fine-mapping arrays coupled with functional genomic data can be used to identify candidate causal variants for functional follow-up.

Controls of Nucleosome Positioning in the Human Genome

Nucleosomes are important for gene regulation because their arrangement on the genome can control which proteins bind to DNA. Currently, few human nucleosomes are thought to be consistently positioned across cells; however, this has been difficult to assess due to the limited resolution of existing data. We performed paired-end sequencing of micrococcal nuclease-digested chromatin (MNase–seq) from seven lymphoblastoid cell lines and mapped over 3.6 billion MNase–seq fragments to the human genome to create the highest-resolution map of nucleosome occupancy to date in a human cell type. In contrast to previous results, we find that most nucleosomes have more consistent positioning than expected by chance and a substantial fraction (8.7%) of nucleosomes have moderate to strong positioning. In aggregate, nucleosome sequences have 10 bp periodic patterns in dinucleotide frequency and DNase I sensitivity; and, across cells, nucleosomes frequently have translational offsets that are multiples of 10 bp. We estimate that almost half of the genome contains regularly spaced arrays of nucleosomes, which are enriched in active chromatin domains. Single nucleotide polymorphisms that reduce DNase I sensitivity can disrupt the phasing of nucleosome arrays, which indicates that they often result from positioning against a barrier formed by other proteins. However, nucleosome arrays can also be created by DNA sequence alone. The most striking example is an array of over 400 nucleosomes on chromosome 12 that is created by tandem repetition of sequences with strong positioning properties. In summary, a large fraction of nucleosomes are consistently positioned—in some regions because they adopt favored sequence positions, and in other regions because they are forced into specific arrangements by chromatin remodeling or DNA binding proteins.

Dissecting the regulatory architecture of gene expression QTLs

BackgroundExpression quantitative trait loci (eQTLs) are likely to play an important role in the genetics of complex traits; however, their functional basis remains poorly understood. Using the HapMap lymphoblastoid cell lines, we combine 1000 Genomes genotypes and an extensive catalogue of human functional elements to investigate the biological mechanisms that eQTLs perturb.ResultsWe use a Bayesian hierarchical model to estimate the enrichment of eQTLs in a wide variety of regulatory annotations. We find that approximately 40% of eQTLs occur in open chromatin, and that they are particularly enriched in transcription factor binding sites, suggesting that many directly impact protein-DNA interactions. Analysis of core promoter regions shows that eQTLs also frequently disrupt some known core promoter motifs but, surprisingly, are not enriched in other well-known motifs such as the TATA box. We also show that information from regulatory annotations alone, when weighted by the hierarchical model, can provide a meaningful ranking of the SNPs that are most likely to drive gene expression variation.ConclusionsOur study demonstrates how regulatory annotation and the association signal derived from eQTL-mapping can be combined into a single framework. We used this approach to further our understanding of the biology that drives human gene expression variation, and of the putatively causal SNPs that underlie it.

Genetic Background Drives Transcriptional Variation in Human Induced Pluripotent Stem Cells

Human iPS cells have been generated using a diverse range of tissues from a variety of donors using different reprogramming vectors. However, these cell lines are heterogeneous, which presents a limitation for their use in disease modeling and personalized medicine. To explore the basis of this heterogeneity we generated 25 iPS cell lines under normalised conditions from the same set of somatic tissues across a number of donors. RNA-seq data sets from each cell line were compared to identify the majority contributors to transcriptional heterogeneity. We found that genetic differences between individual donors were the major cause of transcriptional variation between lines. In contrast, residual signatures from the somatic cell of origin, so called epigenetic memory, contributed relatively little to transcriptional variation. Thus, underlying genetic background variation is responsible for most heterogeneity between human iPS cell lines. We conclude that epigenetic effects in hIPSCs are minimal, and that hIPSCs are a stable, robust and powerful platform for large-scale studies of the function of genetic differences between individuals. Our data also suggest that future studies using hIPSCs as a model system should focus most effort on collection of large numbers of donors, rather than generating large numbers of lines from the same donor.

Epithelial IL-22RA1-Mediated Fucosylation Promotes Intestinal Colonization Resistance to an Opportunistic Pathogen

Summary Our intestinal microbiota harbors a diverse microbial community, often containing opportunistic bacteria with virulence potential. However, mutualistic host-microbial interactions prevent disease by opportunistic pathogens through poorly understood mechanisms. We show that the epithelial interleukin-22 receptor IL-22RA1 protects against lethal Citrobacter rodentium infection and chemical-induced colitis by promoting colonization resistance against an intestinal opportunistic bacterium, Enterococcus faecalis. Susceptibility of Il22ra1−/− mice to C. rodentium was associated with preferential expansion and epithelial translocation of pathogenic E. faecalis during severe microbial dysbiosis and was ameloriated with antibiotics active against E. faecalis. RNA sequencing analyses of primary colonic organoids showed that IL-22RA1 signaling promotes intestinal fucosylation via induction of the fucosyltransferase Fut2. Additionally, administration of fucosylated oligosaccharides to C. rodentium-challenged Il22ra1−/− mice attenuated infection and promoted E. faecalis colonization resistance by restoring the diversity of anaerobic commensal symbionts. These results support a model whereby IL-22RA1 enhances host-microbiota mutualism to limit detrimental overcolonization by opportunistic pathogens.

Functional constraints and frequency of deleterious mutations in noncoding DNA of rodents

Selection against deleterious mutations imposes a mutation load on populations because individuals die or fail to reproduce. In vertebrates, estimates of genomic rates of deleterious mutations in protein-coding genes imply the existence of a substantial mutation load, but many functionally important regions of the genome are thought to reside in noncoding DNA, and the contribution of noncoding DNA to the mutation load has been unresolved. Here, we infer the frequency of deleterious mutations in noncoding DNA of rodents by comparing rates of substitution at noncoding nucleotides with rates of substitution at the fastest evolving intronic sites of adjacent genes sampled from the whole genome sequences of mouse and rat. We show that the major elements of selectively constrained noncoding DNA are within 2,500 bp upstream and downstream of coding sequences and in first introns. Our estimate of the genomic deleterious point mutation rate for noncoding DNA (0.22 per diploid per generation) is similar to that for coding DNA. Mammalian populations therefore experience a substantial genetic load associated with selection against deleterious mutations in noncoding DNA. Deleterious mutations in noncoding DNA have predominantly quantitative effects and could be an important source of the burden of complex genetic disease variation in human populations.