In eukaryotes, gene lengths is extremely variable, ranging from the hundred of base pairs (i.e. the minimal resolution of ChIP-seq) to around a million of base pair.
When performing a metagene analysis, 100-base pair long genes are stretch to the same length as million base pair genes, which can create issues, i.e. promoter signals will "spill out" to gene body and gene ends for short genes, but not for long genes.
For example, when comparing 2 transcription factors usually binding to gene promoters: if the first transcription factor (TF1) is associated with short genes while the second (TF2) is associated to longer genes, metagene profiles can give the impression than TF1 is more associated to gene body and gene ends than TF2.
Metagene profiles are also sometimes affected by the Simpson Paradox. Therefore I would rather recommand avoiding metagene profiles. One could produces several plots with anchors set at transcription start sites, transcription end sites, and middle exon starts, and if possible by filtering short genes as an alternative.
I would recommand EnrichedHeatmap to produce such plots (it can also do metagene profiles, so here you go!), or if you want to take some risks, you can test our recent package epistack.
deepTools is a popular non-R alternative.