Let $(M,g)$ be a compact Riemann surface with no boundary and $u=(u_1,...,u_n)$ be a solution of the following singular Liouville system: \begin{equation*} \Delta_g u_i+\sum_{j=1}^na_{ij}\rho_j(\frac{h_je^{u_j}}{\int_M h_j e^{u_j}dV_g}-\frac{1}{vol_g(M)})=\sum_{t=1}^N4\pi \gamma_t( \delta_{p_t}-\frac{1}{vol_g(M)}), \end{equation*} where $i=1,...,n$, $h_1,...,h_n$ are positive smooth functions, $p_1,...,p_N$ are distinct points on $M$, $\delta_{p_t}$ are Dirac masses, $\rho=(\rho_1,...,\rho_n)$ ($\rho_i\ge 0)$ and $(\gamma_{1},...,\gamma_{N})$ ($\gamma_{t}>-1$ ) are constant vectors. If the coefficient matrix $A=(a_{ij})_{n\times n}$ satisfies standard assumptions we identify a family of critical hyper-surfaces $\Gamma_k$ for $\rho=(\rho_1,..,\rho_n)$ so that a priori estimate of $u$ holds if $\rho$ is not on any of the $\Gamma_k$s. Thanks to the a priori estimate, a topological degree for $u$ is well defined for $\rho $ staying between every two consecutive $\Gamma_k$s. In this article we establish this degree counting formula which depends only on the Euler Characteristic of $M$ and the location of $\rho$. Finally if the Liouville system is defined on a bounded domain in $\mathbb R^2$ with Dirichlet boundary condition, a similar degree counting formula that depends only on the topology of the domain and the location of $\rho$ is also determined.