This is a test version of Biostars. For the public version, visit https://www.biostars.org.
Have bbduk retain quality scores in output after kmer trimming

I am trimming primers off of some .fastq files before processing them through the dada2 pipeline. dada2 wants primers already removed, but also wants quality scores. To trim primers I am using kmer trimming via bbduk and the following command:

bbmap/bbduk.sh in=input.fastq out=output.fastq literal=GTGYCAGCMGCCGCGGTAA ktrim=l k=10

The head of the input file looks like this:

@ERR1873185.1 AV103___MISEQ:377:000000000-ABUYP:1:1101:15575:1896/1
GTGTCAGCAGCCGCGGTAATACGAAGGGGGCTAGCGTTGTTCGGAATTACTGGGCGTAAAGCGCGTGTAGGCGGGTCTTTTATTCAGGGGGGAAATGCCCAGGCTCAACCTTGGAACTGCCTTTGATACTGGAGATCTTGAGTCCGGGGGGGGTGAGTGGATTTTCGAGTGTAGGGGTGAAATTCGTAGATATTCGCAAGAAAACCAGTGGCGAAGGCGGCTCCCTGGCCCGGGACTGACGCTGGGACGCGAAAGCGGGGGGGGCAACGAGGCTTAGCACCCCCTGTGG
+
4FFGGGFFAFGH?EE2BEGHHGHG2EGGGG?EGHF?FA0FHFGGGCFE4FFGF1?/EEDHHHDE/>/EFD3B@////0<11111<<1?F@@--;A9/B99:AE?FFEBF.//9///BF//99B//9999;//9;/9/:B////;FD?B----@-9../B.B/;B/;9.:/9//;/;..-;BB/;;.99..//99/.9..-.//.;.../:;.-;-9..@@>-9.;./.;A9B---;A./9....-;A.9@=;--./.-;==----;A.-..--.//////;:---.//.
@ERR1873185.2 AV103___MISEQ:377:000000000-ABUYP:1:1101:17496:1920/1
GTGCCAGCCGCCGCGGTAAGACGAACCGTGCGAACGTTGTTCGGAATCACTGGGCTTAAAGGGCCCGTAGGCGGGCTGTCAAGTCTGGGGTGAAATCCCGCGGCTCAACCGTGGAACTGCCTTAGATACTGACGGCCTCGAGGGAGGTAGGGGCGAGCGGAACTGTGAGTGGAGCGGTGAAATGCGTTGATATTCACAGGAACTCCGGTGGCGAAGGCGGCTCGCTGGCCCTCTTCTGACGCTGATGCGCGAAAGCTAGGGGAGCAAACGGGATTAGATACCCTGGTAG
+
EG5FGG4FFEEG2E?2E?FDBFEEGEGEFHEEE0EEHF1FGHFCEEFHGHH4B1FEFBBGGH?F/>/EEAEBECGGGGHHBBBGHHFFHCGCHHHHHHHFGG?D@FHGHGBAEGGGGGGGGGG0BCBFGFGBFD@BBFFFAA->.A-.//.AD?BB-;??9BBFBBB/;FFFFFFB<.;9FFF/BD..99BFFFFBFBBEFFFFFFFBBB.9;>B??A>;>9A--..-.AEFFFFFFF/B?DFD;FFF;.9-@-ABF/F???A-:;FFFF...;B//;//9F99BEFF:

The head of the output file looks like this:

>ERR1873185.3001 AV103___MISEQ:377:000000000-ABUYP:1:1101:6839:18042/1
TACGAAGGGGGCTAGCGTTGCTCGGAATCACTGGGCGTAAAGCGCACGTAGGCGGCTTTTTAAGTCAGGG
GTGAAATCCTGGAGCTCAACTCCAGAACTGCCTTTGATACTGGGAAGCTTGAGTCCGGGAGAGGTGAGTG
GAACTGCGAGTGTAGAGGTGAAATTCGTAGATATTCGCAAGAACACCAGTGGCGAAGGCGGCTCACTGGC
CCGGTACTGACGCTGAGGTGCGAAAGCGTGGGGAGCAAACAGGATTAGAAACCCTAGTAG
>ERR1873185.3002 AV103___MISEQ:377:000000000-ABUYP:1:1101:13081:18045/1
CACGTAGGGCGCAAGCGTTGTCCGGAATTATTGGGAGTAAAGAGCTCGTAGGCGGTCCGTCACGTCTGTT
GTGAAAATCCAGGGCTCAACCCTGGACTTGCTGCGGATACGGGCGGACTAGAGGTAGGTAGGGGAGAATG
GAATTCCCGGTGTAGCGGTGAAATGCGCAGATATCGGGAGGAACACCGGTGGCGAAGGCGGTTCTCTGGG
CCTTACCTGCCGCTGAGGAGCGAAAGCGTGGGGAGCGAACAGGATTAGAAACCCGTGTAG

how can I tell bbduk to just trim the primer sequence and retain the quality scores? Furthermore, The output file sequences are in a different order, anyway to have them stay in the same order? I'm worried my current bbduk call is also quality filtering/discarding which I do not want. Its difficult to tell from the file line counts because the sequences in the output of bbduk take up multiple lines in the output, rather than a single line.

trim primer bbmap
• 3,183 views • link
written by caverill • caverill

1 answer

That does not look correct.

BBTools will reformat output sequences/adjust input parameters based on file names (.fa is interpreted as fasta format needed/provided). It appears that your trimmed output is getting converted to fasta format even though the command you show above is for fastq format output. Can you address this discrepancy?

Note: You will need to include N's in literal= option (instead of IUPAC code) to cover all different possibilities.
Note 2: Add ordered=t to keep the output reads in same order as input.

• 1 views • link
written by GenoMax GenoMax

It is up some items: Usage Guide.

File Formats and Extensions

BBTools support most standard sequence formats, including fastq, fasta, fasta+qual, scarf, sam, and (if samtools is installed) bam. They also support gzip and (if bzip2 or pbzip2 is installed) bzip2. The tools are sensitive to file extensions.

• 1 views • link
written by h.mon h.mon

Also- bbduk can handle other degenerate base names than N. From bbduk documentation:

Matching degenerate sequences such as primers: bbduk.sh in=reads.fq out=matching.fq literal=ACGTTNNNNNGTC copyundefined k=13 mm=f

This will clone the reference sequences to represent every possibility for the degenerate bases (Ns and other non-ACGT IUPAC symbols).

• 1 views • link
written by caverill • caverill

One last question- the number of lines in the output file is fewer than the number of lines in the input file, implying some sequences have been removed. Specifically what are the criteria for removing a sequence from the output? As far as I can understand from the documentation this command trims primers, I have no idea under what conditions it removes sequences- if the k-mer sequence is not found within a read? something else?

• 1 views • link
written by caverill • caverill

If you had a primer dimer in the read then that would be almost certainly be thrown out. Other possibility is minlen=10 by default. If your read becomes shorter 10 bp after it gets trimmed then it will be removed. If a k-mer is not found in the read nothing should happen to it.

• 1 views • link
written by GenoMax GenoMax

Log in to answer this question.

More posts like this
  • Paired-end flag for singleton
    written by docdot • docdot
    Dear all, I'm having some trouble to identify singletons in paired-end sequencing data from Hi-C. I have a Hi-C library originated from 150 bp (75x2) …
  • Trouble de-multiplexing single index paired-end illumina library
    written by caverill • caverill
    I recently received a sequencing run from a sequencing center. This is an Illumina paired end Miseq run, Single-index. I usually get per-sample files post …
  • Sorting fq data into groups with missing barcode info?
    written by jorja • jorja
    10/17/18 clean up file names and added reads per file Use "header" instead of reads. each Header means each @M70287:XXXX... Dear all, I am trying …
  • Parsing large sequence file into per-sample files after barcodes/primers removed.
    written by caverill • caverill
    I have a sequence file, called `seqs.fna`. The sequences have had sample names assigned, barcodes and forward primers removed, and have been quality filtered. These …
  • Extract Reverse Reads Based on Header Information of Forward Reads
    written by dennis.guettler • dennis.guettler
    Hi Community. **Context:** I used MiSeq for sequencing a metagenome with paired end setup. But due to overclustering, my reverse reads can't be matched by …
  • Generate an amplicon sequence variant table without quality data using dada2 or similar in R
    written by caverill • caverill
    There are sequences I download from MG-RAST that already have primers removed, quality filtering performed, and forward-reverse reads paired. Each file is a fasta file …
  • Paired-read alignment length from BAM/SAM file
    written by Soumitra Pal • Soumitra Pal
    Hi, I have aligned a set of paired reads using BWA. I want to know the alignment length of each of the pairs. By length …
  • Empty output trimmomatic
    written by luigi.marongiu • luigi.marongiu
    <p>Dear all,<br /> I am applying trimmomatric to trim fastaq files by quality and to remove the adapters. I have two paired files seq1.1.fq and …
  • fastq to fasta perl
    written by cabraham03 • cabraham03
    Hi, I have the next code to convert from fastq to fasta, however I have some problems to be fix. I want from: ``` @YOX24:00004:00021 …
  • searching reads with a certain sequence in fastq file
    written by sumithrasank75 sumithrasank75
    I have a file say master.fastq which looks like : ``` @M00990:202:000000000-ADM27:1:1101:21678:1536 2:N:0:291 CCTTTTACCGACCCGCTCTTTCTCTCCTACGCTTATTTCCGTCTACCCTTCTCTTCACTCGCTATTTCTATTCTTAAAACTATCTTAATGTTCTGCCTTTGCTCTTTTCTTTTTTCTATAACCTCTCTACAGCCAACTCACCCATCTCCTTCCCTGCTACGCTATTCCTCTGTTAGTTTTTCTTCATCATACTTTTCTCATCTCACACTACCTTTGCACTTCTTCCTTTCCACGTCCCCTTTCTCCTACC + -----,,6;6,@+8++6+,,,<C5@,,,,,,+8,,;,,<,,,,,,,;,,;,,,<C,,,;,,8,,8CC,C,,,5,C,,,,99,+,+4,,,3,,9,,6,@<,,,,,,,9,,,,,,,4),,0**,,5),,5**,)59*0),,*)5,)))*,9,3,0++))5D:)))))5;+;+0)*;+*6++++),******3**50,6,+++**,0*,,31)88*0*1*5*1)0*:*7>C;3,035:0)))8).*2**.*:) @M00990:202:000000000-ADM27:1:1101:22685:1539 2:N:0:291 CTTATCACCGACTCTCTCCTTCTCTTCCAAGTTTATTTCCGACTCCCCTTATCTTCACTTGCTATTTCTATTCTTAAAACTATCTCGACCTTTCACCTTTCCCTCTTTCCTTCTTTTCTCTCCTTCTACACTCCCACCCACTCTTACTTCTTTCTTGTCACCGTTTCCATATTATACTTTCTTCTCTTACATAATTTTCTTCCTGCAAACTATTTAAGCAATCTCTTTCTTTCACCCCTTTTATCTCGCC + -----,,<;67@+B,,6,,,,;C5@,,,,,66<,,6,,<,,,66,,;,66,,;;C,,,;,4<,,<CC6C,;,,,;,,6,:,6?A9=,,+++2,5,,,9<?D,,,,,:C?,@,@,,5,?*3*,9**,,0**,*)93))4+))19*0**,,52,56+**5*03*)3)))42+2***5*+=3+,*4*2****,**2*,3,+++*0,50,,**5*****0****)5)***0**,,***3*)0)))3***0*))) ``` I want only …