According to the supplementary methods in the paper you cited:

To control for this effect we performed hierarchical clustering of the RNA expression across all samples based only on the HOX genes (identified as genes with names starting with “HOX”) which revealed four major sample clusters

It sounds like they did not assign specific biopsy sites so much as cluster samples based on HOX gene expression in order to define an additional factor (i.e. the HOX expression cluster assignment of the sample) in their DE model design.

There is a Biostars tutorial Unsupervised clustering on gene expression data that may be a good starting point for this.