5
$\begingroup$

I have a de novo assembly using both multiple SRA and locally sequenced transcriptomes. I started with 270M PE reads from 9 tissues. Here are the assembly stats generated with TrinityStats.pl:

################################
## Counts of transcripts, etc.
################################
Total trinity 'genes':  543668
Total trinity transcripts:      1111089
Percent GC: 33.77

########################################
Stats based on ALL transcript contigs:
########################################

        Contig N10: 2117
        Contig N20: 1324
        Contig N30: 945
        Contig N40: 714
        Contig N50: 561

        Median contig length: 344
        Average contig: 496.17
        Total assembled bases: 551293315


#####################################################
## Stats based on ONLY LONGEST ISOFORM per 'GENE':
#####################################################

        Contig N10: 2172
        Contig N20: 1333
        Contig N30: 925
        Contig N40: 681
        Contig N50: 523

        Median contig length: 323
        Average contig: 473.63
        Total assembled bases: 257499729

As you can see, the total Trinity 'genes' number is very high. I suspect that this is due to a high biological variability of the transcriptomes (coming from organisms sampled in Italy, China, Spain...).

  • How can I assess whether the problem is really the biological variability?
  • How can I reduce this number? I'd rather not downsample the reads. It looks like CD-Hit might do what I'm looking for...
Improve this question
$\endgroup$

2 Answers 2

Reset to default
3
$\begingroup$

As you've suggested, CD-HIT works for reducing transcript numbers. We used a mixture of expression-based filtering and CD-HIT for reducing transcript counts for our genome-guided transcriptome assembly. This reduced numbers by a lot, without much change in BUSCO scores:

  1. Map RNA-seq reads to Trinity-generated transcripts using Salmon
  2. Use the expression of BUSCO genes to set a credible signal cutoff
  3. Subset transcripts based on this threshold [for us it was 50 counts]
  4. Identify the longest Met to Stop open reading frame for each transcript (may not be a good idea in all cases)
  5. Run cdhit to collapse similar transcripts

Here's our CD-HIT command line:

cdhit -T 10 -c 0.98 -i longest_MetStopORF_HC50_TBNOCFED.fasta -o cdhit_0.98_LMOHC50_TBNOCFED.prot.fasta

More details here:

https://bmcbiol.biomedcentral.com/articles/10.1186/s12915-017-0473-4#Sec9

Improve this answer
$\endgroup$
6
  • $\begingroup$ Thank you for your answer. Prior to using CD-HIT, I decided to filter the reads according to BLAST hits with a custom pipeline. Only 150k transcripts had a hit, and of these, only 56K had a hit with the suitable taxon. This means that almost a million transcripts don't match with anything, and probably CD-HIT won't help much. Any idea on why this is happening? Also, given that BUSCO retrieved 975/978 core genes (so my assembly is basically complete), should I care about this number? $\endgroup$
    – LinuxBlanket
    Commented Jul 6, 2018 at 15:16
  • $\begingroup$ CD-HIT does an internal match; it doesn't care about matches to external databases. That BUSCO score is great; even C. elegans isn't 100%. $\endgroup$
    – gringer ♦
    Commented Jul 6, 2018 at 23:30
  • $\begingroup$ Yes, but my question is: should I worry about the fact that I have >1M transcripts, and of these only 5% seem to correspond to a mRNA? What's in the other 950k transcripts reconstructed by trinity? $\endgroup$
    – LinuxBlanket
    Commented Jul 7, 2018 at 0:27
  • $\begingroup$ Ah, and just for the record, today I ran CDHIT and it obtained 780k clusters, so yes, it improved the situation, but not to a great extent. This still looks unrealistic to me. $\endgroup$
    – LinuxBlanket
    Commented Jul 7, 2018 at 0:29
  • $\begingroup$ To better answer that question, I need more information about what you were sequencing, how the sequencing was done, and what sequences are coming up. For example, if there's a high chance of bacterial contamination and sequencing is for total RNA including microRNAs, then 780k is something that could happen. $\endgroup$
    – gringer ♦
    Commented Jul 7, 2018 at 4:54
0
$\begingroup$

Seems to me that your "genes" are very short. What sort organism (prokaryote, lower eukaryote, higher eukaryote) is this? I would look to see if your "genes" contained ORFs, and only keep those with sufficicntly long CDS. You will miss uORFs and ORF fragments this way, but you'd have much more confidence in what you did find.

Improve this answer
$\endgroup$
2
  • $\begingroup$ It's a clam. Aren't these values OK for a transcriptome? Where can I find good reference values? $\endgroup$
    – LinuxBlanket
    Commented Jun 28, 2018 at 16:27
  • 1
    $\begingroup$ Higher eukaryote? Talking about giraffes? scienceblogs.com/evolgen/2006/10/18/… $\endgroup$
    – Kamil S Jaron
    Commented Jun 29, 2018 at 14:10

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.