Hello,

I have 4 independent species pairs (i.e., 4 x Species A vs. Species B contrasts). For each species, I have 4 biological replicates with RNA-seq data.

I am filtering lowly expressed genes, normalizing with edgeR, and using a modeling approach to identify genes with a convergent shift in expression from Species A to Species B across all 4 species pairs. This gives me an observed count of convergent-changing genes. I can provide specific details on the modeling framework if needed.

My objective is to use a permutation/randomization approach to construct a null distribution to test if the number of convergent genes I found is different than what would be expected by chance. I am considering two permutation strategies but they both seem to have some limitations.


Approach 1: Per-gene permutation where the gene read counts are randomly flipped or kept the same between the replicates of Species A and replicates of Species B within a species pair. This approach maintains intact the replicate structure but because it is done on a per-gene basis it breaks the multivariate covariance structure and sample-specific variance across the transcriptome.


Approach 2: Per-sample permutation where each individual replicate’s transcriptome is maintained intact but the group labels (“Species A” or “Species B”) are randomly permuted across replicates within each species pair (maintaining a balanced design of 4 replicates per group). This, however, breaks the replicate structure within each group.


I would run the Permutation thousands of times on the raw count matrix and then go through the filtering, normalization, and modeling to obtain a distribution of expected convergent changing genes to compare with the observed number.

I would appreciate any thoughts on these approaches or any other suggestions would be very welcome!