We have RNA-seq from two species.

We used Kallisto to estimate the gene expression by mapping reads to corresponding references.

After that we use the tximport to get gene level read counts and gene length.

Now we are trying to use DESeq2 to do differential expression analysis and we want to include the gene length for correction between species.

But not sure whether our following codes are correct or not. I know there is "DESeqDataSetFromTximport" function, but since our analysis is between species, we are not sure how to set "tx2gene" to two different species.

Below is briefly our codes. Is this correct and DESeq2 will use the gene length when do differential analysis?

construct the dds_obj: readCounts are at gene level from tximport;

dds_obj <- DESeqDataSetFromMatrix(countData = readCounts,colData = colData,design =~Species+Sex)

add the gene length matrix from tximport corresponding to readCounts matrix

assays(dds_obj)[["avgTxLength"]] <- geneLength

dds <- DESeq(dds_obj)

Thanks

When you quantify, is the reference the same length for both species? What exact steps have you done to quantify?

If you were using DESeqDataSetFromTximport, it would automatically adjust for the differences in gene length across samples, and you wouldn't have to attempt to code it manually. I'd suggest to use that therefore.