In this study, we aimed to investigate the association of four single nucleotide polymorphisms (SNPs) (MTHFR 677 C > T, MTHFR 1298 A > C, MTR 2756 A > G and MTRR 66 A > G), gene-gene interaction and haplotype combination with pulmonary embolism (PE) risk based on Chinese Han population. Logistic regression was performed to investigate association between four SNPs within folate metabolism gene and PE risk, and GMDR model was used to investigate the additional gene-gene interactions among the four SNPs. Logistic analysis showed that rs1801133 and rs1801131 in MTHFR gene were associated with increased PE risk in both additive and dominant models. The carriers with homozygous mutant of rs1801133 polymorphism and homozygous of rs1801131 were associated with increased PE risk, and ORs (95% CI) were 1.71(1.24-2.21) and 1.58 (1.24-2.01), respectively. We also found a significant gene-gene interaction between rs1801133 and rs1801131 on PE. Overall, the cross-validation consistency of this two-locus model was 10/10, and the testing accuracy was 60.72%, after adjusting for covariates. Haplotype containing the rs1801133- T and rs1801131- C alleles were associated with a statistically increased PE risk, OR (95% CI) = 2.68 (1.28-4.13), P < 0.001. We found that rs1801133 and rs1801131 within MTHFR gene, their interaction, and haplotype containing the rs1801133- T and rs1801131- C alleles were all associated with PE risk.