I discuss methods of calculation of propagator diagrams
(massless, those of Heavy Quark Effective Theory, and massive on-shell diagrams)
up to 3 loops.
Integration-by-parts recurrence relations are used to reduce them
to linear combinations of basis integrals.
Non-trivial basis integrals have to be calculated by some other method,
e.g., using Gegenbauer polynomial technique.
Many of them are expressed via hypergeometric functions;
in the massless and HQET cases, their indices tend to integers
at $\varepsilon\to0$.
I discuss the algorithm of their expansion in $\varepsilon$,
in terms of multiple $\zeta$ values.
These lectures were given at Calc-03 school, Dubna, 14–20 June 2003.