I don’t see how this would be possible, but let’s have a loop at the operations:
# sequential
y1 = x * w1 + b1 # linear
y2 = (y1 - running_mean) / sqrt(running_var + eps) * gamma + beta # batchnorm

# replace y1
y2 = (x * w1 + b1 - running_mean) / sqrt(running_var + eps) * gamma + beta

If you multiply gamma into the left brackets, you would mix the running_mean with a trainable parameter, so how should this “fusion” look like?