Case Study One: attempt to reproduce known regulation of NFE2 by GATA1 in bulk RNA-seq

Introduction

The recently published Benchmark and integration of resources for the estimation of human transcription factor activities provides a useful curation of TF/target gene relations divided into five confidence categories. We chose one high confidence (score “A”) relation from this dataset for this case study: the regulation of NFE2 by GATA1 in erythropoiesis. We will use trena to “discover” this relationship, a necessary but not sufficient proof of trena’s value.

Hematopoiesis (the production of blood cells) involves two subprocesses: erythropoiesis (red blood cell biogenesis) and megakaryopoiesis (megakaryocytes and platelets). These processes have been intensively studied for decades, and much has been learned about the complex and hierarchical regulation of these processes of differentiation. GATA1, TALl, KLF1 and the p45 isoform of NFE2 are some of the transcription factors involved.

The Benchmark dataset reports two studies[ref1, ref2] reporting transcriptional regulation of NFE2 by GATA1. In addition, Takayama et al, 2010 report that “In megakaryocytes, GATA1 deficiency reduces p45 expression by approximately 50%, indicating the presence of GATA1-dependent and -independent regulation of the p45 gene.” We lack precise knowledge of the intricacies of regulatory networks in erythropoieis, but the cited research in combination strongly suggests the regulation by GATA1 of NFE2 in erythropoieis.

GTEx Blood

The GTEx project provides RNA-seq count expression data for 25k genes and 755 samples. These data are highly heterogeneous, derived from a mix of red blood cells, white blood cells and platelets from many subjects. It should be no surprise therefore that our trena analysis fails to retrieve a regulatory relationship between GATA1 and NFE2 from this dataset. This exercise introduces the trena technique and the dangers of bulk heterogeneous data.

Expression and all known transcription factors

The GeneOntology project annotates 1663 human genes to the molecular function DNA-binding transcription factor activity:

> tfs.all <- sort(unique(select(org.Hs.eg.db, keys="GO:0003700", keytype="GOALL",  columns="SYMBOL")$SYMBOL))
> tfs <- intersect(tfs.all, rownames(mtx.gtex.blood)
> solvers <-  c("lasso", "lassopv", "pearson", "randomForest", "ridge", "spearman", "xgboost")
> trenaEnsemble <- EnsembleSolver(mtx.gtex.blood, "NFE2", tfs, solvers)
> tbl.model <- run(trenaEnsemble)
> head(tbl.model, n=20)
      gene  betaLasso   lassoPValue pearsonCoeff     rfScore   betaRidge spearmanCoeff      xgboost
1   ZNF467 0.04554979 8.102352e-270    0.8974169 179.4721705 0.017900579     0.9059445 5.378471e-01
2     GAS7 0.04214837 5.185498e-254    0.8928340 138.1390796 0.009325281     0.8954556 1.856330e-01
3   ZNF438 0.12869763 1.108818e-216    0.8868706  81.0674620 0.017991204     0.8905471 6.816107e-02
4     SPI1 0.18316329 1.786604e-114    0.8748366  44.8748867 0.016571769     0.8931940 2.611708e-02
5     MXD3 0.15625475 9.154617e-181    0.8743369  46.8945593 0.010499558     0.8840492 8.380720e-03
6    MLLT1 0.04897373  8.023350e-44    0.8719094  36.9168331 0.007973757     0.8826423 8.462997e-03
7     RFX2 0.03938209  2.781065e-39    0.8611165  32.4493142 0.010776699     0.8703952 1.977503e-03
8     ATF6 0.00000000  1.315896e-01    0.8198707   0.3438544 0.014297394     0.8249904 2.738534e-05
9    CREB5 0.00000000  4.501187e-02    0.8198528   9.6779565 0.016295689     0.8085946 2.863106e-02
10   CEBPD 0.00000000  1.000000e+00    0.8142818   0.1316654 0.008440315     0.8224173 3.893420e-05
11   CEBPA 0.00000000  2.317224e-02    0.8130808   8.6251990 0.016003815     0.8356710 3.202592e-04
12 SUPT4H1 0.00000000  1.000000e+00    0.8105414   3.4865286 0.006986417     0.8311199 1.528261e-03
13    HHEX 0.00000000  1.000000e+00    0.8087241   1.0232650 0.003709475     0.8092507 2.668327e-07
14  ZBTB7B 0.00000000  1.000000e+00    0.8079335   7.7857724 0.011015725     0.8368950 7.776020e-04
15   STK16 0.00000000  7.855677e-01    0.8056880   1.0425907 0.006133720     0.8150958 1.706304e-03
16    TFEB 0.00000000  1.000000e+00    0.7994628   0.2934284 0.004787306     0.8206329 3.367513e-05
17  ZNF787 0.00000000  1.000000e+00    0.7940764   1.2943307 0.011087230     0.8239310 4.486712e-04
18     HLX 0.06593066  3.088806e-41    0.7876786   4.2640642 0.020916575     0.8083032 6.963889e-03
19  ZNF746 0.00000000  1.000000e+00    0.7875233   0.4593815 0.011661213     0.7901268 2.914801e-05
20   MLXIP 0.00000000  1.000000e+00    0.7871229   0.4650622 0.002449240     0.7988580 8.521462e-05

We hope to see the three erythropoiesis transcription factors somewhere in the trena model. In this one, they are not prominent, found at positions 371 725 and 824 in the model as sorted by the absolute value of the pearson coefficient:

> match(c("GATA1", "TAL1", "KLF1"), tbl.lps$gene)   # 
[1] 371 725 824

Expression and only transcription factors with known motifs

The JASPAR 2018 and Hocomoco transcription factor compendia, when combine, identify 780 annotated transcription factor motif. In building the next model, candidate transcription factors are limited to this set.

> tfs.withMotifs <- mcols(query(MotifDb, c("sapiens"), c("jaspar2018", "hocomoo")))$geneSymbol
> tfs <- intersect(tfs.withMotifs, rownames(mtx.gtex.blood)
> solvers <-  c("lasso", "lassopv", "pearson", "randomForest", "ridge", "spearman", "xgboost")
> trenaEnsemble <- EnsembleSolver(mtx.gtex.blood, "NFE2", tfs, solvers)
> tbl.model <- run(trenaEnsemble)
> head(tbl.model, n=20)
     gene  betaLasso   lassoPValue pearsonCoeff     rfScore   betaRidge spearmanCoeff      xgboost
1    SPI1 0.39179921 3.105960e-239    0.8748366 173.1813469 0.038256090     0.8931940 4.484164e-01
2    RFX2 0.13283782 1.236616e-203    0.8611165 113.7966516 0.028434701     0.8703952 2.578767e-01
3   CREB5 0.01020744  4.215279e-27    0.8198528  29.3836851 0.022581439     0.8085946 2.582643e-02
4   CEBPD 0.00000000  7.806114e-01    0.8142818   9.8128799 0.012064964     0.8224173 2.684192e-05
5   CEBPA 0.04856000 4.918991e-124    0.8130808  24.5906707 0.020001647     0.8356710 9.968208e-04
6  ZBTB7B 0.00000000  7.570965e-01    0.8079335  63.6673040 0.019478331     0.8368950 2.942618e-02
7    TFEB 0.00000000  9.128677e-01    0.7994628  11.5517782 0.016902326     0.8206329 3.371606e-05
8    BATF 0.05800616  1.182795e-77    0.7799541   5.0297895 0.022430686     0.7874306 3.195642e-03
9    RARA 0.00000000  1.432931e-04    0.7768149  15.1423384 0.024408500     0.8108477 3.022155e-03
10    SP1 0.00000000  2.669736e-01    0.7766227   1.1118820 0.022444381     0.7866711 1.295055e-04
11  KLF14 0.00000000  9.461310e-01    0.7733890  59.8682708 0.014609629     0.8102489 3.570402e-03
12  NR2E1 0.07466256 1.645394e-128    0.7674846  73.1777412 0.031161125     0.8166762 5.374055e-02
13   IRF2 0.00000000  1.147759e-01    0.7457406   0.2705582 0.022106204     0.7481861 2.877418e-04
14  MEF2A 0.00000000  7.640160e-01    0.7310054   0.4375018 0.016944211     0.7305359 1.044375e-05
15  GLIS2 0.00000000  8.645402e-01    0.7198283   0.9845279 0.007047398     0.7382333 3.208840e-05
16   JDP2 0.00000000  6.085086e-01    0.7195525   0.1852113 0.002444830     0.7269272 1.237316e-04
17   TFE3 0.00000000  3.867438e-08    0.7144317   1.7343149 0.026801329     0.7308633 7.625371e-04
18   BCL6 0.00000000  2.100360e-08    0.7088486   4.5882372 0.020731165     0.6988929 1.135068e-03
19    SRF 0.00000000  8.907811e-01    0.7074436   0.5702123 0.013987309     0.7270949 4.228163e-04
20 PKNOX1 0.00000000  7.175965e-01    0.7052011   0.2088824 0.013225756     0.7120145 8.979636e-05

The three TFs score higher in this model, but can not be said, in any strong sense, to be predicted by it as regulators of NFE2.

> match(c("GATA1", "TAL1", "KLF1"), tbl.model$gene) 
[1]  125 225 252

Expression and highly-conserved, high-scoring transcription factors in a 20kb regulatory region

We hypothesize that transcription factors with well-matched motifs found in highly conserved regulatory regions within +/- 10kb of the target gene’s TSS are more likely than random to be functional binding sites. When found, and when tf/target gene expression is also correlated, these may be viewed as possibly sound trena predictions.

Here we use a precalculated table of FIMO and phast7 scores for 20kb surrounding the NFE2 transcription start site, extracting only those TFs with very high match and conservation. With these data and assumptions, GATA1 rises to rank 15 in the model, but appears as a repressor - contrary to expectation and the findings of the published papers.

> tbl.tfs.elite <- subset(tbl.fimoMotifs, p.value <= fimo.score & phast7 >= phast.score)
> dim(tbl.tfs.elite)
> tfs <- sort(unique(tbl.tbs.elite)$tf)
> length(tfs)  # 52
> match(c("GATA1", "TAL1", "KLF1"), tfs.elite)   # 13 41 16
> solver <- EnsembleSolver(mtx.blood.lps, target.gene, tfs, geneCutoff=1.0, solverNames=solverNames)
> tbl <- run(solver)
> new.order <- order(abs(tbl$pearsonCoeff), decreasing=TRUE)
> tbl <- tbl[new.order,]
> rownames(tbl) <- NULL
> tbl.fimo.phast.stringent <- tbl
> head(tbl.fimo.phast.stringent, n=20)
> head(tbl.fimo.phast.stringent, n=20)
    gene   betaLasso pearsonCoeff    rfScore    betaRidge spearmanCoeff      xgboost
1   SPI1  0.54116978    0.8748366 232.581245  0.174503316     0.8931940 0.7557817023
2  CEBPA  0.23523082    0.8130808 139.720676  0.127110447     0.8356710 0.0281807591
3   RARA  0.00000000    0.7768149  78.055568  0.110157856     0.8108477 0.0551366432
4    SP1  0.06980607    0.7766227  72.914116  0.114400926     0.7866711 0.0026924588
5  KLF16  0.00000000    0.6804686  33.001648  0.065192653     0.7286630 0.0040067050
6    MNT  0.00000000    0.6803979  12.478543  0.056498612     0.6911230 0.0007373286
7  NR6A1  0.00000000    0.6480948   8.466320  0.038737218     0.6580332 0.0013781294
8  THAP1  0.01304063    0.6299785  10.929519  0.071647737     0.6430135 0.0012611141
9  STAT3  0.00000000    0.5609655   4.885006  0.078064628     0.5551787 0.0017422688
10  ELF2  0.00000000    0.5390169   4.306037  0.037256849     0.5308688 0.0007064301
11 TFCP2  0.00000000    0.5278668   2.716275  0.045926049     0.5175253 0.0005061365
12  EGR1 -0.12759104   -0.5131861  32.791422 -0.092334479    -0.5532816 0.0590037734
13  FLI1  0.00000000    0.4952012   1.395661  0.013529524     0.4905536 0.0029439625
14  NFIC  0.00000000    0.4724385   3.651387  0.032486039     0.4540451 0.0007173075
15 GATA1  0.00000000   -0.3971076   4.030233 -0.015668836    -0.4394136 0.0016572605
16   MAZ  0.00000000    0.3159501   2.228739  0.021758139     0.2773423 0.0004919320
17  KLF9  0.00000000    0.3029702   1.480844 -0.004387528     0.3096802 0.0011933040
18  IRF4 -0.01352598   -0.2827062   7.656989 -0.070280645    -0.3225595 0.0008688548
19 PROX1  0.00000000   -0.2742029   1.997726 -0.012063810    -0.3004936 0.0009774531
20   SP3  0.00000000    0.2691114   1.473230 -0.010016557     0.2725386 0.0003386794

This model is built with only the 52 TFs which pass the fimo (sequence match) and phast7 (sequence conservation) thresholds. GATA1, TAL1 and KLF1 are in this much smaller group, so necessarily have a higher rank in the resulting model. Note, however, the negative correlation with the target gene, NFE2.

> match(c("GATA1", "TAL1", "KLF1"), tbl.fimo.phast.stringent$gene)
[1] 15 27 30

Expression and somewhat-conserved, FIMO default matched transcription factors in a 20kb regulatory region

Here we see that loosening the FIMO match threshold to its traditional default value (1e-4) and phast7 conservation to 20% does not increase our ability to predict the regulation of NFE2 by GATA1:

> tbl.tfs.weak <- subset(tbl.fimoMotifs, p.value <= 1e-4 & phast7 > 0.2)
> nrow(tbl.tfs.weak)   # 2846
> tfs.weak <- unique(tbl.tfs.weak$tf)
> length(tfs.weak)  # 525
> match(c("GATA1", "TAL1", "KLF1"), tfs.weak)   # 281 282 330
> 
> solver <- EnsembleSolver(mtx.blood.lps, target.gene, tfs.weak, geneCutoff=1.0, solverNames=solverNames)
> tbl <- run(solver)
> dim(tbl)  # 378
> new.order <- order(abs(tbl$pearsonCoeff), decreasing=TRUE)
> tbl <- tbl[new.order,]
> rownames(tbl) <- NULL
> tbl.fimo.phast.weak <- tbl
> head(tbl.fimo.phast.weak, n=20)
     gene   betaLasso pearsonCoeff     rfScore    betaRidge spearmanCoeff      xgboost
1    SPI1  0.37925044    0.8748366 173.3335617  0.048503015     0.8931940 4.866324e-01
2    RFX2  0.13351525    0.8611165 148.3728199  0.033019631     0.8703952 2.471701e-01
3   CEBPD  0.00000000    0.8142818  14.1532996  0.017885322     0.8224173 6.461158e-05
4   CEBPA  0.06400950    0.8130808  48.0292068  0.030222210     0.8356710 3.197183e-03
5    BATF  0.06772117    0.7799541   8.7910003  0.032610855     0.7874306 3.573105e-03
6    RARA  0.00000000    0.7768149  19.7713182  0.029306674     0.8108477 3.475518e-03
7     SP1  0.00000000    0.7766227   6.5621383  0.026881993     0.7866711 1.750126e-04
8   KLF14  0.00000000    0.7733890  57.5545230  0.015930889     0.8102489 9.025739e-04
9   NR2E1  0.08100175    0.7674846  79.9225143  0.037686012     0.8166762 6.480617e-02
10   IRF2  0.00000000    0.7457406   0.5893363  0.025511229     0.7481861 1.858361e-04
11  MEF2A  0.00000000    0.7310054   0.5754763  0.019650713     0.7305359 1.326880e-04
12  GLIS2  0.00000000    0.7198283   1.6554741  0.006859244     0.7382333 2.377722e-05
13   JDP2  0.00000000    0.7195525   0.3562898  0.001063041     0.7269272 3.631474e-04
14   BCL6  0.00000000    0.7088486   8.9160346  0.029732369     0.6988929 2.176104e-03
15    SRF  0.00000000    0.7074436   0.9801356  0.019549681     0.7270949 2.524539e-04
16 PKNOX1  0.00000000    0.7052011   0.2873513  0.013340846     0.7120145 1.490372e-05
17   RXRA  0.00000000    0.7036052   1.3152228  0.020967989     0.7148850 4.554199e-04
18    JUN -0.08756301   -0.7030002  26.9067254 -0.015677627    -0.6639774 7.007604e-02
19  NR4A1 -0.01683527   -0.7025484   8.3418147 -0.017519536    -0.7020151 1.222394e-03
20 ZNF282  0.00000000    0.6916040   2.5692101  0.025967945     0.7208917 2.445573e-03

With this more inclusive list of TFs, less stringent with respect to sequence match and conservation, our three genes of interest fall to lower positions in the model. For this reason, and for the anti-correlation of GATA1, we decline to call this result a prediction of NFE2 regulation.

> match(c("GATA1", "TAL1", "KLF1"), tbl.fimo.phast.weak$gene)
> [1]  95 168 190

Conclusion

None of the four strategies used above recapitulate the known regulation of NFE2 by GATA1 in erythropoiesis. This probably reflects two factors:

  1. heterogeneous bulk data carries little regulatory signal
  2. erythropoiesis (is this true, Cory et al?) is essentially complete by the time blood is circulating, have reached terminal differentiation states?

References