We investigate pfaffian trial wave functions with singlet and triplet pair orbitals by quantum Monte Carlo methods. We present mathematical identities and the key algebraic properties necessary for efficient evaluation of pfaffians. Following upon our previous study, 1 we explore the possibilities of expanding the wave function in linear combinations of pfaffians. We observe that molecular systems require much larger expansions than atomic systems and linear combinations of a few pfaffians lead to rather small gains in correlation energy. We also test the wave function based on fully-antisymmetrized product of independent pair orbitals. Despite its seemingly large variational potential, we do not observe additional gains in correlation energy. We find that pfaffians lead to substantial improvements in fermion nodes when compared to Hartree-Fock wave functions and exhibit the minimal number of two nodal domains in agreement with recent results on fermion nodes topology. We analyze the nodal structure differences of Hartree-Fock, pfaffian and essentially exact large-scale configuration interaction wave functions. Finally, we combine the recently proposed form of backflow correlations 2,3 with both determinantal and pfaffian based wave functions.