Discrepancy between the number of taggable reads and the number of records in the bam file

[lebelm]

Hello,

I have been trying barQC to recover additional reads with out of phase barcodes and as far as I can tell it has been running well (I processed all the reads with in phase barcodes separately to reduce runtime).
However, I have been surprised to find that the number of taggable reads written at the end of the barQC.log and tha actual number of records in the output BAM files do not match. This is true for all the libraries I have tested

For instance I get in the barQC.log:

...
INFO - Final info:
                        Number of original reads: 20787064
                        Number of taggable reads: 5540482
                        Efficiency: 26.65%
                        Number of unique barcode combinations: 396884
...

But if I run samtools view -c MYRUN_barQC.bam I get only 612,359 records, which is almost 10 times less than I would expect.
There is no error written in the barQC.log or barQC_debug.log or during the run.

The environment has been installed with conda, here are the packages:

Packages in the environment

Name Version Build Channel

libffi 3.4.4 h6a678d5_1
libgcc-ng 11.2.0 h1234567_1
libgomp 11.2.0 h1234567_1
libstdcxx-ng 11.2.0 h1234567_1
libuuid 1.41.5 h5eee18b_0
matplotlib 3.9.3 pypi_0 pypi
ncurses 6.4 h6a678d5_0
numpy 2.2.0 pypi_0 pypi
openssl 3.0.15 h5eee18b_0
packaging 24.2 pypi_0 pypi
pandas 2.2.3 pypi_0 pypi
pillow 11.0.0 pypi_0 pypi
pip 24.2 py312h06a4308_0
psutil 6.1.0 pypi_0 pypi
pyparsing 3.2.0 pypi_0 pypi
pysam 0.22.1 pypi_0 pypi
python 3.12.3 h996f2a0_1
python-dateutil 2.9.0.post0 pypi_0 pypi
pytz 2024.2 pypi_0 pypi
readline 8.2 h5eee18b_0
seaborn 0.13.2 pypi_0 pypi
setuptools 75.1.0 py312h06a4308_0
six 1.17.0 pypi_0 pypi
sqlite 3.45.3 h5eee18b_0
tk 8.6.14 h39e8969_0
tzdata 2024.2 pypi_0 pypi
wheel 0.44.0 py312h06a4308_0
xz 5.4.6 h5eee18b_1
zlib 1.2.13 h5eee18b_1

Do you have any suggestions as to what is happening? Is this an expected behaviour?

Best,

Marie

[m-rossello]

Hi Marie,

First of all, I'm sorry for the delayed response — thank you for your patience, and for taking the time to put together such a clear and detailed report. It really made it much easier to track down the issue.

We identified the likely cause of the discrepancy. The number of taggable reads reported in the log is counted from read2 (where barcodes are extracted and corrected), but the reads actually written to the BAM file are sourced from read1. The matching between the two relies on read names being identical. If your FASTQ files use paired-end suffixes in the read headers — /1 for read1 and /2 for read2 — the names won't match and reads will be silently skipped during the tagging step, producing a much smaller BAM than expected.

You can quickly check whether this is the case by running:

head -1 your_read1.fastq
head -1 your_read2.fastq

If you see something like @READ_NAME/1 and @READ_NAME/2, that confirms the issue.

This is a bug on our end — the read name comparison should be suffix-aware, and we will push a fix. In the meantime, you can work around it by stripping the suffixes from your FASTQ headers before running barQC:

sed 's/\/1$//;s/\/2$//' read1.fastq > read1_fixed.fastq
sed 's/\/1$//;s/\/2$//' read2.fastq > read2_fixed.fastq

Please don't hesitate to follow up here if the workaround doesn't resolve it, or if you run into anything else — I'm happy to keep looking into it with you!

Best,
Maria

[lebelm]

Hi Maria,

Thank you for taking the time to look into it!

I ran pairfq makepairs (https://github.com/sestaton/Pairfq) before barQC, so the read1 and read2 files should contain paired reads only.
Given that, here is what my headers look like:

head -1 R1_notinphase.fq
@A00709:546:HMLY5DSX5:4:1101:30590:1016 1:N:0:CAGATC+GCGATCTA

head -1 R2_notinphase.fq
@A00709:546:HMLY5DSX5:4:1101:30590:1016 2:N:0:CAGATC+GCGATCTA

So the first part of the name is the same, but the second part does contain the read information.
I'll try stripping the second part of the read header and see if it fixes the issue on my end.

What puzzles me is that if the issue came from the read information part, then I would expect an empty output, as no read name in R1 would match with R2 (if they all contain the read information, which they should)...

I'll keep in touch once this has run!

Best,

Marie

[lebelm]

Hi,

I ran barQC again after removing the second part of the names with awk

awk '{FS=" "; if ($0 ~ /^@/) print $1; else print $0}' R1.fq > R1_awk.fq
awk '{FS=" "; if ($0 ~ /^@/) print $1; else print $0}' R2.fq > R2_awk.fq
python /omics-storage/lab-solana/softw/barQC/barQC.py -f1  R1_fixed.fq -f2  R2_fixed.fq -o barQC -b /PATH/TO/barQC/barcode_list -t 32 --stats -v

I still see the discrepancy in the number of "taggable" reads and the total number of reads in the output.

I have checked using

grep "^@" R1_awk.fq > R1_name_awk.lst
grep "^@" R2_awk.fq > R2_name_awk.lst

#count numbers of lines common in file 1 and file 2
comm -12 R1_name_awk.lst R2_name_awk.lst | wc -l

That indeed all corrected files contain all reads with the same name.

So this does not seem to be the source of my issue...

[m-rossello]

Hi Marie,

Thank you for testing this so thoroughly and for the careful follow-up. Your debugging has been really helpful in narrowing things down.

You're right that the read name mismatch is not the source of the issue. We're still investigating the root cause, and we want to be transparent that we don't have a definitive answer yet.

In the meantime, we'd suggest not relying on the "taggable reads" count in the barQC log for your final statistics. Instead, you can get accurate numbers directly from the tagged BAM using barcode_QC_UMI_seq.py, which will give you a more reliable picture of your data.

We haven't forgotten about this — we'll update this issue as soon as we've identified the bug. And please don't hesitate to share any additional observations in the meantime, your reports have been extremely useful!

Sorry again for the slow turnaround, and thanks for your patience.

Best,
Maria