I have a list of 250 genes whose expression was measured in 21 samples. I want to discover which gene pairs are positively correlated. So, I performed Pearson correlation for these genes using the "rcorr" function of the Hmisc package (R). The output gives me the r and p-values. What I want to know is: do I have to perform p-value correction, like Bonferroni or FDR? In the case of an affirmative answer, which R package would be the best option to perform this task?