They're believed to be symmetric because quite often a normal approximation is used. This one works well enough in case p lies around 0.5. binom.test
on the other hand reports "exact" Clopper-Pearson intervals, which are based on the F distribution (see here for the exact formulas of both approaches). If we would implement the Clopper-Pearson interval in R it would be something like (see note):
Clopper.Pearson <- function(x, n, conf.level){
alpha <- (1 - conf.level) / 2
QF.l <- qf(1 - alpha, 2*n - 2*x + 2, 2*x)
QF.u <- qf(1 - alpha, 2*x + 2, 2*n - 2*x)
ll <- if (x == 0){
0
} else { x / ( x + (n-x+1)*QF.l ) }
uu <- if (x == 0){
0
} else { (x+1)*QF.u / ( n - x + (x+1)*QF.u ) }
return(c(ll, uu))
}
You see both in the link and in the implementation that the formula for the upper and the lower limit are completely different. The only case of a symmetric confidence interval is when p=0.5. Using the formulas from the link and taking into account that in this case $n = 2\times x$ it's easy to derive yourself how it comes.
I personally understood it better looking at the confidence intervals based on a logistic approach. Binomial data is generally modeled using a logit link function, defined as:
$${\rm logit}(x) = \log\! \bigg( \frac{x}{1-x} \bigg)$$
This link function "maps" the error term in a logistic regression to a normal distribution. As a consequence, confidence intervals in the logistic framework are symmetric around the logit values, much like in the classic linear regression framework. The logit transformation is used exactly to allow for using the whole normality-based theory around the linear regression.
After doing the inverse transformation:
$${\rm logit}^{-1}(x) = \frac{e^x}{1+e^{x}}$$
You get an asymmetric interval again. Now these confidence intervals are actually biased. Their coverage is not what you would expect, especially at the boundaries of the binomial distribution. Yet, as an illustration they show you why it is logic that a binomial distribution has asymmetric confidence intervals.
An example in R:
logit <- function(x){ log(x/(1-x)) }
inv.logit <- function(x){ exp(x)/(1+exp(x)) }
x <- c(0.2, 0.5, 0.8)
lx <- logit(x)
upper <- lx + 2
lower <- lx - 2
logxtab <- cbind(lx, upper, lower)
logxtab # the confidence intervals are symmetric by construction
xtab <- inv.logit(logxtab)
xtab # back transformation gives asymmetric confidence intervals
note : In fact, R uses the beta distribution, but this is completely equivalent and computationally a bit more efficient. The implementation in R is thus different from what I show here, but it gives exactly the same result.