You might be able to use a hypergeometric test to calculate the odds of seeing co-occurrence of two TFs in promoters, as compared with co-occurrence of those two TFs across the entire genome (perhaps minus non-mappable regions or other regions where you would not expect to have countable TF binding sites).
In R, you could use phyper:
> phyper(a, b, c, d)
I think the following might work:
a = number of observations of co-occurrence of two TFs in 1k window over all promoters (number of times you see a pairing over all proximal promoters)
b = number of observations of co-occurrence of two TFs over 1k windows over whole genome (number of times you see a pairing over all whole-genome 1k windows)
c = number of observations of non-co-occurrence of two TFs over 1k windows over whole genome (number of times you don't see a pairing over all whole-genome 1k windows)
d = number of observations of co-occurrence and non co-occurrence of two TFs in 1k window over all promoters (total observations over promoters, i.e. total number of proximal promoters)
I might use bedops --chop to make 1k windows over the genome, bedops --difference to excise unmappable regions, bedops --range to make 1k windows upstream of the TSS (proximal promoters), and use of bedmap --echo-map-id with some awk scripting could help with counts of TF hit pairs-of-interest over proximal promoter regions, and counts of TF hit pairs-of-interest over the genome. Etc. Bringing counts into R should provide the p-value.
I'm not sure that comparing or counting kmers directly would work here, because DNA binding is imprecise and TFs will bind even if the region of DNA isn't a perfect match with the consensus sequence. That's why sequence logos have some bases at different heights, or information content levels, because the probability of seeing a particular base at a position in a binding site isn't 100%. If it was, we wouldn't need logos.