Read statistics

[johnne]

Hi!
I've run QC, assembly and annotation on my samples and started looking into the assembly stats and annotation coverage files. So far running atlas has been smooth but I find it a bit difficult to follow read numbers through the pipeline, e.g. to make number of reads output from QC add up to fraction assembled/aligned reads at the various steps.

For example, in the stats/combined_contig_stats.tsv file all my samples have 100% assembled reads which surprises me. From what I can tell the stats.sh script from BBTools is used to generate the stats but I'm struggling to understand how this is done. Any help or details on this step would be very much appreciated.

Also, in the minimap2 alignment step, or in some of the postprocessing of the bam files, how are multimapped reads handled? Is there some way to affect this behaviour with config parameters, or to set mapping quality filtering somehow?

[johnne]

Digging a bit more I can see that the combine_contig_stats rule takes output from stats.sh, prodigal gene predictions and pileup.sh logs then uses the Mapped reads and Percent mapped statistics as values for Assembled_Reads and Percent_Assembled_Reads. With no filtering done on map quality (neither in the alignment with minimap2 or pileup rule) those numbers could maybe be a bit misleading, or what are your thoughts on this? One way to introduce quality filtering for alignments could be in the pileup.sh stage using the minmapq parameter.

[johnne]

Testing to run the pileup.sh command as in the pileup_contigs_sample rule but with the minmapq flag set to 10, 20 or 30 shows a drop from 100% to 98.897, 98.357 and 97.700% mapped reads for one sample tested. So not a huge difference though.

[SilasK]

I think I should have a look at this. Yes, I didn't add a filtering step when I replaced bbmap with minimap.
Something else I'm not sure is that the unmapped reads are correctly considered for the percent mapped reads.

[SilasK]

PRs are always welcome.
Do you think filtering based on minimum identity should also be necessary?

[johnne]

Feels like maybe map quality should be enough for this purpose. And also I'm not sure what a resonable identity threshold would be, but maybe you have a better feeling for this? Mostly the map quality will help with removing the very low quality/ambiguous mappings.

What I was also thinking about was whether to allow only using primary alignments in atlas. I see now you pass secondary=t to pileup.sh which could be made user-configurable.

[johnne]

Actually now that I think of it, using secondary alignments might be preferrable in most cases here. If there are similar organisms in multiple samples and samples are cross-mapped e.g. in binning and when quantifying the Genecatalog?

[johnne]

Implemented minmapq and made a PR #647. I've tested using your example data and it works. Setting minimum_map_quality: 10 in a config file resulted in 98.818% and 99.628% Percent_Assembled_Reads reported for sample1 and sample2 megahit assemblies, respectively.

[johnne]

@SilasK , please note that the vamb developers stress that alignments should not be filtered for low-quality mappings prior to binning with vamb. I know that atlas doesn't do that for vamb and not for MetaBat2 either, but as far as I can tell the coverage from pileup.sh is passed to CONCOCT and MaxBin2. It may be best to make the minimum_map_quality parameter specific for the assembly filtering/report, Genecatalog and genome coverage rules. Either by making a separate parameter such as minimum_map_quality_binning or excluding it entirely from the binning rules.

[johnne]

See #650