1 Summary

1.1 ELF3 family and subfamily

load("~/paper_elmer/ESCA_analysis_hg38/escc_vs_esad_hg38/hyper/ELMER_results_hyper.rda")
family <- ELMER:::getdata("TF.family")
family$ELF3_HUMAN.H11MO.0.A %>% sort

##  [1] "EHF"   "ELF1"  "ELF2"  "ELF3"  "ELF4"  "ELF5"  "ELK1"  "ELK3" 
##  [9] "ELK4"  "ERF"   "ERG"   "ETS1"  "ETS2"  "ETV1"  "ETV2"  "ETV3" 
## [17] "ETV3L" "ETV4"  "ETV5"  "ETV6"  "ETV7"  "FEV"   "FLI1"  "GABPA"
## [25] "SPDEF" "SPI1"  "SPIB"  "SPIC"

subfamily <- ELMER:::getdata("TF.subfamily")
subfamily$ELF3_HUMAN.H11MO.0.A

## [1] "EHF"  "ELF3" "ELF5"

1.2 Enriched motifs and probes associated

# which is the enriched motif in the family
en.motif.family <- unlist(sapply(family$ELF3_HUMAN.H11MO.0.A,
                                 function(x) TF$motif[grep(x,TF$top.potential.TF.family)])
)
en.motif.family

##                 SPDEF1                 SPDEF2                 SPDEF3 
## "ELF3_HUMAN.H11MO.0.A"  "EHF_HUMAN.H11MO.0.B" "ELF5_HUMAN.H11MO.0.A" 
##                 SPDEF4                 SPDEF5                 SPDEF6 
## "FLI1_HUMAN.H11MO.1.A" "ETV5_HUMAN.H11MO.0.C" "ETV2_HUMAN.H11MO.0.B" 
##                 SPDEF7                 SPDEF8                 SPDEF9 
## "ETV4_HUMAN.H11MO.0.B" "ETS1_HUMAN.H11MO.0.A" "ELK1_HUMAN.H11MO.0.B" 
##                SPDEF10                SPDEF11                SPDEF12 
##  "ERG_HUMAN.H11MO.0.A" "ETS2_HUMAN.H11MO.0.B" "SPI1_HUMAN.H11MO.0.A" 
##                SPDEF13                SPDEF14                SPDEF15 
## "ELF1_HUMAN.H11MO.0.A" "ELK4_HUMAN.H11MO.0.A"  "FEV_HUMAN.H11MO.0.B" 
##                SPDEF16 
## "ELF2_HUMAN.H11MO.0.C"

en.motif.subfamily <- unlist(sapply(subfamily$ELF3_HUMAN.H11MO.0.A,
                                    function(x) TF$motif[grep(x,TF$top.potential.TF.subfamily)])
)
en.motif.subfamily

##                  ELF31                  ELF32                  ELF33 
## "ELF3_HUMAN.H11MO.0.A"  "EHF_HUMAN.H11MO.0.B" "ELF5_HUMAN.H11MO.0.A"

probes.en.family <- Reduce(union,enriched.motif[en.motif.family])
length(probes.en.family)

## [1] 663

length(Reduce(intersect,enriched.motif[en.motif.family]))

## [1] 16

probes.en.subfamily <- Reduce(union,enriched.motif[en.motif.subfamily])
length(probes.en.subfamily)

## [1] 313

length(Reduce(intersect,enriched.motif[en.motif.subfamily]))

## [1] 141

2 Narrow ELF3 peaks

colnames(peaks) <- c("chrom","Start","End","Name","Width","Score", "fold_enrichment","-log10(pValue)","-log10(qValue)","peak")
peaks$chrom <- paste0("chr",peaks$chrom)
peaks.gr <- makeGRangesFromDataFrame(peaks,keep.extra.columns = T)
DT::datatable(peaks,options = list(scrollX = TRUE),filter = 'top')

ChIPseeker::covplot(peaks.gr, weightCol="fold_enrichment")

3 Comparing paired probes considering family classification

3.1 Get ELF3 probes target regions

pairs <- read_csv("~/paper_elmer/ESCA_analysis_hg38/escc_vs_esad_hg38/hyper/getTFtargets_genomic_coordinates_mapped_to_hg19.hyper.family.csv",
                  col_types = readr::cols())
targets <- pairs[grep("ELF3",pairs$TF),]
DT::datatable(targets,options = list(scrollX = TRUE),filter = 'top')

probes <- makeGRangesFromDataFrame(targets,
                                   keep.extra.columns = T,
                                   seqnames.field = "probe_hg19_seqnames",
                                   start.field = "probe_hg19_start",
                                   end.field =  "probe_hg19_end",
                                   strand.field = "probe_hg19_strand"
)
probes <- unique(probes)
probes

## GRanges object with 663 ranges and 19 metadata columns:
##         seqnames                 ranges strand |          GeneID
##            <Rle>              <IRanges>  <Rle> |     <character>
##     [1]     chr6   [45979349, 45979350]      * | ENSG00000001561
##     [2]     chr6   [45441433, 45441434]      * | ENSG00000001561
##     [3]     chr4   [10447440, 10447441]      * | ENSG00000002587
##     [4]     chr4   [11652505, 11652506]      * | ENSG00000002587
##     [5]    chr17   [37034884, 37034885]      * | ENSG00000002834
##     ...      ...                    ...    ... .             ...
##   [659]    chr17 [ 36070600,  36070601]      * | ENSG00000274620
##   [660]    chr10 [127897162, 127897163]      * | ENSG00000277371
##   [661]    chr17 [ 35242555,  35242556]      * | ENSG00000277589
##   [662]    chr10 [126211892, 126211893]      * | ENSG00000278831
##   [663]    chr10 [126211704, 126211705]      * | ENSG00000278831
##               Probe      Symbol       Sides        Raw.p          FDR
##         <character> <character> <character>    <numeric>    <numeric>
##     [1]  cg16474725       ENPP4          R6 1.270122e-21 1.939775e-20
##     [2]  cg04757168       ENPP4          R5 1.270122e-21 1.939775e-20
##     [3]  cg11578055      HS3ST1          R7 1.237582e-10 3.681389e-10
##     [4]  cg23971987      HS3ST1          L3 1.237582e-10 3.681389e-10
##     [5]  cg27061889       LASP1          R1 7.105825e-20 7.827895e-19
##     ...         ...         ...         ...          ...          ...
##   [659]  cg04136369     MIR378J          L2 7.562441e-05 9.524983e-05
##   [660]  cg17287034 Metazoa_SRP          L7 1.373459e-06 2.310831e-06
##   [661]  cg08659204  AC244093.4          R6 1.506934e-08 3.416747e-08
##   [662]  cg16426764  AL513190.1          R3 9.809833e-06 1.445882e-05
##   [663]  cg03499943  AL513190.1          R3 9.809833e-06 1.445882e-05
##         distNearestTSS DMC_analysis_pvalue DMC_analysis_escc_Minus_esad
##              <integer>           <numeric>                    <numeric>
##     [1]         118379        2.576139e-23                    0.3499644
##     [2]         656295        1.540132e-27                    0.3424920
##     [3]         982780        4.831871e-19                    0.3496876
##     [4]         221115        1.074089e-22                    0.3277020
##     [5]           4741        1.656185e-36                    0.3364163
##     ...            ...                 ...                          ...
##   [659]          95552        4.455193e-24                    0.3275653
##   [660]         339298        1.362491e-20                    0.3099911
##   [661]         113298        1.268293e-13                    0.3257977
##   [662]         100028        6.455562e-32                    0.3729131
##   [663]         100216        9.619440e-24                    0.3195927
##         DMC_analysis_adjust.p external_gene_name original_ensembl_gene_id
##                     <numeric>        <character>              <character>
##     [1]          3.119777e-21              ENPP4        ENSG00000001561.6
##     [2]          3.291965e-25              ENPP4        ENSG00000001561.6
##     [3]          3.257832e-17             HS3ST1        ENSG00000002587.8
##     [4]          1.202478e-20             HS3ST1        ENSG00000002587.8
##     [5]          1.227075e-33              LASP1       ENSG00000002834.16
##     ...                   ...                ...                      ...
##   [659]          5.983252e-22            MIR378J        ENSG00000274620.1
##   [660]          1.145422e-18        Metazoa_SRP        ENSG00000277371.1
##   [661]          4.014973e-12         AC244093.4        ENSG00000277589.1
##   [662]          2.425544e-29         AL513190.1        ENSG00000278831.1
##   [663]          1.242528e-21         AL513190.1        ENSG00000278831.1
##                                                                                                                                                                             TF
##                                                                                                                                                                    <character>
##     [1]                                                 HNF4A;HNF4G;NR2F2;NR2C1;FOXA3;FOXA2;FOXP4;FOXD2;FOXJ1;SPDEF;ELF3;NR1I2;PPARG;RORC;THRA;NR1H3;NR1H4;ZSCAN16;ZNF816;IRF8
##     [2] FOXA3;FOXA2;FOXP4;FOXD2;FOXJ1;TCF7L2;NR1I2;PPARG;RORC;THRA;NR1H3;NR1H4;ZNF620;ZSCAN16;ZNF774;ZNF468;ZNF33A;ZNF816;ZNF619;ZNF765;ZNF799;ZNF763;ZNF823;ZNF124;SPDEF;ELF3
##     [3]                                GATA6;GATA4;NR1I2;PPARG;RORC;THRA;NR1H3;NR1H4;TCF7L2;HNF1A;HNF1B;PDX1;CDX2;HOXB6;HOXB5;MNX1;EVX1;CDX1;HOXA13;HOXB9;SPDEF;ELF3;IRF8;SOX9
##     [4]       GATA6;GATA4;NR1I2;PPARG;RORC;THRA;NR1H3;NR1H4;TCF7L2;SPDEF;ELF3;HNF1A;HNF1B;PDX1;CDX2;HOXB6;HOXB5;MNX1;EVX1;CDX1;HOXA13;HOXB9;FOXA3;FOXA2;FOXP4;FOXD2;FOXJ1;SOX9
##     [5]                                                                                  HNF4A;HNF4G;NR2F2;NR2C1;FOXA3;FOXA2;FOXP4;FOXD2;FOXJ1;SPDEF;ELF3;NR3C2;ZSCAN16;ZNF816
##     ...                                                                                                                                                                    ...
##   [659]                                                         GATA6;GATA4;FOXA3;FOXA2;FOXP4;FOXD2;FOXJ1;SPDEF;ELF3;MECOM;PRDM16;ZNF587;ZNF792;ZNF518A;ZNF587B;ZSCAN16;ZNF816
##   [660]              FOXA3;FOXA2;FOXP4;FOXD2;FOXJ1;HNF4A;HNF4G;NR2F2;NR2C1;SPDEF;ELF3;GATA6;GATA4;HNF1A;HNF1B;PDX1;CDX2;HOXB6;HOXB5;MNX1;EVX1;CDX1;HOXA13;HOXB9;ZSCAN16;ZNF816
##   [661]           GATA6;GATA4;FOXA3;FOXA2;FOXP4;FOXD2;FOXJ1;ATOH1;NEUROG3;SPDEF;ELF3;ZNF620;ZSCAN16;ZNF774;ZNF468;ZNF33A;ZNF816;ZNF619;ZNF765;ZNF799;ZNF763;ZNF823;ZNF124;SOX9
##   [662]                                                          HNF4A;HNF4G;NR2F2;NR2C1;FOXA3;FOXA2;FOXP4;FOXD2;FOXJ1;NR1I2;PPARG;RORC;THRA;NR1H3;NR1H4;SPDEF;ELF3;SOX9;NR3C2
##   [663]                                                                                                                    FOXA3;FOXA2;FOXP4;FOXD2;FOXJ1;SPDEF;ELF3;NR3C2;SOX9
##         probe_hg19_width gene_hg19_seqnames gene_hg19_start gene_hg19_end
##                <integer>        <character>       <integer>     <integer>
##     [1]                2               chr6        46097730      46114436
##     [2]                2               chr6        46097730      46114436
##     [3]                2               chr4        11394774      11431389
##     [4]                2               chr4        11394774      11431389
##     [5]                2              chr17        37026112      37078023
##     ...              ...                ...             ...           ...
##   [659]                2              chr17        35974976      35975084
##   [660]                2              chr10       127557864     127558145
##   [661]                2              chr17        35355897      35357333
##   [662]                2              chr10       126311922     126312648
##   [663]                2              chr10       126311922     126312648
##         gene_hg19_width gene_hg19_strand
##               <integer>      <character>
##     [1]           16707                +
##     [2]           16707                +
##     [3]           36616                -
##     [4]           36616                -
##     [5]           51912                +
##     ...             ...              ...
##   [659]             109                -
##   [660]             282                +
##   [661]            1437                +
##   [662]             727                +
##   [663]             727                +
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths

3.2 Distance between probes and peaks

y <- distanceToNearest(probes,peaks.gr,ignore.strand=T)
qplot(mcols(y)$distance,geom="histogram") + ggthemes::theme_gdocs()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

qplot(log(mcols(y)$distance + 1),geom="histogram") + ggthemes::theme_gdocs()

## `stat_bin()` using `bins = 30`. Pick better value with `binwidth`.

3.3 Results overlap

3.3.1 Overlap between probes (+-250bp) and peaks

probes_extended_500bp <- resize(unique(probes),width = 501,fix = "center")
probes_extended_500bp[,1:2]

## GRanges object with 663 ranges and 2 metadata columns:
##         seqnames                 ranges strand |          GeneID
##            <Rle>              <IRanges>  <Rle> |     <character>
##     [1]     chr6   [45979099, 45979599]      * | ENSG00000001561
##     [2]     chr6   [45441183, 45441683]      * | ENSG00000001561
##     [3]     chr4   [10447190, 10447690]      * | ENSG00000002587
##     [4]     chr4   [11652255, 11652755]      * | ENSG00000002587
##     [5]    chr17   [37034634, 37035134]      * | ENSG00000002834
##     ...      ...                    ...    ... .             ...
##   [659]    chr17 [ 36070350,  36070850]      * | ENSG00000274620
##   [660]    chr10 [127896912, 127897412]      * | ENSG00000277371
##   [661]    chr17 [ 35242305,  35242805]      * | ENSG00000277589
##   [662]    chr10 [126211642, 126212142]      * | ENSG00000278831
##   [663]    chr10 [126211454, 126211954]      * | ENSG00000278831
##               Probe
##         <character>
##     [1]  cg16474725
##     [2]  cg04757168
##     [3]  cg11578055
##     [4]  cg23971987
##     [5]  cg27061889
##     ...         ...
##   [659]  cg04136369
##   [660]  cg17287034
##   [661]  cg08659204
##   [662]  cg16426764
##   [663]  cg03499943
##   -------
##   seqinfo: 23 sequences from an unspecified genome; no seqlengths

unique(width(probes_extended_500bp[,1:2]))

## [1] 501

hits <- suppressWarnings(findOverlaps(peaks.gr,probes_extended_500bp))
probes.hit <- subjectHits(hits)
length(unique(probes_extended_500bp[probes.hit])) # unique regions

## [1] 63

3.3.2 Overlap between peaks and distall probes (+-250bp)

distal.probes <- ELMER::get.feature.probe(feature = NULL,
                                          genome = "hg19", 
                                          met.platform = "450K")

## Downloading transcripts information. Using: Homo sapiens genes (GRCh37.p13)

## Returning distal probes: 163576

distal.probes.extended_500bp <- resize(unique(distal.probes),width = 501,fix = "center")

hits.distal <- suppressWarnings(findOverlaps(peaks.gr,distal.probes.extended_500bp))
distal.probes.hit <- subjectHits(hits.distal)

3.3.3 Regions overlapped

DT::datatable(tibble::as_data_frame(probes[subjectHits(hits)]),
              options = list(scrollX = TRUE),filter = 'top')

DT::datatable(tibble::as_data_frame(peaks.gr[queryHits(hits)]),
              options = list(scrollX = TRUE),filter = 'top')

3.3.4 Enrichement between peaks and all distal probes

x <- matrix(c(
  length(unique(probes[probes.hit])), # all linked probes with hits
  (length(unique(probes)) - length(unique(probes[probes.hit]))), # all linked probes with no hits
  length(unique(distal.probes.hit)), # all distal probes with  hits (4216)
  (length(distal.probes)) - length(unique(distal.probes.hit))), # all distal probes with no hits
  nrow = 2,
  dimnames = list(c("Overlap", "No overlap"),
                  c("Unique Linked probes (family)", "Distal probes"))
)
x <- t(x)
x

##                               Overlap No overlap
## Unique Linked probes (family)      63        600
## Distal probes                    1778     161798

fisher.test(x)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##   7.219744 12.460081
## sample estimates:
## odds ratio 
##   9.554989

3.3.5 Enrichement between peaks and distal probes with signature for family members

motifs.probes <- ELMER:::getdata("Probes.motif.hg38.450K")

probes.family <- rownames(motifs.probes)[rowSums(as.matrix(motifs.probes[,en.motif.family])) > 0]
distal.probes.family.motif <- distal.probes.extended_500bp[names(distal.probes) %in% probes.family]
hits.distal.family.motif <- suppressWarnings(findOverlaps(peaks.gr,distal.probes.family.motif))
distal.probes.family.motif.hit <- subjectHits(hits.distal.family.motif)

x <- matrix(c(
  length(unique(probes[probes.hit])), # all linked probes with hits
  (length(unique(probes)) - length(unique(probes[probes.hit]))), # all linked probes with no hits
  length(unique(distal.probes.family.motif.hit)), # all distal probes with  hits (4216)
  (length(distal.probes.family.motif)) - length(unique(distal.probes.family.motif.hit))), # all distal probes with no hits
  nrow = 2,
  dimnames = list(c("Overlap", "No overlap"),
                  c("Unique Linked probes (family)", "Distal probes with family motif signature"))
)
x <- t(x)
x

##                                           Overlap No overlap
## Unique Linked probes (family)                  63        600
## Distal probes with family motif signature    1613     113000

fisher.test(x)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  5.555149 9.596420
## sample estimates:
## odds ratio 
##   7.356554

3.3.6 Enrichement between peaks and distal probes without signature for family members

# Get probes family/subfamily
motifs.probes <- ELMER:::getdata("Probes.motif.hg38.450K")
probes.family <- rownames(motifs.probes)[rowSums(as.matrix(motifs.probes[,en.motif.family])) > 0]
probes.subfamily <- rownames(motifs.probes)[rowSums(as.matrix(motifs.probes[,en.motif.subfamily])) > 0]


distal.probes.family.nomotif <- distal.probes.extended_500bp[!names(distal.probes) %in% probes.family]
hits.distal.family.nomotif <- suppressWarnings(findOverlaps(peaks.gr,distal.probes.family.nomotif))
distal.probes.family.nomotif.hit <- subjectHits(hits.distal.family.nomotif)

x <- matrix(c(
  length(unique(probes[probes.hit])), # all linked probes with hits
  (length(unique(probes)) - length(unique(probes[probes.hit]))), # all linked probes with no hits
  length(unique(distal.probes.family.nomotif.hit)), # all distal probes with  hits (4216)
  (length(distal.probes.family.nomotif)) - length(unique(distal.probes.family.nomotif.hit))), # all distal probes with no hits
  nrow = 2,
  dimnames = list(c("Overlap", "No overlap"),
                  c("Unique Linked probes (family)", "Distal probes without family motif signature"))
)
x <- t(x)
x

##                                              Overlap No overlap
## Unique Linked probes (family)                     63        600
## Distal probes without family motif signature     165      48798

fisher.test(x)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  22.58514 42.27028
## sample estimates:
## odds ratio 
##   31.03181

4 Comparing paired probes considering subfamily classification

pairs <- read_csv("~/paper_elmer/ESCA_analysis_hg38/escc_vs_esad_hg38/hyper/getTFtargets_genomic_coordinates_mapped_to_hg19.hyper.subfamily.csv",
                  col_types = readr::cols())
targets <- pairs[grep("ELF3",pairs$TF),]
probes <- makeGRangesFromDataFrame(targets,
                                   keep.extra.columns = T,
                                   seqnames.field = "probe_hg19_seqnames",
                                   start.field = "probe_hg19_start",
                                   end.field =  "probe_hg19_end",
                                   strand.field = "probe_hg19_strand"
)
probes <- unique(probes)
probes

## GRanges object with 313 ranges and 19 metadata columns:
##         seqnames                 ranges strand |          GeneID
##            <Rle>              <IRanges>  <Rle> |     <character>
##     [1]     chr6   [45979349, 45979350]      * | ENSG00000001561
##     [2]    chr17   [37034884, 37034885]      * | ENSG00000002834
##     [3]     chrX   [11773466, 11773467]      * | ENSG00000004961
##     [4]     chr7   [27217606, 27217607]      * | ENSG00000005073
##     [5]    chr19   [ 2278618,  2278619]      * | ENSG00000005206
##     ...      ...                    ...    ... .             ...
##   [309]     chr7 [ 42001667,  42001668]      * | ENSG00000261019
##   [310]    chr16 [ 72974588,  72974589]      * | ENSG00000261227
##   [311]    chr17 [ 36070600,  36070601]      * | ENSG00000274620
##   [312]    chr10 [126211892, 126211893]      * | ENSG00000278831
##   [313]    chr10 [126211704, 126211705]      * | ENSG00000278831
##               Probe      Symbol       Sides        Raw.p          FDR
##         <character> <character> <character>    <numeric>    <numeric>
##     [1]  cg16474725       ENPP4          R6 1.270122e-21 1.939775e-20
##     [2]  cg27061889       LASP1          R1 7.105825e-20 7.827895e-19
##     [3]  cg25079603        HCCS          L9 1.248496e-05 1.819473e-05
##     [4]  cg15864691      HOXA11          R4 3.044005e-11 9.708257e-11
##     [5]  cg27341866      SPPL2B          R5 3.102233e-09 7.682987e-09
##     ...         ...         ...         ...          ...          ...
##   [309]  cg02901002  AC010132.4          R6 1.512426e-04 1.803295e-04
##   [310]  cg26909217  AC140912.1          R7 1.431867e-05 2.058754e-05
##   [311]  cg04136369     MIR378J          L2 7.562441e-05 9.524983e-05
##   [312]  cg16426764  AL513190.1          R3 9.809833e-06 1.445882e-05
##   [313]  cg03499943  AL513190.1          R3 9.809833e-06 1.445882e-05
##         distNearestTSS DMC_analysis_pvalue DMC_analysis_escc_Minus_esad
##              <integer>           <numeric>                    <numeric>
##     [1]         118379        2.576139e-23                    0.3499644
##     [2]           4741        1.656185e-36                    0.3364163
##     [3]         643982        5.973122e-28                    0.3930559
##     [4]           7228        1.075144e-29                    0.3729246
##     [5]          49994        1.269338e-11                    0.3082182
##     ...            ...                 ...                          ...
##   [309]         940569        2.481712e-29                    0.3590680
##   [310]         291364        1.498781e-27                    0.3326685
##   [311]          95552        4.455193e-24                    0.3275653
##   [312]         100028        6.455562e-32                    0.3729131
##   [313]         100216        9.619440e-24                    0.3195927
##         DMC_analysis_adjust.p external_gene_name original_ensembl_gene_id
##                     <numeric>        <character>              <character>
##     [1]          3.119777e-21              ENPP4        ENSG00000001561.6
##     [2]          1.227075e-33              LASP1       ENSG00000002834.16
##     [3]          1.348872e-25               HCCS       ENSG00000004961.13
##     [4]          3.036212e-27             HOXA11        ENSG00000005073.5
##     [5]          2.967588e-10             SPPL2B       ENSG00000005206.15
##     ...                   ...                ...                      ...
##   [309]          6.706775e-27         AC010132.4        ENSG00000261019.1
##   [310]          3.213995e-25         AC140912.1        ENSG00000261227.2
##   [311]          5.983252e-22            MIR378J        ENSG00000274620.1
##   [312]          2.425544e-29         AL513190.1        ENSG00000278831.1
##   [313]          1.242528e-21         AL513190.1        ENSG00000278831.1
##                                                                       TF
##                                                              <character>
##     [1]                 HNF4A;HNF4G;ELF3;PPARG;FOXD2;NR1H3;NR1H4;ZSCAN16
##     [2]                 HNF4A;HNF4G;FOXA3;FOXA2;FOXD2;ELF3;NR3C2;ZSCAN16
##     [3]                                                     ELF3;ZSCAN16
##     [4]                             FOXA3;FOXA2;FOXD2;NR2F2;ELF3;ZSCAN16
##     [5]                         HNF4A;HNF4G;PPARG;NR1H3;NR1H4;ELF3;NR2F2
##     ...                                                              ...
##   [309]                             GATA6;GATA4;ELF3;FOXD2;MECOM;ZSCAN16
##   [310]             HNF1A;HNF1B;FOXD2;PPARG;ELF3;NR1H3;NR1H4;NR3C2;NR2F2
##   [311]                 GATA6;GATA4;FOXA3;FOXA2;ELF3;FOXD2;MECOM;ZSCAN16
##   [312] HNF4A;HNF4G;FOXA3;FOXA2;FOXD2;PPARG;ELF3;NR1H3;NR1H4;NR2F2;NR3C2
##   [313]                                     FOXA3;FOXA2;FOXD2;ELF3;NR3C2
##         probe_hg19_width gene_hg19_seqnames gene_hg19_start gene_hg19_end
##                <integer>        <character>       <integer>     <integer>
##     [1]                2               chr6        46097730      46114436
##     [2]                2              chr17        37026112      37078023
##     [3]                2               chrX        11129421      11141198
##     [4]                2               chr7        27221129      27224842
##     [5]                2              chr19         2328614       2340154
##     ...              ...                ...             ...           ...
##   [309]                2               chr7        42940871      42942238
##   [310]                2              chr16        73265954      73267869
##   [311]                2              chr17        35974976      35975084
##   [312]                2              chr10       126311922     126312648
##   [313]                2              chr10       126311922     126312648
##         gene_hg19_width gene_hg19_strand
##               <integer>      <character>
##     [1]           16707                +
##     [2]           51912                +
##     [3]           11778                +
##     [4]            3714                -
##     [5]           11541                +
##     ...             ...              ...
##   [309]            1368                -
##   [310]            1916                +
##   [311]             109                -
##   [312]             727                +
##   [313]             727                +
##   -------
##   seqinfo: 22 sequences from an unspecified genome; no seqlengths

probes_extended_500bp <- resize(unique(probes),width = 501,fix = "center")
probes_extended_500bp[,1:2]

## GRanges object with 313 ranges and 2 metadata columns:
##         seqnames                 ranges strand |          GeneID
##            <Rle>              <IRanges>  <Rle> |     <character>
##     [1]     chr6   [45979099, 45979599]      * | ENSG00000001561
##     [2]    chr17   [37034634, 37035134]      * | ENSG00000002834
##     [3]     chrX   [11773216, 11773716]      * | ENSG00000004961
##     [4]     chr7   [27217356, 27217856]      * | ENSG00000005073
##     [5]    chr19   [ 2278368,  2278868]      * | ENSG00000005206
##     ...      ...                    ...    ... .             ...
##   [309]     chr7 [ 42001417,  42001917]      * | ENSG00000261019
##   [310]    chr16 [ 72974338,  72974838]      * | ENSG00000261227
##   [311]    chr17 [ 36070350,  36070850]      * | ENSG00000274620
##   [312]    chr10 [126211642, 126212142]      * | ENSG00000278831
##   [313]    chr10 [126211454, 126211954]      * | ENSG00000278831
##               Probe
##         <character>
##     [1]  cg16474725
##     [2]  cg27061889
##     [3]  cg25079603
##     [4]  cg15864691
##     [5]  cg27341866
##     ...         ...
##   [309]  cg02901002
##   [310]  cg26909217
##   [311]  cg04136369
##   [312]  cg16426764
##   [313]  cg03499943
##   -------
##   seqinfo: 22 sequences from an unspecified genome; no seqlengths

unique(width(probes_extended_500bp[,1:2]))

## [1] 501

hits <- suppressWarnings(findOverlaps(peaks.gr,probes_extended_500bp))
probes.hit <- subjectHits(hits)
length(unique(probes_extended_500bp[probes.hit])) # unique regions

## [1] 45

4.1 Enrichement between peaks and distal probes

x <- matrix(c(
  length(unique(probes[probes.hit])), # all linked probes with hits
  (length(unique(probes)) - length(unique(probes[probes.hit]))), # all linked probes with no hits
  length(unique(distal.probes.hit)), # all distal probes with  hits (4216)
  (length(distal.probes)) - length(unique(distal.probes.hit))), # all distal probes with no hits
  nrow = 2,
  dimnames = list(c("Overlap", "No overlap"),
                  c("Unique Linked probes (subfamily)", "Distal probes"))
)
x <- t(x)
x

##                                  Overlap No overlap
## Unique Linked probes (subfamily)      45        268
## Distal probes                       1778     161798

fisher.test(x)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  10.84653 21.09489
## sample estimates:
## odds ratio 
##   15.28118

4.2 Enrichement between peaks and distal probes with signature for subfamily members

distal.probes.subfamily.motif <- distal.probes.extended_500bp[names(distal.probes.extended_500bp) %in% probes.subfamily]
hits.distal.subfamily.motif <- suppressWarnings(findOverlaps(peaks.gr,distal.probes.subfamily.motif))
distal.probes.subfamily.motif.hit <- subjectHits(hits.distal.subfamily.motif)

x <- matrix(c(
  length(unique(probes[probes.hit])), # all linked probes with hits
  (length(unique(probes)) - length(unique(probes[probes.hit]))), # all linked probes with no hits
  length(unique(distal.probes.subfamily.motif.hit)), # all distal probes with  hits (4216)
  (length(distal.probes.subfamily.motif)) - length(unique(distal.probes.subfamily.motif.hit))), # all distal probes with no hits
  nrow = 2,
  dimnames = list(c("Overlap", "No overlap"),
                  c("Unique Linked probes", "Distal probes with subfamily motif signature"))
)
x <- t(x)
x

##                                              Overlap No overlap
## Unique Linked probes                              45        268
## Distal probes with subfamily motif signature     981      41065

fisher.test(x)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  4.975279 9.736710
## sample estimates:
## odds ratio 
##   7.027134

4.3 Enrichement between peaks and distal probes without signature for subfamily members

distal.probes.subfamily.nomotif <- distal.probes.extended_500bp[!names(distal.probes) %in% probes.subfamily]
hits.distal.subfamily.nomotif <- suppressWarnings(findOverlaps(peaks.gr,distal.probes.subfamily.nomotif))
distal.probes.subfamily.nomotif.hit <- subjectHits(hits.distal.subfamily.nomotif)

x <- matrix(c(
  length(unique(probes[probes.hit])), # all linked probes with hits
  (length(unique(probes)) - length(unique(probes[probes.hit]))), # all linked probes with no hits
  length(unique(distal.probes.subfamily.nomotif.hit)), # all distal probes with  hits (4216)
  (length(distal.probes.subfamily.nomotif)) - length(unique(distal.probes.subfamily.nomotif.hit))), # all distal probes with no hits
  nrow = 2,
  dimnames = list(c("Overlap", "No overlap"),
                  c("Unique Linked probes", "Distal probes without subfamily motif signature"))
)
x <- t(x)
x

##                                                 Overlap No overlap
## Unique Linked probes                                 45        268
## Distal probes without subfamily motif signature     797     120733

fisher.test(x)

## 
##  Fisher's Exact Test for Count Data
## 
## data:  x
## p-value < 2.2e-16
## alternative hypothesis: true odds ratio is not equal to 1
## 95 percent confidence interval:
##  17.97910 35.30655
## sample estimates:
## odds ratio 
##   25.41648

4.4 Summary

summary.df <- as.data.frame(unique(summary.df))
summary.df$ratio <- paste0(
  round(((summary.df$Overlap/apply(summary.df,1,
                                   function(y) {
                                     sum(as.integer(y[1]),as.integer(y[2]))
                                   })) * 100),
        digits = 2),"%")
summary.df

5 Enrichement plot

names(probes_extended_500bp) <- probes_extended_500bp$Probe
p <- probes_extended_500bp[targets$Probe]

mat1 <- normalizeToMatrix(peaks.gr, 
                          p, 
                          value_column = "peak",
                          extend = 5000, 
                          mean_mode = "w0", 
                          w = 50)


set.seed(1234)
chip.heatmap <- EnrichedHeatmap(mat1,
                                #km = 2, # Kmeans cluster
                                #row_title_rot = 0,
                                #cluster_rows = TRUE, 
                                #top_annotation = HeatmapAnnotation(enriched = anno_enriched(gp = gpar(col = 2:4, lty = 1:3))),
                                #split = sample(c("A", "B"), 455, replace = TRUE),
                                name = "ELF3 TF", 
                                col = c("white", "red"))


exp <- read_tsv("/hdd/ESCA/RNA-seq/ELF3/results_hg19/stringtieFPKM/Eso26_shELF3-2_FT-SA08056_S11_L001_1Aligned.sortedByCoord.out.gene_abund.txt",col_types = readr::cols())
exp$`Gene ID` <- sapply(exp$`Gene ID`, function(x) unlist(stringr::str_split(x,"\\."))[1])
exp <- exp[match(targets$GeneID,exp$`Gene ID`),]

exp.heatmap <- Heatmap(log2(exp$TPM + 1), 
                       col = c("white", "orange"), 
                       name = "Targets log2(TPM + 1) - Eso26shELF3", 
                       show_row_names = FALSE, 
                       width = unit(5, "mm")
) 

exp2 <- read_tsv("/hdd/ESCA/RNA-seq/ELF3/results_hg19/stringtieFPKM/Eso26_shCtrl_FT-SA08055_S10_L001_1Aligned.sortedByCoord.out.gene_abund.txt",col_types = readr::cols())
exp2$`Gene ID` <- sapply(exp2$`Gene ID`, function(x) unlist(stringr::str_split(x,"\\."))[1])
exp2 <- exp2[match(targets$GeneID,exp$`Gene ID`),]

exp.heatmap2 <- Heatmap(log2(exp2$TPM + 1), 
                        col = c("white", "orange"), 
                        name = "Targets log2(TPM + 1) - Eso26shCtrl", 
                        show_row_names = FALSE, 
                        width = unit(5, "mm")
) 


meth_col_fun = circlize::colorRamp2(c(min(p$DMC_analysis_escc_Minus_esad), 
                                      max(p$DMC_analysis_escc_Minus_esad)), 
                                    c("white", "red"))
dna.met.heatmap <- Heatmap(p$DMC_analysis_escc_Minus_esad,
                           name = "Diff DNA met (escc - esad)",
                           col = meth_col_fun,
                           width = unit(5, "mm"))

pair.pvalue.heatmap <- Heatmap(-log10(p$FDR),
                               name = "-log10(Pair FDR)",
                               width = unit(5, "mm"))

chip.heatmap + exp.heatmap + exp.heatmap2 + dna.met.heatmap + pair.pvalue.heatmap

6 RNA-seq aligned to hg19

6.1 Subfamily targets

load("ELF3_RNA_seq_DEA.rda")
dea <- results(dds)
dea

## log2 fold change (MLE): sample.group Eso26 shELF vs Eso26 shCtrl F 
## Wald test p-value: sample.group Eso26 shELF vs Eso26 shCtrl F 
## DataFrame with 22213 rows and 6 columns
##                   baseMean log2FoldChange     lfcSE        stat     pvalue
##                  <numeric>      <numeric> <numeric>   <numeric>  <numeric>
## ENSG00000000003  787.21799      0.9688876 0.5538892  1.74924435 0.08024879
## ENSG00000000419 1108.33237      0.1636841 0.4527124  0.36156314 0.71767851
## ENSG00000000457  229.31015     -0.2494749 0.5657434 -0.44096839 0.65923588
## ENSG00000000460  172.24987      0.6017720 0.6313255  0.95318808 0.34049478
## ENSG00000000938    1.00223     -0.1926380 5.4033485 -0.03565161 0.97156016
## ...                    ...            ...       ...         ...        ...
## ENSG00000273483   1.002230     -0.1926380  5.403349 -0.03565161  0.9715602
## ENSG00000273486  17.171535      0.7459568  1.449597  0.51459588  0.6068354
## ENSG00000273488  11.893125      0.2927859  1.687949  0.17345663  0.8622925
## ENSG00000273492   1.069045     -3.6351986  5.850612 -0.62133643  0.5343783
## ENSG00000273493   1.403122      3.8348712  5.237585  0.73218311  0.4640568
##                      padj
##                 <numeric>
## ENSG00000000003  0.998371
## ENSG00000000419  0.998371
## ENSG00000000457  0.998371
## ENSG00000000460  0.998371
## ENSG00000000938  0.998371
## ...                   ...
## ENSG00000273483  0.998371
## ENSG00000273486  0.998371
## ENSG00000273488  0.998371
## ENSG00000273492  0.998371
## ENSG00000273493  0.998371

x.cut <- 1
y.cut <- 0.05

dea$group <- "Not Significant"
dea[which(dea$pvalue < y.cut & dea$log2FoldChange < -x.cut ),"group"] <- "Downregulated"
dea[which(dea$pvalue < y.cut & dea$log2FoldChange > x.cut ),"group"] <- "Upregulated"

pairs <- read_csv("~/paper_elmer/ESCA_analysis_hg38/escc_vs_esad_hg38/hyper/getTFtargets_genomic_coordinates_mapped_to_hg19.hyper.subfamily.csv",
                  col_types = readr::cols())
targets <- pairs[grep("ELF3",pairs$TF),]

idx <- rownames(dea) %in% unique(targets$GeneID)
dea[idx,"group"] <- paste0(dea[idx,"group"]," - ELMER infered")
table(as.data.frame(dea)$group)

## 
##                   Downregulated   Downregulated - ELMER infered 
##                              91                               6 
##                 Not Significant Not Significant - ELMER infered 
##                           21179                             692 
##                     Upregulated     Upregulated - ELMER infered 
##                             236                               9

f <- list(
  family = "Courier New, monospace",
  size = 18,
  color = "#7f7f7f"
)
x <- list(
  title = "log2FoldChange",
  titlefont = f
)
y <- list(
  title = "-log10(raw p-value)",
  titlefont = f
)

library(plotly)
p <- plot_ly(data =  as.data.frame(dea), 
             x = dea$log2FoldChange, 
             y = -log10(dea$pvalue), 
             text = rowData$Gene_Name[match(rownames(dea),rowData$Gene_ID)], 
             mode = "markers", 
             color = dea$group) %>% 
  layout(title ="Volcano Plot") %>%
  layout(xaxis = x, yaxis = y)  %>%
  layout(shapes=list(list(type='line', 
                          x0 = x.cut, 
                          x1 = x.cut, 
                          y0 = 0, 
                          y1 = max(-log10(dea$padj),na.rm = T), 
                          line=list(dash='dot', width=1)),
                     list(type='line', 
                          x0 = -x.cut, 
                          x1 = -x.cut, 
                          y0 = 0, 
                          y1 = max(-log10(dea$padj),na.rm = T), 
                          line =list(dash='dot', width=1)),
                     list(type ='line', 
                          x0 = min(dea$log2FoldChange), 
                          x1 = max(dea$log2FoldChange), 
                          y0 =  -log10(y.cut), 
                          y1 =  -log10(y.cut), 
                          line = list(dash='dot', width=1))
  ) 
  )
p

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

6.2 Family targets

dea$group2 <- "Not Significant"
dea[which(dea$pvalue < y.cut & dea$log2FoldChange < -x.cut ),"group2"] <- "Downregulated"
dea[which(dea$pvalue < y.cut & dea$log2FoldChange > x.cut ),"group2"] <- "Upregulated"

pairs <- read_csv("~/paper_elmer/ESCA_analysis_hg38/escc_vs_esad_hg38/hyper/getTFtargets_genomic_coordinates_mapped_to_hg19.hyper.family.csv",
                  col_types = readr::cols())
targets <- pairs[grep("ELF3",pairs$TF),]

idx <- rownames(dea) %in% unique(targets$GeneID)
dea[idx,"group2"] <- paste0(dea[idx,"group2"]," - ELMER infered")
DT::datatable(as.data.frame(dea),options = list(scrollX = TRUE),filter = 'top')

## Warning in instance$preRenderHook(instance): It seems your data is too
## big for client-side DataTables. You may consider server-side processing:
## https://rstudio.github.io/DT/server.html

table(as.data.frame(dea)$group2)

## 
##                   Downregulated   Downregulated - ELMER infered 
##                              87                              10 
##                 Not Significant Not Significant - ELMER infered 
##                           20684                            1187 
##                     Upregulated     Upregulated - ELMER infered 
##                             229                              16

f <- list(
  family = "Courier New, monospace",
  size = 18,
  color = "#7f7f7f"
)
x <- list(
  title = "log2FoldChange",
  titlefont = f
)
y <- list(
  title = "-log10(raw p-value)",
  titlefont = f
)

library(plotly)
p <- plot_ly(data =  as.data.frame(dea), 
             x = dea$log2FoldChange, 
             y = -log10(dea$pvalue), 
             text = rowData$Gene_Name[match(rownames(dea),rowData$Gene_ID)], 
             mode = "markers", 
             color = dea$group) %>% 
  layout(title ="Volcano Plot") %>%
  layout(xaxis = x, yaxis = y)  %>%
  layout(shapes=list(list(type='line', 
                          x0 = x.cut, 
                          x1 = x.cut, 
                          y0 = 0, 
                          y1 = max(-log10(dea$padj),na.rm = T), 
                          line=list(dash='dot', width=1)),
                     list(type='line', 
                          x0 = -x.cut, 
                          x1 = -x.cut, 
                          y0 = 0, 
                          y1 = max(-log10(dea$padj),na.rm = T), 
                          line =list(dash='dot', width=1)),
                     list(type ='line', 
                          x0 = min(dea$log2FoldChange), 
                          x1 = max(dea$log2FoldChange), 
                          y0 =  -log10(y.cut), 
                          y1 =  -log10(y.cut), 
                          line = list(dash='dot', width=1))
  ) 
  )
p

## No trace type specified:
##   Based on info supplied, a 'scatter' trace seems appropriate.
##   Read more about this trace type -> https://plot.ly/r/reference/#scatter

Evaluation enrichment of ELF3 targets

1 Summary

1.1 ELF3 family and subfamily

1.2 Enriched motifs and probes associated

2 Narrow ELF3 peaks

3 Comparing paired probes considering family classification

3.1 Get ELF3 probes target regions

3.2 Distance between probes and peaks

3.3 Results overlap

3.3.1 Overlap between probes (+-250bp) and peaks

3.3.2 Overlap between peaks and distall probes (+-250bp)

3.3.3 Regions overlapped

3.3.4 Enrichement between peaks and all distal probes

3.3.5 Enrichement between peaks and distal probes with signature for family members

3.3.6 Enrichement between peaks and distal probes without signature for family members

4 Comparing paired probes considering subfamily classification

4.1 Enrichement between peaks and distal probes

4.2 Enrichement between peaks and distal probes with signature for subfamily members

4.3 Enrichement between peaks and distal probes without signature for subfamily members

4.4 Summary

5 Enrichement plot

6 RNA-seq aligned to hg19

6.1 Subfamily targets

6.2 Family targets