In the past I've run separate multiple regression models for many correlated independent variables and one dependent variable. For this I've been using the R package multtest (http://www.bioconductor.org/packages/release/bioc/html/multtest.html). This allowed me to compute adjusted p-values that took the correlation structure of my matrix of independent variables into account.
Now, I want to do the same thing but for several dependent variables as well. Put differently, I have a matrix Y with my dependent variables and a matrix X with my independent variables. From this I want to estimate X*Y regression models. Importantly, I want the adjusted p-values to account for the correlation structure of both the independent and the dependent variables. I'm looking forward to your suggestions. I've suggested to the authors of multtest before to extend their library to accomodate this case but this hasn't happened yet.
Example (added): Let's say I have gene expression data from 10 different tissues. Now I want to know if gene expression is correlated with a 100 different SNPs. This means I'm effectively testing 100*10 = 1000 hypothesis. However, all these hypothesis are not independent of each other. The SNPs might be correlated to each other due to linkage disequilibrium and gene expression might also be correlated accross different tissues, depending on their similarity. Therefore a Bonferroni correction of my p-values for this 1000 statistical tests would be too conservative. I'm looking for a way to derive adjusted p-values that accounts for the above described dependencies within both the independent and the dependent variables.