FASTQ and BAM give discordant results

[vincentwalter]

Expected Result

Starting an analysis from FASTQ files and a BAM file generated from these FASTQ files should produce the same result.

Actual Result

A small fraction of reads is displayed to be unique to either the FASTQ or BAM input.

Exact MiXCR commands

vwalter@work:~/projects/playground$ R1=path/to/R1.fastq.gz
vwalter@work:~/projects/playground$ R2=path/to/R2.fastq.gz
vwalter@work:~/projects/playground$ OUT=path/to/output
vwalter@work:~/projects/playground$ samtools import -i -1 $R1 -2 $R2 --order ro -O bam,level=0 | \
    samtools sort -@10 -M -o "${OUT}.bam" -
vwalter@work:~/projects/playground$ mixcr analyze cellecta-human-rna-xcr-umi-drivermap-air -Xmx50g     $R1 $R2 "${OUT}_fastq"
vwalter@work:~/projects/playground$ mixcr analyze cellecta-human-rna-xcr-umi-drivermap-air -Xmx50g     "${OUT}.bam" "${OUT}_bam"
vwalter@work:~/projects/playground$ mixcr clonesDiff -c "${OUT}_fastq.clns" "${OUT}_bam.clns"
Unique clones in cloneset 1: 439 (0.19%)
Reads in unique clones in cloneset 1: 2584 (0.05%)
Unique clones in cloneset 2: 131 (0.06%)
Reads in unique clones in cloneset 2: 1025 (0.02%)

Software Versions

vwalter@work:~/projects/playground$ mixcr --version
MiXCR v4.6.0 (built Sat Dec 09 20:48:42 CET 2023; rev=c9fafa41fe; branch=no_branch; host=fv-az520-378)
RepSeq.IO v2.4.0 (rev=5b5b4475b8)
MiLib v3.4.0 (rev=d42f072b23)
Built-in V/D/J/C library: repseqio.v4.0

Library search path:
- built-in libraries
- /home/vwalter/projects/playground/.

vwalter@work:~/projects/playground$ samtools --version
samtools 1.20
Using htslib 1.20
Copyright (C) 2024 Genome Research Ltd.

Samtools compilation details:
    Features:       build=configure curses=yes 
    CC:             gcc-11
    CPPFLAGS:       
    CFLAGS:         -Wall -g -O2
    LDFLAGS:        
    HTSDIR:         
    LIBS:           
    CURSES_LIB:     -lncursesw

HTSlib compilation details:
    Features:       build=configure libcurl=yes S3=yes GCS=yes libdeflate=no lzma=yes bzip2=yes plugins=no htscodecs=1.6.0
    CC:             gcc-11
    CPPFLAGS:       
    CFLAGS:         -Wall -g -O2 -fvisibility=hidden
    LDFLAGS:        -fvisibility=hidden 

HTSlib URL scheme handlers present:
    built-in:    preload, data, file
    S3 Multipart Upload:         s3w, s3w+https, s3w+http
    Amazon S3:   s3+https, s3+http, s3
    Google Cloud Storage:        gs+http, gs+https, gs
    libcurl:     imaps, pop3, gophers, http, smb, gopher, sftp, ftps, imap, rtmpte, smtp, smtps, rtsp, rtmpe, scp, ftp, telnet, mqtt, rtmp, ldap, https, ldaps, rtmps, rtmpt, pop3s, rtmpts, tftp, smbs, dict
    crypt4gh-needed:     crypt4gh
    mem:         mem

[mizraelson]

Hi, thank you for reporting this! Can you please try the latest develop version and compare the outputs. We have identified and fixed one bug that have potentially caused the issue. Let me know if it helped!

[vincentwalter]

It's better but not completely solved. I reran with the release (c9fafa41fe) and the develop version (b0900b576f):

vwalter@work:~/projects/playground$ R1=path/to/R1.fastq.gz
vwalter@work:~/projects/playground$ R2=path/to/R2.fastq.gz
vwalter@work:~/projects/playground$ OUT=path/to/output
vwalter@work:~/projects/playground$ samtools import -i -1 $R1 -2 $R2 --order ro -O bam,level=0 | \
    samtools sort -@10 -M -o "${OUT}.bam" -
vwalter@work:~/projects/playground$ mixcr analyze cellecta-human-rna-xcr-umi-drivermap-air -Xmx100g $R1 $R2 "${OUT}_fastq_c9fafa41fe"
vwalter@work:~/projects/playground$ ./mixcr analyze cellecta-human-rna-xcr-umi-drivermap-air -Xmx100g $R1 $R2 "${OUT}_fastq_b0900b576f"
vwalter@work:~/projects/playground$ mixcr analyze cellecta-human-rna-xcr-umi-drivermap-air -Xmx100g "${OUT}.bam" "${OUT}_bam_c9fafa41fe"
vwalter@work:~/projects/playground$ ./mixcr analyze cellecta-human-rna-xcr-umi-drivermap-air -Xmx100g "${OUT}.bam" "${OUT}_bam_b0900b576f"

The differences get smaller:

vwalter@work:~/projects/playground$ ./mixcr clonesDiff -c "${OUT}_fastq_c9fafa41fe.clns" "${OUT}_bam_c9fafa41fe.clns"
Unique clones in cloneset 1: 439 (0.19%)
Reads in unique clones in cloneset 1: 2584 (0.05%)
Unique clones in cloneset 2: 131 (0.06%)
Reads in unique clones in cloneset 2: 1025 (0.02%)
vwalter@work:~/projects/playground$ ./mixcr clonesDiff -c "${OUT}_fastq_b0900b576f.clns" "${OUT}_bam_b0900b576f.clns"
Unique clones in cloneset 1: 155 (0.07%)
Reads in unique clones in cloneset 1: 962 (0.02%)
Unique clones in cloneset 2: 111 (0.05%)
Reads in unique clones in cloneset 2: 721 (0.01%)

However, is it expected, that the clns files differ this drastically between versions ?

vwalter@work:~/projects/playground$ ./mixcr clonesDiff -c "${OUT}_fastq_c9fafa41fe.clns" "${OUT}_fastq_b0900b576f.clns"
Unique clones in cloneset 1: 226519 (99.98%)
Reads in unique clones in cloneset 1: 4817261 (99.99%)
Unique clones in cloneset 2: 226488 (99.98%)
Reads in unique clones in cloneset 2: 4817734 (99.99%)
vwalter@work:~/projects/playground$ ./mixcr clonesDiff -c "${OUT}_bam_c9fafa41fe.clns" "${OUT}_bam_b0900b576f.clns"
Unique clones in cloneset 1: 226222 (99.99%)
Reads in unique clones in cloneset 1: 4798599 (99.99%)
Unique clones in cloneset 2: 226455 (99.99%)
Reads in unique clones in cloneset 2: 4817078 (99.99%)

[mizraelson]

Hi,

If you could share the input files, I can investigate further. Generally, for reads that have low similarity to the reference, the randomized part of the alignment algorithm can play a bigger role. There is a seed that is calculated from the ID and the header of every read; this seed will differ in BAM and FASTQ formats. For reads that align well with the reference, this does not affect the outcome.

Regarding the difference between the two versions, some algorithms have improved, and the output results can slightly differ and be more reliable. From the numbers you have shared, the difference is less than 0.5%.

[vincentwalter]

I sent a link to the files to support@milaboratories.com.

About the difference in clonesets between the two versions. Where do you get the less than 0.5%? The output says Unique clones in cloneset 1: 226222 (99.99%) so I was under the impression that almost no overlap exists between the clonesets generated with the different versions.

[mizraelson]

Hi, thank you for sharing the data. We have identified the reason for the remaining difference between BAM and FASTQ. It is due to the way nucleotide wildcards are resolved. This algorithm requires randomization which comes from different seeds in bam and fastq files. We will fix this issue in the upcoming release but it only affects a very small fraction of reads.

As for the difference between the versions, the clonesDiff function relies on full sets of gene hits. For example: 1) IGHD2-2*00(40),IGHD3-22*00(36),IGHD4-11*00(35) 2) IGHD2-2*00(40),IGHD3-22*00(36),IGHD1-26*00(35)

The clones above will be considered different because of the last D gene hit. But if you compare the samples using the best hit for every gene (IGHD2-2*00 in the example above), all clones will be identical.

[vincentwalter]

Thank you for looking into it!