Differential expression analysis using two different de novo assembled transcriptomes

[Birdman]

Hi,

I'm trying to do DE analysis using DEseq (or DEseq2) to find species differences between RNA-seq samples (n=8) of two very closely related bird species.

In summary, I assembled two 'reference' transcriptomes out of all the Illumina reads (100bp PE) from every species using Trinity. I then blasted the two transcriptomes to obtain gene names for the assembled contig. Finally, I aligned the samples individually with their corresponding transcriptomes using BWA and obtained raw counts with eXpress. I'm now analyzing expression in R with DeSEQ.

The problem is that I sometimes have very similar sequences between both species that are have different gene names (according to a uniprot_sprot + nr blast, evalue < 1e-5, 1 result per query - I then blasted transcriptomes together to see if genes were corresponding to each other). Interestingly, it's usually only a subset of sequences (within the same "compxxxxx") that give different gene names, the other sequences have the same names. I'm not sure why this happens (alternative splicing/sequencing errors/other?) but it surely is a problem when I try to analyze gene expression. Some genes that appear to be 'overexpressed' in one species, when blasted against the other species' transcriptome, are very similar to other genes that, not surprisingly, are 'overexpressed' in the other species. It is most likely that these are in fact the same genes that are similarly expressed in both species.

I was thinking to remove all those ambiguous contigs, but it represents a lot of them (~50% of all contigs). Another idea I had was to blast a transcriptome with the other to obtain gene names instead of blasting against whole databases for every transcriptome.

What do you think would be the best approach to deal with this?

[TiborNagy]

You can not compare the names of two different Trinity output, because it is a de-novo tool. You can use Pasa to remove some of the duplications or use Trinotate to keep the best candidate ORFs.

[rskr]

I just want to emphasize that blasting transcriptomes against each other is especially naive, since it is fundamentally a local alignment, and any results generated thusly would be highly suspect. BLAST is fine for manually blasting a few transcripts then looking at the results, or for input into a more sophisticated transcript modeling program. However, it isn't going to chain transcripts together when you have a decent spanning match, compared to some small exact match somewhere else.

[Birdman]

Thank you both for your input. I will try using Trinotate's candidate ORFs.

[Wallysb01]

It probably is alternative splicing that that is shifting some isoforms to find better matches in other reference species or to other genes in the same family. Trinity may be assembling transcripts that are currently not annotated in some or all reference species in uniprot or nr. To make things easier for you I think you should make your blast searches more stringent and limit them to one reference species. If this is in birds, you might as well go with chicken refseq proteins, or something similar. Or if you’re following the trinioate package, you could use chicken uniprot.

If you keep following the trinity package for DE testing, it will give gene level tests that are basically pooling all your isoforms (in those compXXXXX) into a single gene. From there you could just take the annotation assignment from a single compXXXX_seqYY for each gene to compare between species. How you do this could be tricky, but there might be clever options like taking the highest blast score within the entire compXXXX set, or something similar.... Some of your own bash/python scripting would be needed for that. But getting your DE tests out of all those trinity isoforms and into the gene level instead is probably preferable.

This is not a simple thing you’re trying to do, so any route you take will have caveats.