Treatment-resistant depression (TRD) occurs in ~30% of patients with major depressive disorder (MDD) but the genetics of TRD was previously poorly investigated. Whole exome sequencing and genome-wide genotyping were available in 1209 MDD patients after quality control. Antidepressant response was compared to non-response to one treatment and non-response to two or more treatments (TRD). Differences in the risk of carrying damaging variants were tested. A score expressing the burden of variants in genes and pathways was calculated weighting each variant for its functional (Eigen) score and frequency. Gene-based and pathway-based scores were used to develop predictive models of TRD and non-response using gradient boosting in 70% of the sample (training) which were tested in the remaining 30% (testing), evaluating also the addition of clinical predictors. Independent replication was tested in STAR*D and GENDEP using exome array-based data. TRD and non-responders did not show higher risk to carry damaging variants compared to responders. Genes/pathways associated with TRD included those modulating cell survival and proliferation, neurodegeneration, and immune response. Genetic models showed significant prediction of TRD vs. response and they were improved by the addition of clinical predictors, but they were not significantly better than clinical predictors alone. Replication results were driven by clinical factors, except for a model developed in subjects treated with serotonergic antidepressants, which showed a clear improvement in prediction at the extremes of the genetic score distribution in STAR*D. These results suggested relevant biological mechanisms implicated in TRD and a new methodological approach to the prediction of TRD.