It's this, my purpose is to find some genes having correlation to the gleason scores in prostate cancer, especially the positive correlation between the grade of gleason scores and gene expression values. I download TCGA prostate cancer data through tcgabiolinks package in R, and I retain only tumor samples. I created a group variables via gleason scores. first grade for <=6, second grade for 3+4, third grade for 4+3, fourth grade for 4+4, 3+5, 5+3, fifth grade for 9~10, the code is below:

gleason_group <- ifelse(gleason_data$subtype_Clinical_Gleason == "3+3"| gleason_data$subtype_Clinical_Gleason == "2+4",
                    "I", 
                    ifelse(gleason_data$subtype_Clinical_Gleason == "3+4", 
                           "II", 
                           ifelse(gleason_data$subtype_Clinical_Gleason == "4+3", 
                                  "III", 
                                  ifelse(gleason_data$subtype_Clinical_Gleason == "4+4"| gleason_data$subtype_Clinical_Gleason == "3+5" | gleason_data$subtype_Clinical_Gleason == "5+3", 
                                         "IV", 
                                         "V"))))

and for the dependant variable is a ordinal discrete data, so I think the ordinal logit model can help me. I'm not familiar with this, Am I right? and I find some packages, example for "ordinal" and "rms". which is the best suitable for my conditions? Or I shouldn't use logit model and there are a better solution to me?

I appreciate you in advance.