Multiqc: BBmap BBsplit output not found

[brucemoran]

Hi Phil,

MultiQC v1.5 is not finding BBsplit output, as defined here. No errors thrown etc, file format is same, see attached. Not sure what else I can report, apologies.

Bruce.

test.statsfile.txt

[ewels]

Thanks @brucemoran - I'll take a look when I get a chance.

[brucemoran]

Hi Phil,

had an hour to dig around in this today. Found a workaround by reformatting BBsplit refstats output into 'disambiguate'-style *_summary.txt. My use case is to QC number of reads per each genome (reasonable assumption that this will be desirable for most using BBsplit): simply specifying total reads to genomeX, genomeY, and those that are ambiguous. Unsure what is involved (I have no Python) but may be simpler to parse refstats; statsfile is 52 lines, 5 cols; refstats n lines for n genomes (+ header), 7 cols.

Really enjoying MultiQC btw, this is the most tempted I've been to learn Python yet! Also kudos on your Nextflow work.

Bruce.

bbsplit.refstats.txt

[ewels]

Thanks Bruce! I don't know the BBsplit code myself very well, I think it was written by @boulund, but I'll take a look when I get a chance. Not much time these days and a fairly big MultiQC backlog currently, but I'll do my best.

Really enjoying MultiQC btw, this is the most tempted I've been to learn Python yet!

Brilliant! I've found that MultiQC is a really nice project to get people started on Python actually (hackathons, new people in the office). Plenty of existing code to work from and nice small modular tasks to do.

Phil

[ewels]

Aha! Wait, I've found the error. The code for BBsplit output is not yet written 😆

The BBmap module doesn't have complete support for all tools yet, see http://multiqc.info/docs/#bbmap

I think this was intended but not yet done.

Phil

[brucemoran]

Hah, that's a problem alright. That doc specifies refstats is available BTW, but in bbmap_filetypes.py there is a gap before rpkm where refstats should probably go, and no entry in bbmap.py.

Bruce

[boulund]

Hi!

How fun that someone is actually using the BBMap MultiQC module! :D
The BBMap suite has so many different types of output that we just didn't have the time to implement all of them when we wrote the module in the first place about a year ago. I wanted to something to give me basic QC summaries (something similar to FastQC output) and we just never got around to implementing support for everything, hoping that users eventually would complain about what was missing, and here we are 👍.

I haven't used BBSplit for any real work yet, so I'm not entirely sure what information users of that typically would want to summarize in a MultiQC report. You mentioned just summarizing the number of reads per reference as one thing. Do you have any further suggestions for how to visualize the output in the report? I think it'll be fairly straightforward to throw something together to do what you want. We never settled on how to best visualize e.g. statsfile output in the MultiQC report.

It looks like we just completely missed refstats. It should be fairly straightforward to implement support for the refstats 7-column format anyway. Unfortunately I think the BBMap module is the least beginner-friendly module to start with for a budding pythonista. The parser is quite esoteric...

[brucemoran]

Yes, always seems like nobody will see your work, then they do!

I use BBsplit rarely (xenograft mouse models), but I use BBDuk for adapter and quality trimming, it's really fast, and flexible with inputs/outputs formats etc. Couldn't really get my head around the BBMap module setup. A shame, as the individual module for disambiguate which I hacked could quite easily be adapted, I just have no idea where to import and how to define what is required unless at a very basic level. Also limited time for this!

If you can get refstats working that would be great, really want to be able to see how much of each reference is mapped to, and then also how much is ambiguous, e.g. unambiguousReads per reference (each on its own line), and then ambiguousReads (same for each reference).

Thanks,

Bruce.

[boulund]

I can have a quick look to get an estimate for how much work it would entail.

[boulund]

Ok, now I've had a look. There are a few things required in order to implement support for bbsplit's refstats output:

• A custom parser for the (by BBMap standards) strange output format. BBMap outputs typically contain histogram data, and the non-key/value parser used is not up to the task at the moment. It would need some extensions to handle the more varied content of the refstats table.
• A custom function to plot or produce a table for the data you want to include in the report. The primary thing here would be to decide what to include in the report and in what format (table, plot, plots in separate tabs?), as that will dictate our data structure a bit.

There is one strange thing with the example refstats file you have supplied. It looks like it uses spaces to separate columns. Do you know if this is always so for the refstats output from bbsplit.sh? And are the number of spaces in the header constant, or do they vary with size of the values printed?

[brucemoran]

Apologies, the refstats file I uploaded was required to be .txt by Github which I think must have converted whitespace to spaces. The original file is definitely tab-delimited.

For my use case, output to general stats and bargraph.plot would be best. Personally I like to line up all the stats per sample, so we can see M Seqs, M Aligned, % Aligned, % Dups etc., and then include M genomeX, M genomeY, M ambiguous, % genomeX, % genomeY, % ambiguous columns from BBsplit refstats. For plotting, Picard MarkDuplicates or Star genecount_chart are good examples, so something like:

keys = OrderedDict()
     keys['M_genomeX'] = { 'color': '#2f7ed8', 'name': 'M Seqs Mapping to GenomeX' }
     keys['M_genomeY'] = { 'color': '#0d233a', 'name': 'M Seqs Mapping to GenomeY' }
     keys['M_ambiguous'] = { 'color': '#492970', 'name': 'M Seqs Mapping Ambiguously' }

Idea is to show how much of sequence went to each genome, and what was ambiguous. Does that make sense?

Really appreciate your help on this.

Bruce.

[boulund]

I think that sounds reasonable. I can probably have a look at this sometime during the weekend. I'll report back as soon as I've had time to have a proper sit down with this.

[boulund]

I'm a bit concerned about how this will scale if BBSplit is used for something like binning metagenomic reads to a very large set of reference genomes (e.g. if someone wants to bin to NCBI RefSeq, which contains the genomes of some 77k organisms, you could start to imagine how badly this could go for a MultiQC report). Have any suggestions for how to handle this extreme variability in the number of references used?

[ewels]

+1 for considering edge cases! 🎉

Could just add a conditional that if there are more than x genomes (20? depends on the plot I guess), that we take the top x represented ones, as determined by calculating a total count per reference across all samples... And add a warning to the description. eg. a bootstrap alert box:

<div class="alert alert-info">
    <strong>Note:</strong> Only the top 20 reference genomes are shown, sorted by total total count across all samples.
</div>

[boulund]

We'll probably have to use something much lower than 20. If we're potentially putting three columns of data for every reference sequence into the main stats table, then a limit of 3 or something might be more reasonable. It might make more sense to put a set of histograms in the BBSplit section of the report. For low numbers of reference sequences you'll probably instantly see which one is which anyway, and for large numbers you'll at least get a sense of the distribution of binned reads. I guess most people will want to dig into the raw bins anyway, so just presenting some summary data should be enough. Giving min, mean, median, and max values could be nice as well.

[ewels]

If we're potentially putting three columns of data for every reference sequence into the main stats table

This is way too many I think. If they need to be in a table then that's fine, but it should probably be a separate table in its own section.

We could also have a stacked barplot with different tabs for read counts, % aligned, dups etc. Each bar stacked per genome. I think this is pretty much what you're suggesting, right..?

[boulund]

I like the idea of the stacked barplots. I'll probably try that first! I never had the time to look at this like I intended yesterday, hopefully maybe later next week...