<p>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:</p>

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

<p>Here's our CD-HIT command line:</p>

<pre><code>cdhit -T 10 -c 0.98 -i longest_MetStopORF_HC50_TBNOCFED.fasta -o cdhit_0.98_LMOHC50_TBNOCFED.prot.fasta
</code></pre>

<p>More details here:</p>

<p><a href="https://bmcbiol.biomedcentral.com/articles/10.1186/s12915-017-0473-4#Sec9" rel="nofollow noreferrer">https://bmcbiol.biomedcentral.com/articles/10.1186/s12915-017-0473-4#Sec9</a></p>

CC BY-SA 4.0
2018-07-06T15:16:21.690Z
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?
LinuxBlanket
CC BY-SA 4.0
2018-07-06T23:30:23.237Z
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%.
gringer
CC BY-SA 4.0
2018-07-07T00:27:32.747Z
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?
LinuxBlanket
CC BY-SA 4.0
2018-07-07T00:29:27.827Z
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.
LinuxBlanket
CC BY-SA 4.0
2018-07-07T04:54:00.060Z
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.
gringer
CC BY-SA 4.0
2018-07-07T08:00:17.130Z
Mind if we talk about it in [chat](https://chat.stackexchange.com/?tab=site&host=bioinformatics.stackexchange.com)? I'm afraid this is getting long for comments :)
LinuxBlanket
CC BY-SA 4.0
2018-06-29T00:03:45.953Z
2018-06-29T00:03:45.953Z
gringer
<p>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. </p>

CC BY-SA 4.0
2018-06-28T16:27:27.713Z
It's a clam. Aren't these values OK for a transcriptome? Where can I find good reference values?
LinuxBlanket
CC BY-SA 4.0
2018-06-29T14:10:57.830Z
Higher eukaryote? Talking about giraffes? http://scienceblogs.com/evolgen/2006/10/18/do-higher-eukaryotes-smoke-mor/
Kamil S Jaron
CC BY-SA 4.0
2018-06-28T16:25:51.520Z
2018-06-28T16:25:51.520Z
Ian Sudbery
<p>I have a <em>de novo</em> 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 <code>TrinityStats.pl</code>:</p>

<pre><code>################################
## 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
</code></pre>

<p>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...). </p>

<ul>
<li>How can I assess whether the problem is really the biological variability?</li>
<li>How can I reduce this number? I'd rather not downsample the reads. It looks like <a href="https://github.com/weizhongli/cdhit/wiki" rel="nofollow noreferrer">CD-Hit</a> might do what I'm looking for...</li>
</ul>

<assembly><transcriptome>
Reduce number of transcripts in a highly variable de novo transcriptome assembly
# Reduce number of transcripts in a highly variable de novo transcriptome assembly

(Asked by: username_0 on 2018-06-28)

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...

**Tags:** assembly, transcriptome

---

## 2 answers

### Answer by username_1 on 2018-06-29. Score: 3 | Accepted answer

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:

Map RNA-seq reads to Trinity-generated transcripts using Salmon
Use the expression of BUSCO genes to set a credible signal cutoff
Subset transcripts based on this threshold [for us it was 50 counts]
Identify the longest Met to Stop open reading frame for each transcript (may not be a good idea in all cases)
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

> **Comments:**
> **username_0**: 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?  
> **username_1**: 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%.  
> **username_0**: 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?  
> **username_0**: 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.  
> **username_1**: 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.  
> **username_0**: Mind if we talk about it in [chat](https://chat.stackexchange.com/?tab=site&host=bioinformatics.stackexchange.com)? I'm afraid this is getting long for comments :)  

---

### Answer by username_2 on 2018-06-28. Score: 0

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.

> **Comments:**
> **username_0**: It's a clam. Aren't these values OK for a transcriptome? Where can I find good reference values?  
> **username_3**: Higher eukaryote? Talking about giraffes? http://scienceblogs.com/evolgen/2006/10/18/do-higher-eukaryotes-smoke-mor/  

---

