ReadLen=36

[PanosFirmpas]

Hi,

You have set your biasFree script analysis to "ReadLen=36".
Does that not mean that if someone uses this with a bam file of different read length, the analysis will be wrong ?

[HouyuZhang]

Hi Panos,

Thank you for this question.

Actually, the ReadLen parameter is used to calculate genome mappability. This parameter can be trace-backed to the seqOutBias tallymer --read-size parameter, and further to the gt tallymer mkindex -mersize parameter. The reasons why we use k=36 as the default value are as follows:

• A benchmark study tested different read lengths and reported that when using k=36 compared to k=75, the difference of uniquely mapped positions is ~2% for smaller genome size like C.elegans and D.melanogaster, ~7% for larger genome size like human and mouse. But the time cost is about 6-fold.
• Another study made a more detailed benchmark and found the ratio of uniquely mapped positions reaches the plateau after k>20.
• So, this is a trade-off. Considering our pipeline is intended to generalize for many species, we choose the least acceptable length k=36 to pre-calculate the mappability files. That's also the default values seqOutBias chose in their analysis.

Yes, you are right, using the specific read length from your experiment (usually 150 bp nowadays) will lead to a more accurate result, but may need more computational cost.
Maybe it's a good practice to modify the ReadLen parameter as a user-defined parameter in a new version.

Glad to discuss more.

Houyu Zhang

[PanosFirmpas]

Thank you for the swift answer !

So, while computing the genome mappability for the exact readLength would be more accurate, the 36bp calculation is still a "good-enough" approximation ?

Thank you also for the offer to discuss more, I will probably come back with further questions :)

Best Regards,
Panos

[HouyuZhang]

Hi Panos,
Thank you for the further question.
I have checked our previous analytic script. We used k=51 for the mESC ATAC-seq data in our paper since the exact read length is 51bp. We have not tested other k values.

Tn5InsertPrefer/5_BiasFreeATAC_mESC.sh

Line 20 in 518572f

seqOutBias ${FASTA_FILE} ${name}_plus.bam --kmer-mask ${plus_mask} --bw=${name}_plus_${plus_mask}-mer.bigWig --shift-counts --read-size=51

As for your question about the "good-enough" approximation, I think maybe a titration using varied k values will tell us whether this parameter affects the final results of our pipeline. However, based on the theories from the papers I listed, I personally think using k=36 or a larger value will give pretty similar results. I will update the testing result here when possible.

Anyway, following your suggestions, I will open the ReadLen parameter to users in the next version.

Best wishes,
Houyu Zhang