Formatting plot of symmetric points differently to asymmetric points

Question

From an array of numbers, I would like to determine (and then colour appropriately) which points are symmetric about a given axis (here, $x = 1/2$), and which are not.

For example, form this:

ListPlot[{{1.5,1.5},{1,1},{1.5,-0.5},{-0.5,-0.5},{-0.5,1},{-0.5,1.5}}, 
PlotStyle -> Red]

I would like to achieve this:

Plot[{0}, {x, -2, 2}, PlotRange -> {-2, 2}, Epilog -> 
{PointSize[Medium], Red,Point[{{{1.5, 1.5}, {-0.5, 1.5}, {1.5, -0.5}, {-0.5, -0.5}}}],
PointSize[Medium], Blue, Point[{{-0.5, 1}, {1, 1}}]}]

(Please excuse the clumsiness with which I code: I am very new to Mathematica.)

I have looked at Symmetric in Wolfram Documentation, but it doesn't seem to cover this. I don't know of a method for selecting certain arrays and formatting them with one colour, then selecting a different set and formatting that with another colour, other than by doing it manually. Is there a way of automating this process? (My goal is to do this for arrays of over 1000 points, so sorting manually is not really an option!)

Update2

Data: array link

Update3

@Kuba, Using the code you provided in answer to my other post:

sol = Solve[N[Table[BernoulliB[n, z], {n, 1000, 1000}] == 0]];

would it be possible to find symmetric points?

(And also format the plot: AspectRatio -> Automatic, Plot Axis @ x = 1/2, etc.?)

Update4

Note: So as not to mislead, I have included a corrected plot. The strange pattern produced without precision adjustment:

is purely caused by numerical error (thanks to Antonio Vargas for pointing this out) - the corrected plot for n = 171 (n = 1000 is too much for my machine!), calculated to greater precision, looks like this:

Plotted using Kuba's code:

sol = Solve[N[Table[BernoulliB[n, z], {n, 171, 171}] == 0, 50]];
point[x_?VectorQ] := {PointSize@.01, Red, Point[x]}
point[x : (_?VectorQ ..)] := {PointSize@.01, Blue, Point /@ {x}}
axis = 1/2;
point @@@ GatherBy[{Re@z, Im@z} /. sol,(*Round added*)
Round[{Abs[#[[1]] - axis], #[[2]]}, 0.001] &] // 
Graphics[#, Frame -> True, GridLines -> {{axis}, None}] &

Answer 1

list = DeleteDuplicates@RandomInteger[{-10, 10}, {155, 2}];
axis = 1/2;

First let's gather points with the same distance form your axis and the same second coordinate.

Those which have simmetric partner will appear in pair so we can gather by length to distinguish singles from pairs.

pts = Join @@@ GatherBy[
               GatherBy[list, {Abs[#[[1]] - axis], #[[2]]} &]
               , Length]

Graphics[{PointSize[0.05],
          Table[{RGBColor @@ RandomReal[1, 3], Point@pts[[ i]]}, {i, Length@pts}]
         }, GridLines -> {{axis}, None}, Frame -> True]

---

Edit respond to OP's edit. I've also included b.gatessucs suggestion to skip second GatherBy. This is slightly different method.

sol = Solve[N[Table[BernoulliB[n, z], {n, 1000, 1000}] == 0]];

point[x_?VectorQ] := {PointSize@.01, Red, Point[x]}
point[x : (_?VectorQ ..)] := {PointSize@.03, Blue, Point /@ {x}}

axis = 1/2;

point @@@ GatherBy[{Re@z, Im@z} /. sol, {Abs[#[[1]] - axis], #[[2]]} &
                  ] // Graphics[#, Frame -> True, GridLines -> {{axis}, None}] &

There are not so many symmetric pairs because points have real coordinates, you can add Round if you can/want, in order to make the criterium weaker:

point @@@ GatherBy[{Re@z, Im@z} /. sol, 
 (*Round added*)   Round[{Abs[#[[1]] - axis], #[[2]]}, .1] & 
                  ] // Graphics[#, Frame -> True, GridLines -> {{axis}, None}] &

Look out with Round to not catch nearby points as symmetric.