A new generation of neutrinoless double beta decay (0νββ-decay) experiments with improved sensitivity is currently under design and construction. They will probe inverted hierarchy region of the neutrino mass pattern. There is also a revived interest to the resonant neutrinoless doubleelectron capture (0νECEC), which has also a potential to probe lepton number conservation and to investigate the neutrino nature and mass scale. The primary concern are the nuclear matrix elements. Clearly, the accuracy of the determination of the effective Majorana neutrino mass from the measured 0νββ-decay half-life is mainly determined by our knowledge of the nuclear matrix elements. We review recent progress achieved in the calculation of 0νββ and 0νECEC nuclear matrix elements within the quasiparticle random phase approximation. A considered self-consistent approach allow to derive the pairing, residual interactions and the two-nucleon short-range correlations from the same modern realistic nucleon-nucleon potentials. The effect of nuclear deformation is taken into account. A possibility to evaluate 0νββ-decay matrix elements phenomenologically is discussed.