I am a complete noob to RNA-seq, but I'm trying to use DESeq2 to do differential expression analysis in a rather unusual system. I would like to compare gene expression between two closely related species (same genus). I have reference genomes for both. The RNA is just from control conditions, no treatment. I'm aware that a common approach for cross-species RNA-seq is to analyze only conserved exons, but because these species are so closely related--most genes in species A have a 1:1 ortholog in species B--that I'm hoping to just use a gene-level approach. Currently, I'm analyzing the data in two ways:

(1) Align reads from both species to the reference genome of species A, and proceed with downstream analysis as if my samples were two different treatments on species A.

(2) Align reads from each species to its own reference genome. Discard all genes that don't have a 1:1 orthology relationship. Change the gene names in the species B raw counts files to the names of their species A orthologs. Proceed with downstream analysis as if these were two different treatments on species A, analyzing only these "shared" genes.

I like this second approach because it takes advantage of both reference genomes (which my collaborators would like me to use), but I'm not sure if I'm making any truly bad assumptions, for example if it's bad that gene X is a different length in species A from its ortholog in species B. With either approach (1) or approach (2) I get a similar number of DEGs (a bit more for (1), as we would expect), and the majority of the DEGs from either analysis are called significant in both. Does anyone have opinions on this? Am I completely crazy for even attempting strategy (2)? Thank you so much in advance!!!

I don't have a solution per se, I think this kind of analysis can be difficult to remove all "technical" issues to get to cross-species gene expression changes.

(1) has a systematic bias, because you're favoring the reads from species A, due to alignment problems for reads from species B.

(2) sounds better. You can use the gene lengths as a normalization factor. Make a matrix which gives the gene length for each gene and each sample (so the same size as the count matrix). Then before you run DESeq(), you can store this matrix to:

assays(dds)[["avgTxLength"]] <- mat

This will be picked up by DESeq() and used for normalization.