Looks like FeynCalc is the way to go.

Here's an example problem from an online course I'm taking:

Recalling that [a^,a^†]=1, calculate the commutator:

[a^,(a^†)n]

a is an annihilation operator and a dagger is a creation operator.

Mathematica 10.0 for Linux x86 (64-bit)
Copyright 1988-2014 Wolfram Research, Inc.

In[1]:= !!fordan
(* FeynCalc can be installed by
        Import["http://www.feyncalc.org/install.m"]
*)
(* load FeynCalc *)

Quiet@Needs@"HighEnergyPhysics`FeynCalc`";
(* FI just sets the output to FeynCalcInternalForm, kind of FullForm *)
FI;

(* remember in FeynCalc:  .  is used as non-commutative operator *)

DeclareNonCommutative[a,b,c];  (* declares a,b,c to be non commutative *)

test = a . ( 2*b - z c).a ;

(* while now we add a commutator  [a,c] = 1 *)

Commutator[a, c] =1;

Expand @ DotSimplify @ test


In[2]:= <<fordan
Loading FeynCalc from /home/rolfm/.Mathematica/Applications/HighEnergyPhysics
FeynCalc 8.2.0 Type ?FeynCalc for help or visit http://www.feyncalc.org/
$PrePrint is set to FeynCalcForm. Use FI and FC to change the display format.
Loading FeynArts, see www.feynarts.de for documentation
FeynArts 3.7 patched for use with FeynCalc

Out[2]= -(a*z) + 2*a . b . a - z*c . a . a

One possibility is to use DotSimplify from FeynCalc. Or you look at the source and use the relevant parts.