(* how Mathematica stores data *)
a+b // FullForm

3a // FullForm

3a+4(a+b) // FullForm

{a,3,4,b} // FullForm

Head[a]
Head[a+b]

(* rules act on FullForm *)
3a+4(a+b) /. 2a->b
3a+4(a+b) /. a->b/2
3a+4(a+b) /. x_ :>s /; x>3
List[a,3,4,b] /. List:> Plus

(* Metric : metric Tensor,
   GammaProd: head for product of gamma matrices
 Ga: head for gamma matrix 
 GammaProd[Ga[mu],Ga[nu]] means (γ^μ)(γ^ν) *)

(* metric is symmetric *)
SetAttributes[Metric,Orderless];

(* function that applies replacement rules repeatedly to its argument *)
GammaContract[arg_]:= arg //. {
(* contraction of two gamma matrices in a gamma product *)
GammaProd[a___,Ga[mu_],Ga[mu_],b___]:>4 GammaProd[a,b],

(* empty list: identity matrix *)
GammaProd[] :> GI,

(* order gamma matrices canonically *)
GammaProd[a___,Ga[mu_],Ga[nu_],b___] :> -GammaProd[a,Ga[nu],Ga[mu],b]+2 Metric[mu,nu] GammaProd[a,b] /;!OrderedQ[{mu,nu}],	
(* contract metric with gamma product *)
Metric[mu_,nu_]GammaProd[a___,Ga[mu_],b___]:> GammaProd[a,Ga[nu],b],

(* contract metric with metric *)
Metric[mu_,nu_] Metric[nu_,rh_] :> Metric[mu,rh],

(* contract metric with itself *)
Metric[mu_,mu_]:> 4
};

(* test *)
GammaProd[Ga[mu],Ga[mu]] // GammaContract

GammaProd[Ga[mu],Ga[mu],Ga[nu],Ga[nu]]// GammaContract

GammaProd[Ga[mu],Ga[nu],Ga[mu],Ga[nu]]// GammaContract

GammaProd[Ga[mu],Ga[nu],Ga[mu]]// GammaContract

GammaProd[Ga[mu],Ga[al],Ga[be],Ga[mu]] // GammaContract // Expand

GammaProd[Ga[mu],Ga[be],Ga[al],Ga[mu]] // GammaContract// Expand

GammaProd[Ga[mu],Ga[al],Ga[be],Ga[ga],Ga[mu]] // GammaContract// Expand

GammaProd[Ga[mu],Ga[al],Ga[be],Ga[ga],Ga[de],Ga[eps],Ga[mu]] // GammaContract// Expand