Hi Gordon, I'm sorry for any confusion! I'm looking at an allopolyploid (which has undergone a recent whole-genome duplication event via inter-specific hybridization) and am particularly interested in comparing gene expression between each gene pair. These genes are organized in large, highly similar syntenic regions throughout each ancestral genome (subgenomes) such that chromosomes copies have retained a fairly consistent gene order. For the sake of this analysis, I'm looking only at paralogous gene pairs which are single-copy genes in their respective subgenome (this constitutes the majority of genes).

Replicates 1, 2, and 3 refer to the biological replicates (each represents RNA extracted from a standardized tissue volume harvested from 10 individuals).

The duplicate genes (homeologs, in this case) aren't identical and have a number of polymorphisms that can aid in distinguishing subgenome-specific reads during read mapping (this is part of the motivation for gene-length normalization, although the vast majority are nearly identical in sum exon length). Reads which are best hits with a particular gene are kept, and the minority of identically (or somewhat ambiguously) mapped reads are subsequently discarded. There are other post-mapping tools which can do more comprehensive read-sorting but aligning alone is usually sufficient given the correlation between read-sorting outputs and aligners like STAR.

Given these counts I'd like to compare gene expression between subgenomes, but comparing subgenome expression using default edgeR/DESeq2 (which I believe include library corrections by default) would erroneously underestimate the difference in subgenome-mapped read counts since the reads originated from the same sample, and therefore, differences in mapped totals between subgenomes are biologically informative.

I hope that was less confusing but I'll admit this is a pretty non-standard use case for either of these tools.

RE my library size query: I'm just looking to see if manually including a library size (sum of A subgenome genes + sum of B subgenome genes) can be done in either program as the internal library correction would be problematic.

Edit: Here's an example of what I mean in edgeR:

library_size = function(x){
  # Library size for polyploid comparison:
  samples = gsub("\\..*$", "", colnames(x))  # remove subgenome IDs from colnames
  lapply(samples, function(sample){
    lib = x[,grep(sample, colnames(x))]
    sum(lib[,1],lib[,2])
  }) %>% unlist
}

lib_sizes = library_size(counts_merged)
DGE = DGEList(counts=counts_merged, group = group, remove.zeros = T, lib.size = lib_sizes)

Sorry, by RPK/FPK I'm referring normalization of my raw counts by gene length in kilobases (raw counts/(gene length/1000)). Since edegR and DESeq2 internally correct for library size I figured representing my data by FPK/RPK would be acceptable but a colleague brought up the fact that I shouldn't use this method as it would decrease the power.