Skip to content

Removing adapters from Trinity Transcriptome assemblies

Jump to bottom
mathog edited this page Jan 25, 2019 · 12 revisions

When Trinity 2.8.3 was run with or without --trimmomatic on SRR531964 there were Illumina adapters present in the resulting transcript file. This was discovered when running the recommended NCBI UniVec screen for a TSA submission, like so:

blastn -task blastn -reward 1 -penalty -5 \
    -gapopen 3 -gapextend 3 -dust yes -soft_masking true\
    -evalue 700 -searchsp 1750000000000 -outfmt 6 \
    -db $PATH_TO_UNIVEC/UniVec \
    -query 283_trinity_SRR531964.Trinity.fasta \
    -num_threads 40 -out /tmp/283_trinity_SRR531964.fmt6.out

blastn -task blastn -reward 1 -penalty -5 \
    -gapopen 3 -gapextend 3 -dust yes -soft_masking true \
    -evalue 700 -searchsp 1750000000000 -outfmt 6 \
    -db $PATH_TO_UNIVEC/UniVec \
    -query 283_trinity_trim_SRR531964.Trinity.fasta \
    -num_threads 40 -out /tmp/283_trinity_trim_SRR531964.fmt6.out

grep NGB0036[02] /tmp/283_trinity_SRR531964.fmt6.out | head -5
TRINITY_DN1596_c0_g2_i6 uv:NGB00360.1:1-58      100.000 37      0       0       1440    1476    58      22      6.03e-11        74.6
TRINITY_DN1544_c0_g1_i10        uv:NGB00360.1:1-58      100.000 37      0       0       1       37      22      58      6.03e-11        74.6
TRINITY_DN1503_c0_g1_i6 uv:NGB00362.1:1-61      100.000 24      0       0       1       24      38      61      0.004   48.5
TRINITY_DN1503_c0_g1_i8 uv:NGB00362.1:1-61      100.000 24      0       0       1       24      38      61      0.004   48.5
TRINITY_DN1532_c0_g1_i5 uv:NGB00360.1:1-58      100.000 38      0       0       723     760     58      21      1.50e-11        76.6

grep NGB0036[02] /tmp/283_trinity_trim_SRR531964.fmt6.out | head -5
TRINITY_DN1596_c0_g1_i1 uv:NGB00360.1:1-58      100.000 23      0       0       2684    2706    58      36      0.017   46.5
TRINITY_DN1559_c0_g1_i4 uv:NGB00362.1:1-61      97.436  39      1       0       2060    2098    61      23      1.57e-08        66.6
TRINITY_DN1516_c0_g1_i1 uv:NGB00362.1:1-61      100.000 39      0       0       3       41      23      61      3.74e-12        78.6
TRINITY_DN1516_c0_g1_i2 uv:NGB00362.1:1-61      100.000 39      0       0       3       41      23      61      3.74e-12        78.6
TRINITY_DN1516_c0_g1_i3 uv:NGB00362.1:1-61      100.000 39      0       0       3       41      23      61      3.74e-12        78.6

Adapters NGB00360 and NGB00362 are both found with perfect (or nearly so) matches extending to their 3' ends. There are many other matches which are closely related to these adapters, but these are the longest matches. When --trimmomatic was employed the log file contained:

Using PrefixPair: 'TACACTCTTTCCCTACACGACGCTCTTCCGATCT' and 'GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT'

The adapters found in the assembly by blastn are these two Univec entries:

>gnl|uv|NGB00360.1:1-58 Illumina PCR Primer (Oligonucleotide sequence copyright 2007-2012 Illumina, Inc. All rights reserved.)
AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT

>gnl|uv|NGB00362.1:1-61 Illumina Paired End PCR Primer 2.0 (Oligonucleotide sequence copyright 2007-2012 Illumina, Inc. All rights reserved.)
CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCT

The first trimmomatic sequence matches the 3' end of NGB00360 while the second Trimmomatic sequence differs a lot from NGB00362.

                        TACACTCTTTCCCTACACGACGCTCTTCCGATCT  Trimmomatic sequence 1
AATGATACGGCGACCACCGAGATCTACACTCTTTCCCTACACGACGCTCTTCCGATCT  NGB00360 

                           GTGACTGGAGTTCAGACGTGTGCTCTTCCGATCT Trimmomatic sequence 2
CAAGCAGAAGACGGCATACGAGATCGGTCTCGGCATTCCTGCTGAACCGCTCTTCCGATCT NGB00362

A small script "trim_adapters.pl" was written to identify and remove known adapters. The way it works is that first the blastn screen against UniVec is run. Identify any adapters which are present and then remove them like so:

./trim_adapters.pl  \
   -infile testing.fasta \
   -outfile testing_out.fasta \
   -scanfile univec_vs_transcriptome_reduced.fmt6.out \
   -adapters "NGB00360:58 NGB00362:61" \
   -pepfile testing.pep \
   -minlen 200 \
   -relaxed \
   -drop_nopep

(To reduce the run time restrict the search space. If the adapters are already known put those in a fasta file "subset". If the adapters are not known first do an initial search of Univec with a small subset of the mRNA sequences. Identify the adapters in those results. Then extract those Univec adapter entries into the "subset" file. Create a blast database from it like so:

   makeblastdb -in subset -dbtype nucl -parse_seqids

Lastly, use the "subset" database instead of "UniVec" in the blastn search.)

Use the -help option to see the parameter descriptions and other information. It runs reasonably quickly, even though it is single threaded, and emits a series of messages like this:

2018-10-04 14:07:22 Starting
Command line arguments were infile: testing.fasta outfile: testing_out.fasta scanfile: univec_vs_transcriptome_reduced.fmt6.out adapters: NGB00360:58 NGB00362:61 pepfi
le: testing.pep minlen: 200 drop_nopep: yes
2018-10-04 14:07:22 Obtaining sequence lengths from: testing.fasta
2018-10-04 14:07:22 Processing pepfile: testing.pep
2018-10-04 14:07:24 CDS ranges:110582, 2nd or higher CDS range for transcript:0
2018-10-04 14:07:24 Processing scanfile: univec_vs_transcriptome_reduced.fmt6.out
2018-10-04 14:07:24 Processing testing.fasta
2018-10-04 14:07:26 Done  Read:110582 Wrote:110578 Trimmed:9914 newlen<200:0 deltaLen>100:56 cds_cut:6939 cds_drop:1 len_drop:3 noprot_drop:0

In this instance 1 transcript was dropped because the entire CDS was on the wrong side of the adapter (cds_drop), 3 transcripts were dropped because they were less than 200bp after trimming, none were dropped because there was no corresponding protein for the transcript (ie, here they were 1:1). Of the remaining 110578 emitted transcripts 9914 where trimmed, none were <200bp (because those were dropped) and 56 had their length change by more than 100bp.

To see what remains of the trimmed adapters, run the blastn screen again and then check the status of the specified adapters:

blastn -task blastn -reward 1 -penalty -5 \
    -gapopen 3 -gapextend 3 -dust yes -soft_masking true\
    -evalue 700 -searchsp 1750000000000 -outfmt 6 \
    -db $PATH_TO_UNIVEC/UniVec \
    -query testing_out.fasta \
    -num_threads 40 -out testing_out.fmt6.out

grep NGB0036[02] testing_out.fmt6.out
(1) SRR531843_TRINITY_DN16330_c0_g1_i4      uv:NGB00362.1:1-61      100.000 16      0       0       1724    1739    6       21      287     32.5
(2) SRR531843_TRINITY_DN5608_c0_g2_i1       uv:NGB00360.1:1-58      96.154  26      1       0       10      35      29      54      1.1     40.5
(1) SRR531949_TRINITY_DN13099_c0_g1_i1      uv:NGB00362.1:1-61      100.000 16      0       0       1340    1355    43      28      287     32.5
(1) SRR531954_TRINITY_DN2435_c0_g2_i5       uv:NGB00362.1:1-61      100.000 21      0       0       32      52      1       21      0.27    42.5
(3) SRR531955_TRINITY_DN33086_c0_g3_i1      uv:NGB00360.1:1-58      90.566  53      5       0       10      62      1       53      0.017   46.5
(2) SRR531955_TRINITY_DN5640_c0_g1_i5       uv:NGB00360.1:1-58      100.000 24      0       0       7       30      35      58      0.004   48.5
(2) SRR531955_TRINITY_DN5640_c0_g1_i6       uv:NGB00360.1:1-58      100.000 24      0       0       7       30      35      58      0.004   48.5
(1) SRR531956_TRINITY_DN23092_c0_g1_i1      uv:NGB00360.1:1-58      100.000 16      0       0       254     269     24      9       287     32.5

Where (1) are short sequences which are similar to the adapter, (2) are instances where the adapter was present twice, and (3) is ambiguous. The (2) issue arises when an adapter is present multiple times in the transcript like this:

   -------------------------->  transcript
   -A> -B>                      two copies of the adapter

The blastn command only returns a single alignment. If that is the A alignment the B adapter will still be present. If this happens run another cycle of this script. Note, that usually happens, expect to run two cycles. (The example above was from an earlier version of the script, current versions remove all of the type (3) contaminants in the first pass, and type (2) contaminants are rare.)

The usual recommendation is to trim the reads before assembly. Currently the built in --trimmomatic does not remove all adapters, perhaps the stand alone trimmomatic program would resolve most of these. In any case, whatever slips through to the final assembly may be addressed with this script.

The Transcriptome assembly may of course contain contaminants other than Illumina adapters. Those will also be listed in the output of the blastn screen against Univec, but this tool will not act on them.

Note added 01/22/2019. Updated script so that it handles different blast database formats better. Depending on how the database was built the name of the vector match will change format. Added -inranges -outranges so that on multiple cycles of running this program the part that remains is annotated while the original length is also noted.

Note added 1/24/2019. Added -relaxed switch, which should generally be used. Makes a more conservative cut to remove bases at the 3' end which might not have been included in the adapter alignment found by blastn, and trims a few more adapters which aligned but not quite as well. With -relaxed employed it was observed that in the second of two cycles with trim_adapters.pl on 23 libraries only a single mRNA family had an adapter contaminant. So the first cycle was sufficient for 22 libraries and nearly sufficient for 1 other.

Clone this wiki locally