A computer-implemented digital image analysis method for co-expression analysis of multiple markers in a tissue sample comprising: accessing a plurality of single marker channel images, each of the plurality of single marker channel images being a digital image whose pixel intensity values correlate with the amount of a respective single marker in the tissue sample; computing a heat map of marker expression for each of the single marker channel images, wherein computing the heat map comprises applying a low pass filter to each of the plurality of single marker channel images; identifying one or more candidate regions of interest, ROIs, in each heat map of marker expression, wherein the identifying of the one or more candidate ROIs comprises: applying a local maximum filter to each heat map of marker expression for identifying the highest ranked local pixel intensity maxima in said heat map, selecting a predetermined number of sub-regions of the heat map surrounding a respective one of identified local intensity maxima by delineating a field of view having a set of N × N pixels around each of the selected predetermined number of highest ranked local intensity maxima, and using the delineated fields of view as the candidate ROIs of said heat map; computing overlay masks comprising the identified one or more candidate ROIs from each heat map of marker expression, an overlay mask being an image mask masking out all pixels not belonging to one of the identified candidate ROIs of the heat map the overlay mask was derived from; identifying one or more co-localized ROIs from the overlay masks by: identifying intersection areas of the candidate ROIs of two or more of the heat maps, mapping the identified intersection areas back to the two or more heat maps; and using the candidate ROIs or parts thereof which correspond to said mapped-back intersection areas as the one or more co-localized ROIs; and selectively analyzing pixel intensities of the co-localized ROIs of two or more of the heat maps for determining if the markers corresponding to the two or more heat maps are coexpressed in the tissue sample. accessing a plurality of single marker channel images, each of the plurality of single marker channel images being a digital image whose pixel intensity values correlate with the amount of a respective single marker in the tissue sample; computing a heat map of marker expression for each of the single marker channel images, wherein computing the heat map comprises applying a low pass filter to each of the plurality of single marker channel images; identifying one or more candidate regions of interest, ROIs, in each heat map of marker expression, wherein the identifying of the one or more candidate ROIs comprises: applying a local maximum filter to each heat map of marker expression for identifying the highest ranked local pixel intensity maxima in said heat map, selecting a predetermined number of sub-regions of the heat map surrounding a respective one of identified local intensity maxima by delineating a field of view having a set of N × N pixels around each of the selected predetermined number of highest ranked local intensity maxima, and using the delineated fields of view as the candidate ROIs of said heat map; applying a local maximum filter to each heat map of marker expression for identifying the highest ranked local pixel intensity maxima in said heat map, selecting a predetermined number of sub-regions of the heat map surrounding a respective one of identified local intensity maxima by delineating a field of view having a set of N × N pixels around each of the selected predetermined number of highest ranked local intensity maxima, and using the delineated fields of view as the candidate ROIs of said heat map; computing overlay masks comprising the identified one or more candidate ROIs from each heat map of marker expression, an overlay mask being an image mask masking out all pixels not belonging to one of the identified candidate ROIs of the heat map the overlay mask was derived from; identifying one or more co-localized ROIs from the overlay masks by: identifying intersection areas of the candidate ROIs of two or more of the heat maps, mapping the identified intersection areas back to the two or more heat maps; and identifying intersection areas of the candidate ROIs of two or more of the heat maps, mapping the identified intersection areas back to the two or more heat maps; and using the candidate ROIs or parts thereof which correspond to said mapped-back intersection areas as the one or more co-localized ROIs; and selectively analyzing pixel intensities of the co-localized ROIs of two or more of the heat maps for determining if the markers corresponding to the two or more heat maps are coexpressed in the tissue sample.