MW 12:50-2:05pm in Beckman B302 Profs: Serafim Batzoglou & Gill Bejerano TAs: Harendra Guturu & Panos A

1. CS273A Lecture 16: Functional Genomics MW  12:50-2:05pm in Beckman B302 Profs: Serafim Batzoglou & Gill Bejerano TAs: HarendraGuturu & PanosAchlioptas http://cs273a.stanford.edu [BejeranoFall13/14]

2. Gene set enrichment analysis: The gene regulatory version http://cs273a.stanford.edu [BejeranoFall13/14]

3. Cluster all genes for differential expression Experiment Control (replicates) (replicates) Most significantly up-regulated genes genes Unchanged genes Most significantly down-regulated genes http://cs273a.stanford.edu [BejeranoFall13/14]

4. The Gene Sets to test come from the Ontologies TJL-2004

5. Gene set 3 up regulated Gene set 2 down regulated Ask about whole gene sets Gene Set 1 Gene Set 2 Gene Set 3 Exper. Control + ES/NES statistic - http://cs273a.stanford.edu [BejeranoFall13/14]

6. Combinatorial Regulatory Code 2,000 different proteins can bind specific DNA sequences. Proteins DNA Protein binding site Gene DNA A regulatory region encodes 3-10 such protein binding sites. When all are bound by proteins the regulatory region turns “on”,and the nearby gene is activated to produce protein. http://cs273a.stanford.edu [BejeranoFall13/14]

7. ChIP-Seq: first glimpses of the regulatory genome in action Peak Calling Cis-regulatory peak http://cs273a.stanford.edu [BejeranoFall13/14] 7

8. What is the transcription factor I just assayed doing? • Collect known literature of the form • Function A: Gene1, Gene2, Gene3, ... • Function B: Gene1, Gene2, Gene3, ... • Function C: ... • Ask whether the binding sites you discovered are preferentially binding (regulating) any one or more of the functions listed above. • Form hypothesis and perform further experiments. Gene transcription start site Cis-regulatory peak http://cs273a.stanford.edu [BejeranoFall13/14] 8

9. Example: inferring functions of Serum Response Factor (SRF) from its ChIP-seq binding profile Gene transcription start site SRF binding ChIP-seq peak • ChIP-seq identified 2,429 SRF binding peaks in human Jurkat cells1 • SRF is known as a “master regulator of the actin cytoskeleton” • In the ChIP-Seq peaks, we expect to find binding sites regulating (genes involved in) actin cytoskeleton formation. http://cs273a.stanford.edu [BejeranoFall13/14]

10. Example: inferring functions of Serum Response Factor (SRF) from its ChIP-seq binding profile Gene transcription start site SRF binding ChIP-seq peak Ontology term (e.g. ‘actin cytoskeleton’) π π π π • Existing, gene-based method to analyze enrichment: • Ignore distal binding events. • Count affected genes. • Rank by enrichment hypergeometric p-value. N = 8 genes in genome K = 3 genes annotated with n = 2 genes selected by proximal peaks k = 1 selected gene annotated with π π π π P = Pr(k ≥1 | n=2, K =3, N=8) π π http://cs273a.stanford.edu [BejeranoFall13/14]

11. We have (reduced ChIP-Seq into) a gene list!What is the gene list enriched for? Pro: A lot of tools out there for the analysis of gene lists. Cons: These tools are built for microarray analysis. Does it matter ?? Microarray data Microarray data Generegulation data Microarray tool http://cs273a.stanford.edu [BejeranoFall13/14]

12. SRF Gene-based enrichment results • Original authors can only state: “basic cellular processes, particularly those related to gene expression” are enriched1 SRF acts on genes both in nucleus and cytoplasm, that are involved in transcription and various types of binding SRF SRF Z ~ Where’s the signal? Top “actin” term is ranked #28 in the list. ~ [1] Valouev A. et al., Nat. Methods, 2008 http://cs273a.stanford.edu [BejeranoFall13/14] 12

13. Associating only proximal peaks loses a lot of information Relationship of binding peaks to nearest genes for eight human (H) and mouse (M) ChIP-seq datasets Restricting to proximal peaks often leads to complete loss of key enrichments http://cs273a.stanford.edu [BejeranoFall13/14]

14. Bad Solution: Associating distal peaks brings in many false enrichments π π π Why bad? 14% of human genes tagged ‘multicellular organismal development’. But 33% of base pairs have such a gene nearest upstream/downstream. SRF ChIP-seq set has >2,000 binding events. Throw a random set of 2,000 regions at the genome. What do you get from a gene list analysis? Term Bonferroni corrected p-value nervous system development 5x10-9 system development 8x10-9 anatomical structure development 7x10-8 multicellular organismal development 1x10-7 developmental process 2x10-6 Large “gene deserts” are oftennext to key developmental genes http://cs273a.stanford.edu [BejeranoFall13/14]

15. Real Solution: Do not convert to gene list.Analyze the set of genomic regions Gene regulatory domain Genomic region (ChIP-seq peak) Gene transcription start site Ontology term ( ‘actin cytoskeleton’) π π π π GREAT = Genomic RegionsEnrichment of Annotations Tool π p = 0.33 of genome annotated with π n = 6 genomic regions P = Prbinom(k ≥5 | n=6, p =0.33) k = 5 genomic regions hit annotation π π π Since 33% of base pairs are near a ‘multicellular organismal development’ gene, we now expect 33% of genomic regions to hit this term by chance. => Toss 2,000 random regions at genome, get NO (false) enrichments. http://cs273a.stanford.edu [BejeranoFall13/14]

16. How does GREAT know how to assign distal binding peaks to genes? Future: High-throughput assays based on chromosome conformation capture (3C) methods will elucidate complex regulation mechanisms • Currently: Flexible computational definitions allow assignment of peaks to nearest gene, nearest two genes, etc. • Default: each gene has a “basal regulatory domain” of 5 kb up- and 1kb downstream of transcription start site, extends to basal domain of nearest genes within 1 Mb • Though some associations may be missed or incorrect, in general signal richness and robustness is greatly improved by associating distal peaks http://cs273a.stanford.edu [BejeranoFall13/14]

17. GREAT infers many specific functions of SRF from its binding profile Top GREAT enrichments of SRF Ontology Term # Genes Binomial Experimental P-value support* Top gene-based enrichments of SRF 30 31 7x10-9 5x10-5 Gene Ontology actin cytoskeleton actin binding Miano et al. 2007 Miano et al. 2007 32 26 5x10-7 2x10-6 Bertolotto et al. 2000 Poser et al. 2000 Pathway Commons TRAIL signaling Class I PI3K signaling 5 1x10-8 TreeFam Chai & Tarnawski 2002 FOS gene family 84 28 44 23 5x10-76 4x10-9 1x10-6 2x10-4 Positive control ChIp-Seq support Natesan & Gilman 1995 TF Targets (top actin-related term 28th in list) Targets of SRF Targets of GABP Targets of YY1 Targets of EGR1 * Known from literature – as in function is known, SOME of the genes are known, and the binding sites highlighted are NOT. Similar results for GABP, NRSF, Stat3, p300 ChIP-Seq [McLean et al., Nat Biotechnol., 2010] http://cs273a.stanford.edu [BejeranoFall13/14]

18. Limb P300: I was blind and I can see Gene List http://cs273a.stanford.edu [BejeranoFall13/14]

19. GREAT works with ANY cis-regulatory rich setExample: GWAS Compendium set Height-associated unlinked SNPs http://cs273a.stanford.edu [BejeranoFall13/14]

20. GREAT analysis of histone mark combinations http://cs273a.stanford.edu [BejeranoFall13/14]

21. GREAT includes multiple ontologies • Twenty ontologies spanning broad categories of biology • 44,832 total ontology terms tested in each GREAT run (2,800 terms) (6,700) (5,215) (3,079) (834) (911) (5,781) (615) (427) (19) (456) (222) (9) (150) (1,253) (6,857) (288) (8,272) (706) (238) Michael Hiller http://cs273a.stanford.edu [BejeranoFall13/14]

22. Advantages of the GREAT approach Tailored to the biology of gene regulation: • Distal sites are incorporated, not ignored • Variable length gene regulatory domains • Multiple bindings next to same target gene rewarded • Binding sites associated to (both) TSS, not gene body • Extensive ontologies, some home-made http://cs273a.stanford.edu [BejeranoFall13/14]

23. Algorithmic Optimization: A it works; B make it efficient http://cs273a.stanford.edu [BejeranoFall13/14]

24. enter GREAT.stanford.edu Choose genome Input peak list Hit submit! http://cs273a.stanford.edu [BejeranoFall13/14]

25. GREAT web app: (Optional): alter association rules http://great.stanford.edu Three association rule choices Lnp Evx2 HoxD cluster Literature-curated domains for a small subset of genes [adapted from Spitz, Gonzalez, & Duboule, Cell, 2003] http://cs273a.stanford.edu [BejeranoFall13/14]

26. GREAT web app: output summary Additional ontologies, term statistics, multiple hypothesis corrections, etc. Ontology-specific enrichments Cool visualization opportunities! http://cs273a.stanford.edu [BejeranoFall13/14]

27. GREAT web app: term details page Genes annotated as “actin binding” with associated genomic regions Genomic regions annotated with “actin binding” Drill down to explore how a particular peak regulates Plectin and its role in actin binding http://cs273a.stanford.edu [BejeranoFall13/14]

28. You can also submit any trackstraight from the UCSC Table Browser A simple, well documented programmatic interface allows any tool to submit directly to GREAT. (See our Help / Inquiries welcome!) http://cs273a.stanford.edu [BejeranoFall13/14]

29. GREAT web app: export data HTML output displays all user selected rows and columns Tab-separated values also available for additional postprocessing http://cs273a.stanford.edu [BejeranoFall13/14]

30. GREAT Web Stats 200-400 job submissions per day, from 8,000 IP addrs Over 175,000 jobs served, hundreds of citations http://cs273a.stanford.edu [BejeranoFall13/14]

31. Adding a new species to GREAT We need: A good assembly A high quality gene set Good gene annotations* *Most valuable for species with independent annotations! http://cs273a.stanford.edu [BejeranoFall13/14]

32. Test case: early neocortex development E12.5 E14.5 E16.5 Dissected the dorsal cerebral wall We Performed p300, H3K27Ac ChIP-seq. What have we learned? http://cs273a.stanford.edu [BejeranoFall13/14]

33. Collect, Call Peaks: Mostly Distal Did the experiment work? http://cs273a.stanford.edu [BejeranoFall13/14]

34. In silico validation Most enriched terms: what you want, where & when you want it http://cs273a.stanford.edu [BejeranoFall13/14]

35. Are all cell populations represented? SVZ-IZ CP VZ [Ayoubet al., PNAS. 2011] http://cs273a.stanford.edu [BejeranoFall13/14]

36. What are different enhancer groups doing? http://cs273a.stanford.edu [BejeranoFall13/14]

37. Which genes are most regulated? [Ayoub et al, 2011] http://cs273a.stanford.edu [BejeranoFall13/14]

38. The most densely regulated genes http://cs273a.stanford.edu [BejeranoFall13/14]

39. Mn1 – Cryba4 gene desert http://cs273a.stanford.edu [BejeranoFall13/14]

40. ChIP-seq is a powerful technique crosslink ChIP-seq lessons: Transcription factors (TF) bind nextto many relevant target genes. Example: SRF regulates actin cytoskeleton.Binds near 40 actin genes in immune cells. Most binding is distal. Ex: 75% of SRF binding sites near actin genes well over 1kb from TSS (upto 300kb). Binding is context specific. Ex: SRF also regulates muscle development.SRF is not enriched near muscle genes,when assayed in immune cells. Many TF functions yet undiscovered. fragment antibody immunoprecip. high throughput sequencing map to genome + find peaks http://cs273a.stanford.edu [BejeranoFall13/14]

41. How can we discover new TF functions? • Exciting new technologies: designer genome editing. • Which TF to knock out next? Where to look for its effects? • ChIP-seq requires • Quality antibody • Enough cells of right context • Sequencing costs •  Seldom used in exploratory mode • (ENCODE sampled 100 mostly basal TFs in very few contexts) • TF function exploration: expression perturbations, genetics screens. • Which TF next? Where to look? Educated guesses. • Goal: Devise a rapid method for TF function prediction http://cs273a.stanford.edu [BejeranoFall13/14]

42. How can we discover new TF functions? • Observations: • Gene function annotation is very rich.(Examples: GO, MGI expression & phenotype, etc. etc.) • Conserved non-coding, likely cis-regulatory, sequencemakes up 5-10% of the mammalian genomes. • Conserved binding site prediction has become accurate.(Predict only with strong evidence, minimize false positives) • Idea: • Use transcription factor tendency to bind next to manycontext-relevant genes to make function predictions. http://cs273a.stanford.edu [BejeranoFall13/14]

43. Transcription Factor Function Prediction 1. Curate a rich library of TF binding site motifs.Pro: hundreds are now known. (inc. from Tim, Jussi’s work) 2. Predict cross-species conserved binding sites.Pro: careful conserved binding site prediction is accurate. 3. Search for extreme binding site concentration next to genes of any particular function.Pro: Leverage the observed phenomenon, and the very large body of knowledge about all (target) genes in the genome. Con: You will not get all binding sites of any function. Pro: You may predict many diverse TF functions. Predict SRF to regulate: Actin cytoskeleton, Muscle development, ... SRF http://cs273a.stanford.edu [BejeranoFall13/14]

44. This is what we did: 1. Build motif library 300 motifs 332 motifs for 300 transcription factors (updating for many more) http://cs273a.stanford.edu [BejeranoFall13/14]

45. 2. Predict conserved binding sites . = same as human TTTCCCTTAAAAGGCTTAAATAAACTCACCAGTGTTTAATT ......................................... ......................................... ...T..................................... ...T..................................... ...T..............G...........G......T... ...T............C........AT..TG.....C.... ...T................C....AT...G.....C.... ...T.........................TG.......... ...T.....................C...TG.....G.G.. ...T.........................TG........CG ...T.........................TG.......... ...T......................T...G......G... ...T..................................... ...T........................TTG.......... ...T..........................GAC.......A ...T..........................C.......... ...T......................T.-.CA......... ...T...............G..................... ...T..................................... ...T..................................... ...T.........................TG.......... ...T.........................TG.......... TTTCCCTTAAAAGGCTTAAATAAACTCACCAGTGTTTAATT Human Chimp Gorilla Orangutan Rhesus Tarsier Mouse lemur Bushbaby Tree shrew Mouse Guinea pig Squirrel Rabbit Alpaca Cow Cat Microbat Megabat Hedgehog Rock hyrax Tenrec Armadillo Sloth Human • We in fact allow: • Imperfect motif matches • Binding site / alignment wobble • Subset of species support • Guard against alignment fragmentations • Predict efficiently • Improve state of the art using“Excess conservation” scoring http://cs273a.stanford.edu [BejeranoFall13/14]

46. Compare to shuffled motifs & weed out! SRF motif shuffle #1 shuffle #2 shuffle #3 shuffle #4 shuffle #10 … http://cs273a.stanford.edu [BejeranoFall13/14]

47. Use three reference genomesrich in gene function annotations Predicted for human, for mouse and for zebrafish http://cs273a.stanford.edu [BejeranoFall13/14]

48. 3. Predict binding site/TF functions MYOG BMP4 MYF6 NKX2-5 gene transcription start site SRF binding site CAV3 ACTA1 Enhancer to gene association  SRF must regulate muscle structure development: predicted to bind next to 157 genes (p=7.43×10-41) http://cs273a.stanford.edu [BejeranoFall13/14]

49. Conservative multiple testing correction • Control false positives when predicting function for hundreds of transcription factors: SRF motif shuffle #1 shuffle #2 shuffle #3 shuffle #4 Remove terms with E[occurrences] ≥ 1 kidney morphogenesis absent semicircular canals dorsal spinal cord development … shuffle #10 Retains only 5% of GREAT predictions. 2,543 transcription factor to function links (16% FDR) http://cs273a.stanford.edu [BejeranoFall13/14]

50. PRISM vs. ChIP-seq Actin cytoskeleton SRF T cell TermPRISMChIP-seq actin cytoskeleton Known structural constituent of muscle Known dilated heart ventricles Known regulation of insulin secretion Novel muscle SRF heart http://cs273a.stanford.edu [BejeranoFall13/14]