Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-07-06T15:16:21.690Z
Answers.Text: 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?
Answers.UserDisplayName: LinuxBlanket
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-07-06T23:30:23.237Z
Answers.Text: 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%.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-07-07T00:27:32.747Z
Answers.Text: 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?
Answers.UserDisplayName: LinuxBlanket
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-07-07T00:29:27.827Z
Answers.Text: 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.
Answers.UserDisplayName: LinuxBlanket
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-07-07T04:54:00.060Z
Answers.Text: 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.
Answers.UserDisplayName: gringer
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-07-07T08:00:17.130Z
Answers.Text: 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 :)
Answers.UserDisplayName: LinuxBlanket
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-06-29T00:03:45.953Z
Answers.LastActivityDate: 2018-06-29T00:03:45.953Z
Answers.UserDisplayName: gringer
Answers.Body: <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>
  
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-06-28T16:27:27.713Z
Answers.Text: It's a clam. Aren't these values OK for a transcriptome? Where can I find good reference values?
Answers.UserDisplayName: LinuxBlanket
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-06-29T14:10:57.830Z
Answers.Text: Higher eukaryote? Talking about giraffes? http://scienceblogs.com/evolgen/2006/10/18/do-higher-eukaryotes-smoke-mor/
Answers.UserDisplayName: Kamil S Jaron
Answers.ContentLicense: CC BY-SA 4.0
Answers.CreationDate: 2018-06-28T16:25:51.520Z
Answers.LastActivityDate: 2018-06-28T16:25:51.520Z
Answers.UserDisplayName: Ian Sudbery
Body: <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>
  
Tags: <assembly><transcriptome>
Title: Reduce number of transcripts in a highly variable de novo transcriptome assembly
ThreadText: # 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](https://github.com/weizhongli/cdhit/wiki) 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:
  
  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>
  
  > **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/  
  
  ---
  
  