[lac302]

Trying out bbmap for the first time...I'm aligning illumina 2x275 reads to a fungal genome. Error and sub rates for the first few samples are high, around 70%. I'm assuming the high error rate is due to the library prep (nextera xt). Insert size ranges from 400-2k. Any explanation for the high sub rate or are they both tied together. Thanks.

[Shini Sunagawa]

Just to make sure we are on the same page, I attach the toy input file I described earlier.

Running:

dedupe.sh in=in.fa out=out.fa outd=outd.fa am=t ac=t fo=f dot=dot.txt

will give the desired output:
seq1 and seq4 in out.fa
seq2 and seq3 in outd.fa (since they are contained in seq1)

but no dot.txt file.

Running:

dedupe.sh in=in.fa out=out.fa outd=outd.fa am=f ac=f fo=t dot=dot.txt

as you suggest, will output all 4 sequences in out.fa and none in outd.fa (again no dot.txt). To get the dot.txt file written, I have to add pc=t, but still out.fa contains all 4 sequences, including the two 100% substrings (seq2,seq3) of seq1.

Would be great if you found time to look into this!

[Brian Bushnell]

Hi Shini,

am=t (100% identical sequences SHOULD be absorbed)
ac=t (100% contained sequences SHOULD be absorbed)
fo=f (I do not want to cluster/assemble sequences)

Even though you want to absorb them, things that are absorbed will not show up in the dot output. Also, I think you have to set "fo=t" because the dot file prints out all the overlaps that are found (which include containments). Dedupe runs in 3 phases:

1) Load reads and find exact matches.
2) Find containments.
3) Find overlaps.

Phase 1 and 2 do not actually generate the information needed for the dot file, only phase 3. So, you need "am=f ac=f fo=t". Sorry!

As for the speed, yep, I put a lot of effort into making it very fast. I wrote it mainly because we used Minimus2 for deduplication our merged assemblies, which for large metagenomes would take days and then often crash or run out of time.

Anyway, I'll take a look at it this weekend and see if that's something that's really easy to add.

[Shini Sunagawa]

Hi Brian,

Thanks for you quick response, really appreciate you took some of your limited time!

I hoped the dot file might be just what I was looking for, so I ran

dedupe.sh in=in.fa out=out.fa outd=outd.fa am=t ac=t fo pc dot=dot.txt

For clarification: I set
am=t (100% identical sequences SHOULD be absorbed)
ac=t (100% contained sequences SHOULD be absorbed)
fo=f (I do not want to cluster/assemble sequences)
dot=dot.txt

However, the dot file appears to contain the indices of the unique input sequences only, but what is missing are the indices of all contained sequences, that is, you get the same result irrespective of the number of different (or same) substring sequences of a given sequence (while I would be interested in ALL the indices/ids of the different sequences).

To give you an idea on the speedup, I tested two programs that would output the information of interest and they both took about 1 hour to dereplicate 1 M sequences, while dedupe does it in 25 seconds (!). To be fair, I am not sure how much overhead it is to keep track of the ids/indices and to print them.

Anyway, the results seem equivalent, so it would be fantastic if this became a feature one day!

Best,
Shini

[Brian Bushnell]

Hi Shini,

Dedupe technically has the ability to do that, although currently it doesn't actually do it. Hmmm...

It does have a "uniqueonly" flag; if you use that, it will only send sequences to "out" if they are unique, and anything nonunique gets sent to "outd" (so if sequences A and B are duplicates, both get sent to outd). At a minimum, that will make the problem much smaller so that maybe it can be dealt with by a slower program that has the desired output format, if there is one. For mostly non-redundant data, that should accelerate things by a huge factor.

I don't have very much free time these days, so while I agree that this would be useful, I'm not sure when I'd have an opportunity to implement it. Therefore - I've added it to my TODO list, but that list is pretty long. Do you happen to have a programmer at your disposal? If so, Dedupe's dot output (for example, "dot=overlaps.dot") does contain all of the information of which sequences overlap which other sequences (if you run Dedupe with the correct flags - "am=f ac=f fo dot=dot.txt"), and can be transformed into the text output you're looking for. It would be very difficult for a non-programmer to do the transformation. Though, it is readable by non-programmers, so you might take a look.

Sorry I can't help further, though I will try to add that functionality when possible.

[Shini Sunagawa]

dedupe.sh

Dear Brian,

I have been looking for a tool that would quickly dereplicate (100% containments) nucleotide sequences and track for each unique sequence the identifiers of the removed duplicates.

Something like:

dedupe.sh in=in.fa out=out.fa outd=outd.fa mid=100 mop=100

where:

in.fa:
seq1
seq2 (contained in seq1)
seq3 (contained in seq1)
seq4

out.fa:
seq1
seq4

outd.fa:
seq2
seq3

I am interested in:
seq1<tab>seq2,seq3
seq4

dedupe.sh does a fantastic job in returning out and outd, but I cannot find any option that would return the information I am interested in. Is this something that I am missing? Otherwise, I believe this could be a great feature, since compared to other tools that return this information, dedupe is so much faster.

Best,
Shini

[Brian Bushnell]

For now I only want an estimation of duplicates for some HS4000 paired end data so I will just omit "out=". I assume a single "out=" will result in an interleaved file (for PE input) that would have to be resolved afterwards?

"out" is synonymous with "out1". So if the input is paired, "out=x.fq" or "out1=x.fq" will produce interleaved output; "out=x1.fq out2=x2.fq" or "out1=x1.fq out2=x2.fq" would both produce dual-file output (once Dedupe parses "out1" correctly).

[GenoMax]

Thanks for finding that! Looks like a typo crept in; I'll fix that in the next release. In the mean time, you can use "out=" instead of "out1=".

For now I only want an estimation of duplicates for some HS4000 paired end data so I will just omit "out=". I assume a single "out=" will result in an interleaved file (for PE input) that would have to be resolved afterwards?

[Brian Bushnell]

Hello, I am working with genomic data belonging to mammals. I have been aligning raw reads to all the genes within an organism's chromosomes. While working with this data set, I have noticed some odd values when calculating coverage with BBTool's built-in pileup feature. I would really appreciate it if anyone knew how to interpret this.

I don't see anything odd in general - high variances are not really unusual. But the -1 values are definitely problematic and should not happen. I'll investigate and get back to you.

[Brian Bushnell]

Thanks for finding that! Looks like a typo crept in; I'll fix that in the next release. In the mean time, you can use "out=" instead of "out1=".

[GenoMax]

@Brian: I am having trouble using dedupe.sh with paired end reads (using BBMap v.35.59).

Unknown parameter out1=sampleID_L001_R1_001.fastq.gz

     at jgi.Dedupe.<init>(Dedupe.java:383)

     at jgi.Dedupe.main(Dedupe.java:80)

command I am using is in the format you had posted

$ dedupe.sh in1=read1.fq in2=read2.fq out1=x1.fq out2=x2.fq ac=f

Single-end files work fine.

[JulesWinchester]

Here are the files I promised. I also attached the whole statistical output, in case anyone wants to see.

Thanks in advance!

[JulesWinchester]

Hello, I am working with genomic data belonging to mammals. I have been aligning raw reads to all the genes within an organism's chromosomes. While working with this data set, I have noticed some odd values when calculating coverage with BBTool's built-in pileup feature. I would really appreciate it if anyone knew how to interpret this.

Here is the issue:

Data being used:
-CDS (W/ Introns, not spliced) of genes extracted from a Chromosome Genbank file (found @ NCBI ftp://ftp.ncbi.nlm.nih.gov/genomes/P...odytes/CHR_01/).
-Raw reads collected from ENA (Adapter free WGS sequences)

BBMap Commandline used to calculate coverage:
bbmap.sh in1=Chimp_1.fastq in2=Chimp_2.fastq ref=ChimpChrom1.fa local=t nodisk covstats=Chrom1Stats.txt covhist=Chrom1Hist.txt

Results for the chimp. The results for all other mammals I've worked with are very similar.
Coverage - 43.97
Standard Deviation - 507.47

I proceeded to look at the statistical output produced by BBMap to see if the SD values were caused by specific gene sequences; there I saw that a lot of genes had median folds equal to 0, others that had surpassed 1,000, and even some with negative values (-1). Why does this happen, is this normal? If not, is there a way to fix this?


I have attached the stats file produced by BBMap on this post if anyone would like to see how the output looks like. I took the liberty of parsing the file based on median fold values (equal to or more than 100), only keeping the columns belonging to gene ID and Median Fold. On a follow-up post, I will add attach a text file containing the low median fold values (anything below 100,).
Note: The coverage and SD tendency remains, even when using a reference fasta containing only the longest isoform per gene (1 isoform per gene).

[Brian Bushnell]

Hi everyone,

I wanted to mention that there is now some additional documentation for BBTools in the form of guides/tutorials, in the /docs/ folder. Currently there are guides for BBDuk, BBMerge, Seal, Tadpole, Reformat, Dedupe, BBNorm, and Taxonomy, and I plan to add more in the near future. There's also an overview of the general usage of all BBTools (UsageGuide.txt) and a list of all the commonly-used tools with a brief description (ToolDescriptions.txt).

I hope these will be useful, and please let me know if anything is unclear or needs to be expanded, or there is a common use that the guides don't address.

-Brian

[GenoMax]

@babine: Are you only "borrowing" the reference i.e. you do not have raw/original data for the "not-so-accurate" reference?