[duane.hassane]

Thanks. To be more clear, I am after removing the primers themselves (either from the fastq or BAM) and preserving the intervening sequence only. In this use case, there are hundreds of different primer pairs with known mapping locations (because they are used for targeted enrichment of different exons). Thus, for each primer pair, the intervening sequence (target being enriched) is different.

[GenoMax]

I think @duane.hassane is asking to remove primers from amplicons with a common core sequence but with variable (length and/or sequence?) primers on the ends of that core sequence for each group.

[Brian Bushnell]

BBTools has is a pair of programs for cutting (and outputting to a file) the sequence between primer pairs, but not for cutting out primers themselves. Is that what you are after?

The usage is like this:

Assume your set of primers for one end is AAAAAA and AAAAAT, and for the other end is GGGGGG and GGGCCC:

msa.sh in=amplicons.fq out=primer1.sam literal=AAAAAA,AAAAAT
msa.sh in=amplicons.fq out=primer2.sam literal=GGGGGG,GGGCCC

That will find the optimal alignment for the optimal primer, and output one line per amplicon (twice). Then run this:

cutprimers.sh in=amplicons.fq out=middle.fq sam1=primer1.sam sam2=primer2.sam

"middle.fq" will contain the sequence between the primers, one per amplicon, using the best alignment. I designed this to cut V4 regions from full-length 16s.

Notes: msa.sh stands for MultiStateAligner, not for Multiple Sequence Alignment. All sequences can be in any orientation (forward or reverse complement).

Please let me know if that doesn't address your question; I may have misunderstood it.

[duane.hassane]

Is there a way of efficiently using the BBMap tools to cut variable length and variable sequence primer pairs from amplicon sequencing data either pre-alignment (fastq) or post-alignment (BAM)? I am thinking in regards to RainDance and Fluidigm data where there are hundreds of amplicons, many of which overlap, each with different forward/reverse primers.

[Brian Bushnell]

Those flags work in BBDuk as well, except for "untrim", which is unique to BBMap.

A good variant caller will take the quality into account when weighing possible variations. Quality-trimming before mapping gives the mapper less information, and thus a higher probability of ambiguously-mapped or incorrectly-mapped reads (though adapter-trimming before mapping is always good). Even if the trimmed bases are low quality - say, 50 bases at Q13 - that's still an expected 45 matches and 5 mismatches, which can greatly increase mapping confidence in a repetitive area, and thus improve the ability to call the variations from the high-quality part of the read.

On the other hand, quality-trimming before assembly is often a good idea, though the exact threshold depends on the assembler and dataset.

[fdts]

Aah, I see. Thanks to both of you. Would these flags work with BBDuk then?

Data quality is decent. I will map to a de novo metagenome assembly generated from the same data. I am kind of hopeful to be able to do a little variant calling after metagenome binning and mapping. Given this situation, would you still rather not quality trim?

Thanks

[Brian Bushnell]

To expand upon that, BBMap has the options for internally quality-trimming (and untrimming) reads, but I only use them in special cases. For example, BBMap is used a lot for filtering out contamination. To ensure high specificity, I only want to remove reads that map to contaminant sequences with high identity; but I don't want the identity calculation to include low-quality bases. So, I map with "qtrim=rl trimq=15 untrim" or similar, which will trim before mapping but then output the untrimmed read.

For normal cases such as variant calling or coverage calculation, don't quality-trim unless you have really low quality data. I would not expect that to be the case for a 100bp library.

[GenoMax]

Do you recommend quality trimming before mapping using BBMap? I am talking illumina 100bp PE.
I have been searching the forum for a while but could not find a hint.

Thanks

Unless you have a lot of data that has Q scores of 10 or less you can get away without quality trimming.

[fdts]

quality filtering

Do you recommend quality trimming before mapping using BBMap? I am talking illumina 100bp PE.
I have been searching the forum for a while but could not find a hint.

Thanks

[Brian Bushnell]

I would also really like such a tool! Yes, I DO have plans to do something like that, but I am not sure about the timing, since I have had these plans for almost a year now. But hopefully I will make it this year, as Dedupe is really close to being an assembler, just not quite there.

[Zapages]

Hi Brian,

Thank you for creating these awesome applications. I was wondering if you have any plans to create an OLC aligner or something that improves Dedupe so that it can merge multiple assemblies together while taking care of any duplicate sequences based on a certain or user specified percentage. Then finally it outputs a consensus sequence based on a user specified percent similarity and some other variables.

I believe this will be very helpful with working with viral or bacterial de-novo genome (DNA-Based) Assemblies. Some ideas.

[sdriscoll]

Excellent. My PI will be stoked.

[Brian Bushnell]

I was thinking about an automatic paper generator, which would be more useful to me. But I guess grant proposals are even more standardized. I'll see if I can roll that into Reformat.

[sdriscoll]

Hey Brian,
A feature many of the bench scientists I work with would love is an automated grant proposal generator. Any chance of something like that in future updates? The input could be some images/figures and a basic subject...

[Brian Bushnell]

Oh, that's good to know. My exposure to plant RNA-seq is mainly Chlamy and Arabidopsis - and you know what they say about assumptions.