MOTIVATION: Recent advances in DNA sequencing technologies have allowed the detailed characterization of genomes in large cohorts of tumors, highlighting their extreme heterogeneity, with no two tumors sharing the same complement of somatic mutations. Such heterogeneity hinders our ability to identify somatic mutations important for the disease, including mutations that determine clinically relevant phenotypes (e.g. cancer subtypes). Several tools have been developed to identify somatic mutations related to cancer phenotypes. However, such tools identify correlations between somatic mutations and cancer phenotypes, with no guarantee of highlighting causal relations. RESULTS: We describe ALLSTAR, a novel tool to infer reliable causal relations between somatic mutations and cancer phenotypes. ALLSTAR identifies reliable causal rules highlighting combinations of somatic mutations with the highest impact in terms of average effect on the phenotype. While we prove that the underlying computational problem is NP-hard, we develop a branch-and-bound approach that employs protein-protein interaction networks and novel bounds for pruning the search space, while properly correcting for multiple hypothesis testing. Our extensive experimental evaluation on synthetic data shows that our tool is able to identify reliable causal relations in large cancer cohorts. Moreover, the reliable causal rules identified by our tool in cancer data show that our approach identifies several somatic mutations known to be relevant for cancer phenotypes as well as novel biologically meaningful relations. AVAILABILITY AND IMPLEMENTATION: Code, data, and scripts to reproduce the experiments available at https://github.com/VandinLab/ALLSTAR.