There is a field of research that can be identified by the keywords "mimetic discretizations" and "compatible spatial discretizations" that starts from your same assumption: algebraic equations arising from the discretization of partial differential equations should mimic the properties of the continuum operators.
A good starting point could be Principles of Mimetic Discretizations of Differential Operators by Bochev and Hyman, which appeared in this book.
It is an interesting framework, but I would not focus to much on "matrix properties" per se: finally what matters is obtaining a "good solution" (robust, accurate, not to much expensive, ...) and not having a nice matrix at hand.
I would suggest you to start asking more narrow questions: e.g. how can I deal efficiently with the unsymmetric terms? Can I trade accuracy with speed by using a symmetric solver?
Addendum
Jed Brown correctly points out that sometimes lack of symmetry is only apparent.
Let me explain this point with a very trivial example. Suppose that $\mathcal{L}u=f$ + b.c (continuous, with $\mathcal{L}$ self adjoint) becomes $Ax=b$ (discrete, with $A=A^T$) where the b.c. are expressed as $u = Bv$ (see this answer for a concrete example, sorry for the shameless self citation). Now by substitution we have $ABv = b$ where $AB$ is not symmetric, but symmetry can be restored by premultiplication by $B^T$ to obtain $B^TABv = B^Tf$. If matrix $B$ is square (and note that usually it is an identity matrix plus some low rank correction), problem $ABv=f$ is perfectly equivalent to $B^TABv = B^Tf$. If premultiplication by $B^T$ is not feasible, sometimes a symmetric problem can be obtained by augmenting the system, i.e. by introducing extra unknown variables.
[disclaimer: what follows is only a conjecture since I had no time to work out your problem in detail]
In the case at hand b.c. are introduced by adding ghost points outside the problem domain: unknown fields at the ghost points are obtained by extrapolation from the domain points, taking into account boundary geometry and b.c. (The algorithm actually goes the other way: from the ghost point an image point in the problem domain is derived, and its value is obtained by interpolation of the surrounding problem points. For this discussion however the approaches are equivalent.) The source of "unsymmetry" is clear, and the situation is slightly different from my trivial example above, so that there is no easy recipe for fixing the symmetry.
My suggestion is to investigate in these two directions:
Is this a sensible way of introducing the boundary conditions? If the answer is yes, lack of symmetry is acceptable at the formulation level.
Is lack of symmetry acceptable at the computational level? I.e. are you able to solve the unsymmetric problem? If no, try to modify the system for restoring symmetry in a way that does not introduces further approximations (like the trivial example above) or trade accuracy with speed by introducing some extra approx.
