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
Thanks for your reply. "is the reference the same length for both species?" Do you mean genes? No, since they are from different species, they have different reference transcripts (of course with different length). So for each species I ran Kallisto against their own reference transcripts. Then I get the gene-level readCounts and Length matrix for both species through tximport. Now I am trying to do differential analysis between orthologs between species using DESeq2 and I want to make use of the Length matrix along with the read counts.
If you can find a way to import the quants altogether, you could straightaway use tximport and DESeqDataSetFromTximport which controls for gene length across samples automatically.
I want to try your suggestion, but for the two species I have two different "tx2gene" (i.e., mapping file from transcripts to genes), can I just concatenate the two file into one? But in this case, for a gene each transcript will only express in one species, which is different from analysis between conditions for the same species. I am not sure how can I make use of DESeqDataSetFromTximport for different species. Many thanks
I see, so you can create two txi lists, and then manually combine each of the matrices, so you get a new list of matrices which are wider than previously. We don’t have any code for this but you can probably hack it with some base R.
Thanks. I will try it.