We have established a quantum mechanics (QM)/molecular mechanics (MM) hybrid method for calculating the Raman spectra of protein-bound cofactors using the alpha-subunit of C-phycocyanin containing a phycocyanobilin (PCB) chromophore as a test case. The PCB cofactor was described with density functional theory, whereas the protein matrix was treated with the CHARMM force field. The Hessian matrix of the QM region was built by taking into account bonded and nonbonded interactions with the protein environment and projected onto the internal coordinate space. Force constants were scaled with a global set of scaling factors, and the Raman intensities were computed using a finite-field method combined with a fourth-order differentiation algorithm for the calculation of the polarizability derivatives. In general, the QM/MM results provided a substantially improved description of the experimental resonance Raman (RR) spectra of the protein-bound cofactor compared to QM calculations of isolated PCB models in vacuo. The results allow the assessment of the effect of the protein-cofactor interactions on the RR spectra and reveal the potential and limitations of QM calculations on isolated tetrapyrroles for determining the chromophore structures in the various species and states of phytochromes for which three-dimensional structures are not available.