CallDuplexConsensusReads Raw Reads Filtered Due to Not Enough Reads

[ARE4S]

Hello everyone,

with this I'm referring to issue #556. I'd like to call Duplex Consensus reads using fgbio as well. Unfortunately, I'm facing the exact same issue as mentioned above.

The samples I'm using are paired end sequenced and tagged with 6bp UMI and 11bp stem resulting in a cigar like 6M11S+T 6M11S+T.
Afterwards, I'm extracting UMIs with fgbio-2.0.2 ExtractUmisFromBam, followed by GroupReadsByUmi and CallDuplexConsensusReads
java -jar fgbio-2.0.2.jar ExtractUmisFromBam \ -i /path_unmapped.bam \ -o /path_unmapped_umi_extracted.bam \ -r 6M11S+T 6M11S+T \ --molecular-index-tags=ZA ZB \ --single-tag=RX

java -jar fgbio-2.0.2.jar GroupReadsByUmi \ --input=/path_umi_extracted_aligned_merged_filtered.bam --output=/path_umi_grouped.bam \ --strategy=paired --edits=1

java -jar fgbio-2.0.2.jar CallDuplexConsensusReads \ --input=/path_umi_grouped.bam --output=/path_consensus_unmapped.bam \ --error-rate-post-umi 30 --error-rate-pre-umi 45 \ --min-input-base-quality=20

This leads to the following outputs (examplary Read-ID extracted)

original fastq read 1:

@M01950:199:000000000-CN4MR:1:1106:27160:16985 1:N:0:2 CACAGGAGTAGCTCAAGAGCTCGTCTGGTAACTATCCATATATGAATAGCCCATTATGGACTCGGGTGAGCTGGTATAGGGGTCTGGGTACTCAGACTTGATGGCCCGGCTAGGAAAGTGGCCATATGTTTGGTAACCTTGCAGGCTGCCG + CCCCCCFCFFFFGGGGGGGGGGGGGHGHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGGGHGHGHHHHGHHHHHHGGGHHHHHHHHHHHHHHHHHHHHHHHGGGGGHHHHHHHHHGHHFHHHHHHHGFFHHHHHHHHHHGHHGHG

original fastq read 2:

@M01950:199:000000000-CN4MR:1:1106:27160:16985 2:N:0:2 TTGAACTCAGTAGCTCACAAAACATCCACTCTGCCTCCAAAGGCCTACCTCTGAACCATGCTGCCTTGCCTCCTACAGACTATGACAGAAGTCCCTTTGTAACATCCCCCATTAGCATGACAATGCCCCCTCACGGCAGCCTGCAAGGTT + CCCCCFFFFFFFGGGGGGGGGGHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGHHHHHHHHHHHHHHHHHHGHHHHHHHHHHHHHHGHHHGFGHHHHHHHHGGGGHHHHHGHHHHHHGGHHGGGGHHHGGGGGHHHHHGHGHH

bam-file reads after UMI extraction:

M01950:199:000000000-CN4MR:1:1106:8824:16985+CACAGGTTGAAC 163 chr1 200048348 60 133M = 200048464 250 CAAAACATCCACTCTGCCTCCAAAGGCCTACCTCTGAACCATGCTGCCTTGCCTCCTACAGACTATGACAGAAGTCCCTTTGTAACATCCCCCATTAGCATGACAATGCCCCCTCACGGCAGCCTGCAAGGTT GGGGGHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGHHHHHHHHHHHHHHHHHHGHHHHHHHHHHHHHHGHHHGFGHHHHHHHHGGGGHHHHHGHHHHHHGGHHGGGGHHHGGGGGHHHHHGHGHH ZA:Z:CACAGG ZB:Z:TTGAAC MC:Z:134M MD:Z:133 RG:Z:A NM:i:0 MQ:i:60 UQ:i:0 AS:i:133RX:Z:CACAGG-TTGAAC

M01950:199:000000000-CN4MR:1:1106:8824:16985+CACAGGTTGAAC 83 chr1 200048464 60 134M = 200048348 -250 CGGCAGCCTGCAAGGTTACCAAACATATGGCCACTTTCCTAGCCGGGCCATCAAGTCTGAGTACCCAGACCCCTATACCAGCTCACCCGAGTCCATAATGGGCTATTCATATATGGATAGTTACCAGACGAGCT GHGHHGHHHHHHHHHHFFGHHHHHHHFHHGHHHHHHHHHGGGGGHHHHHHHHHHHHHHHHHHHHHHHGGGHHHHHHGHHHHGHGHGGGHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGHGGGGGGGG ZA:Z:CACAGG ZB:Z:TTGAAC MC:Z:133M MD:Z:134 RG:Z:A NM:i:0 MQ:i:60 UQ:i:0 AS:i:134RX:Z:CACAGG-TTGAAC

grouped reads:

M01950:199:000000000-CN4MR:1:1106:8824:16985+CACAGGTTGAAC 83 chr1 200048464 60 134M = 200048348 -250 CGGCAGCCTGCAAGGTTACCAAACATATGGCCACTTTCCTAGCCGGGCCATCAAGTCTGAGTACCCAGACCCCTATACCAGCTCACCCGAGTCCATAATGGGCTATTCATATATGGATAGTTACCAGACGAGCT GHGHHGHHHHHHHHHHFFGHHHHHHHFHHGHHHHHHHHHGGGGGHHHHHHHHHHHHHHHHHHHHHHHGGGHHHHHHGHHHHGHGHGGGHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGHGGGGGGGG ZA:Z:CACAGG ZB:Z:TTGAAC MC:Z:133M MD:Z:134 RG:Z:A MI:Z:6968/B NM:i:0 MQ:i:60 UQ:i:0 AS:i:134 RX:Z:CACAGG-TTGAAC

M01950:199:000000000-CN4MR:1:1106:8824:16985+CACAGGTTGAAC 163 chr1 200048348 60 133M = 200048464 250 CAAAACATCCACTCTGCCTCCAAAGGCCTACCTCTGAACCATGCTGCCTTGCCTCCTACAGACTATGACAGAAGTCCCTTTGTAACATCCCCCATTAGCATGACAATGCCCCCTCACGGCAGCCTGCAAGGTT GGGGGHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHHGHHHHHHHHHHHHHHHHHHGHHHHHHHHHHHHHHGHHHGFGHHHHHHHHGGGGHHHHHGHHHHHHGGHHGGGGHHHGGGGGHHHHHGHGHH ZA:Z:CACAGG ZB:Z:TTGAAC MC:Z:134M MD:Z:133 RG:Z:A MI:Z:6968/B NM:i:0 MQ:i:60 UQ:i:0 AS:i:133 RX:Z:CACAGG-TTGAAC

For this exemplary MI-tag, I don't get any /A reads at all. This is the case for 99% of the reads, where there are either only /A or only /B, but not both at the same time. I think therefore I'm not able to call Duplex Consensus reads. After going through the entire data, I'm still not able to figure out where the issue is located.

I would appreciate any help solving this!
Thank you in advance
Adrian

[nh13]

I can reproduce with:

java -jar fgbio-2.0.3-da9ecbcc-SNAPSHOT.jar FastqToBam -i fastq1.fq fastq2.fq \
    --sample sample --library sample -r "6M11S+T" "6M11S+T" -o unmapped.bam

samtools fastq -N unmapped.bam \
    | bwa mem -Y -p ~hs38DH.fa - \
    | java -jar fgbio-2.0.3-da9ecbcc-SNAPSHOT.jar ZipperBams -u unmapped.bam -r hs38DH.fa \
    > mapped.bam

java -jar fgbio-2.0.3-da9ecbcc-SNAPSHOT.jar GroupReadsByUmi -i mapped.bam -o grouped.sam \
    --strategy=paired --edits=1

[ARE4S]

Thanks a lot for your quick reply!

I've tested the data using your commands with same results.
But Duplex Consensus calling failed due to Not Enough Reads:

[2022/10/27 09:36:50 | CallDuplexConsensusReads | Info] Total Raw Reads Considered: 2.331.558.
[2022/10/27 09:36:50 | CallDuplexConsensusReads | Info] Raw Reads Filtered Due to Being Fragment/Non-Paired Reads: 0 (0,0000).
[2022/10/27 09:36:50 | CallDuplexConsensusReads | Info] Raw Reads Filtered Due to Not Enough Reads (Either Total, AB, or BA): 2.331.374 (0,9999).
[2022/10/27 09:36:50 | CallDuplexConsensusReads | Info] Raw Reads Filtered Due to Mismatching Cigars: 1 (0,0000).
[2022/10/27 09:36:50 | CallDuplexConsensusReads | Info] Consensus reads emitted: 12.
[2022/10/27 09:36:50 | FgBioMain | Info] CallDuplexConsensusReads completed. Elapsed time: 0,44 minutes.

Am I missing something basic?

[nh13]

Nothing obvious that I can see. Please be patient as it may a little while before we can volunteer time to look into this as we’re prioritizing client work.

[ARE4S]

This is understandable. Thank you for your effort!

[BoudSTER]

Hello,

I have the same issue ...

[2023/02/27 15:32:14 | FgBioMain | Info] Executing CallDuplexConsensusReads from fgbio version 1.0.0 as maxime.bucci@hpc-ucs13 on JRE 11.0.1-internal+0-adhoc..src with snappy, JdkInflater, and JdkDeflater
[2023/02/27 15:32:42 | CallDuplexConsensusReads | Info] processed     1,000,000 records.  Elapsed time: 00:00:28s.  Time for last 1,000,000:   28s.  Last read position: chr2:198,265,400
[2023/02/27 15:33:09 | CallDuplexConsensusReads | Info] processed     2,000,000 records.  Elapsed time: 00:00:54s.  Time for last 1,000,000:   26s.  Last read position: chr4:153,249,495
[2023/02/27 15:33:33 | CallDuplexConsensusReads | Info] processed     3,000,000 records.  Elapsed time: 00:01:19s.  Time for last 1,000,000:   24s.  Last read position: chr7:92,761,629
[2023/02/27 15:34:00 | CallDuplexConsensusReads | Info] processed     4,000,000 records.  Elapsed time: 00:01:46s.  Time for last 1,000,000:   26s.  Last read position: chr9:139,401,633
[2023/02/27 15:34:29 | CallDuplexConsensusReads | Info] processed     5,000,000 records.  Elapsed time: 00:02:15s.  Time for last 1,000,000:   29s.  Last read position: chr12:6,128,523
[2023/02/27 15:35:02 | CallDuplexConsensusReads | Info] processed     6,000,000 records.  Elapsed time: 00:02:47s.  Time for last 1,000,000:   32s.  Last read position: chr14:96,180,287
[2023/02/27 15:35:34 | CallDuplexConsensusReads | Info] processed     7,000,000 records.  Elapsed time: 00:03:19s.  Time for last 1,000,000:   31s.  Last read position: chr17:29,663,664
[2023/02/27 15:36:06 | CallDuplexConsensusReads | Info] processed     8,000,000 records.  Elapsed time: 00:03:52s.  Time for last 1,000,000:   32s.  Last read position: chr19:17,954,628
[2023/02/27 15:36:35 | CallDuplexConsensusReads | Info] processed     9,000,000 records.  Elapsed time: 00:04:20s.  Time for last 1,000,000:   28s.  Last read position: chrX:44,733,151
[2023/02/27 15:36:47 | CallDuplexConsensusReads | Info] processed     9,529,476 records.  Elapsed time: 00:04:33s.  Time for last 529,476:   12s.  Last read position: hs37d5:35,474,270
[2023/02/27 15:36:47 | CallDuplexConsensusReads | Info] Total Raw Reads Considered: 9,529,476.
[2023/02/27 15:36:47 | CallDuplexConsensusReads | Info] Raw Reads Filtered Due to Mismatching Cigars: 8,063 (0.0008).
[2023/02/27 15:36:47 | CallDuplexConsensusReads | Info] Raw Reads Filtered Due to Not Enough Reads (Either Total, AB, or BA): 7,923,946 (0.8315).
[2023/02/27 15:36:47 | CallDuplexConsensusReads | Info] Raw Reads Filtered Due to Being Fragment/Non-Paired Reads: 0 (0.0000).
[2023/02/27 15:36:47 | CallDuplexConsensusReads | Info] Consensus reads emitted: 591,646.
[2023/02/27 15:36:47 | FgBioMain | Info] CallDuplexConsensusReads completed. Elapsed time: 4.61 minutes.
[bam_sort_core] merging from 0 files and 2 in-memory blocks...

Do you have a solution ?
Thanks in advance.

Boudster

[nh13]

Sorry, we have had time to take a look. Feel free to reach out to us if you'd like to sponsor this work.

[BoudSTER]

Hi :),

On a new run the problem persists, do you have any clues?

[2023/07/03 20:12:58 | CallDuplexConsensusReads | Info] Total Raw Reads Considered: 200,702,824.
[2023/07/03 20:12:58 | CallDuplexConsensusReads | Info] Raw Reads Filtered Due to Mismatching Cigars: 463,122 (0.0023).
[2023/07/03 20:12:58 | CallDuplexConsensusReads | Info] Raw Reads Filtered Due to Not Enough Reads (Either Total, AB, or BA): 158,664,354 (0.7905).
[2023/07/03 20:12:58 | CallDuplexConsensusReads | Info] Raw Reads Filtered Due to Being Fragment/Non-Paired Reads: 0 (0.0000).
[2023/07/03 20:12:58 | CallDuplexConsensusReads | Info] Consensus reads emitted: 1,639,244.
[2023/07/03 20:12:58 | FgBioMain | Info] CallDuplexConsensusReads completed. Elapsed time: 21.54 minutes.

Thanks in advance
And thank you very much for this valuable tool

BoudSTER

[nh13]

Sorry we haven’t looked into this. We provide this software and support on a volunteer basis so I don’t know when we’ll look into this. Have you had a chance to obtain funds for sponsorship?