Raw mean coverage difference with salting

[abolia]

Hi Daniel,

Can you please explain me how is the raw mean coverage is calculated with BMFtools depth. When I change the salting parameter, the raw mean coverage also changes. If its just the mean coverage of total number of raw reads (from fastq) it should be same irrespective of salting used. Therefore, I wanted to check if there is a different logic behind calculating the raw mean coverage.

Thanks,
Ashini

[dnbaker]

That could depend on filter settings. How do those compare to the actual numbers of raw reads?

[abolia]

Sorry, I asked the wrong question. I meant to ask why is raw mean coverage changing based on different salting parameter used. Shouldn't that be same as its raw reads mean coverage.

The raw reads are consistent. Sorry about the confusion.

Thanks,
Ashini

[dnbaker]

As salting increases, so spurious unique observations increase due to errors in barcode reading, which should have been collapsed into another family.

Also, depending on how salting is performed affects how collapsed reads will align, so that could also be affecting it.

[abolia]

I understand that this could affect collapsed mean coverage, but I am still not getting how raw mean coverage gets affected by this. Isn't the raw coverage before collapsing so shouldn't it be consistent, irrespective of salting used during collapsing.

Thanks,
Ashini

[dnbaker]

Raw mean coverage is repeating the collapsed coverage calculation with weighting by family size. We haven't actually aligned the raw dataset.

[abolia]

Hey Daniel,

Thanks for the reply. That makes sense now. I have another question on the total founding reads in BMFtools famstats. Does these represent raw reads ? When I compare them with a simple bwa aligned bam file (without any collapsing) I get 60,488,497 reads. However, the total founding reads in BMFTools famstats are only 25,396,920. So, I am wondering why are these so different?

Thanks again,
Ashini

[dnbaker]

Great question. What other processing have those files had? Checking the large validation dataset we used for 1.0/1.1 ensured that reads were all accounted for, so at least at that point, we weren't losing any records in the process.

In particular, was there any filtering done on the dataset (minimum family size or mapping quality)?

[abolia]

No, I did not filter any reads before calculating the coverage stats. The steps follow BMFtools, then skewer to mask adapters, bwa alignment and BMFtools coverage calculations. No filtering at any step.

Ashini

[dnbaker]

Thank you for reporting the issue, and I'm happy to help.

I'm a little unsure about how you got an uneven number of raw reads for paired-end sequencing, but it seems like you're definitely losing reads along the way.

I have a couple of questions to ask so I can help better.

First, can you reproduce the issue on a small (or smaller?) subset?

Second, would you be willing/able to provide a script listing commands performed and a dataset which reproduces the issue? (If you'd prefer it not be publicly available, feel free to email me it.)

[abolia]

Can you provide me your email address?

Thanks,
Ashini