How to group reads from same duplex DNA fragment (Sureselect XT HS) grouping by UMI

[gcarbajosa]

Hello everybody,

I am trying to analyse a dataset generated using Sureselect XT-HS.

I have three fastq files: forward, reverse and UMIs.

For example:

# FORWARD
(base) [k1506263@login3(rosalind) HUMfeqzR]$ zcat 001MS/210601_I100400190019_V350004379_L1_001MS_1.fq.gz |head
@V350004379L1C001R0010000005:0:0:0:0 1:N:0:TCTCTTCA
TCCCTTTCTCTCCTTGGTACTTCTCTGCCAAATCCCTCCCAGAGAAACAGAGAACTTCTTTCTCTACCTGAACCAGGTTTTTAGGAAATGGGCCTTTTTC
+
GG?G?F@FGGGFGFB;FGGGGEAEGGEFG7GF>EGFD2DBFGCGBGGGAEBF?@?FGGGFFGFFDGFGEFFE@GDCF?GGFGEFG<EGEGFECFF76=AG
@V350004379L1C001R0010000006:0:0:0:0 1:N:0:TCTCTTCA
CCTACCAGACCCTTCAGTGAAGCTTTGATGATGGGCTTACTGACCAACCTGGCAGACAGGGAGCTGGTTCACATGATCAACTGGGCGAAGAGGGTGCCAG
+
FGE<FFFB1F;G=@FDFGD-F8@GB(F>>9F?CFCGFDFGCFD9EFGGEC6FBEB2*EEF(;7FGEF@F?EAEAF8CFB4?FD>FC9FFD7EFF?FF@F@
@V350004379L1C001R0010000028:0:0:0:0 1:N:0:TCTCTTCA
CGCCAATTTTAGCCTCTTTGCCCCCTTTGGCAGTGGCTGCCCAGTCAATGGGCAGAGCCAGCTGAGGGGGGCCTTTGGAAGTGGGGCGCTGCCCACTGGC

# REVERSE
(base) [k1506263@login3(rosalind) HUMfeqzR]$ zcat 001MS/210601_I100400190019_V350004379_L1_001MS_2.fq.gz | head
@V350004379L1C001R0010000005:0:0:0:0 2:N:0:TCTCTTCA
GGCCTCAGAAGACAGAAGAGAGACTGGGTTATTCCTCCCATCAGCTGCCCAGAAAATGAAAAAGGCCCATTTCCTAAAAACCCGGGTCAGGGAGAGAAAG
+
EDDDFD;B02/8F;F@DE:FFFEA<CAFDF?FFCF>FFF9FFA89@C=FF7CFFF77>-F4;FEF:FFFCBF,F93CA4FEC'97(G9CB;49*-F'=E9
@V350004379L1C001R0010000006:0:0:0:0 2:N:0:TCTCTTCA
ACTATTTCTCCCATGACATCACAACAAGTTCTTGAAAAGCTATTGACTCTTAAAAGCAGCGCTTCGCATTCTTACCTGGCACCCTCTTCGCCCAGGTGAT
+
FDAEFDFFFFFF;8FEED;CGAFDFDDFFFFDFD<F7AFFFEF=F@?FF=>AFF>FA'F>?FEA<?ECFF?CGDE@/2AD,F@C@FCF5?FF6EE.;FF/
@V350004379L1C001R0010000028:0:0:0:0 2:N:0:TCTCTTCA
CCCAAAAAATATTGCCATTTTTCTGACCTTGGTAAGCAGCTGGGAATAGTAGGCAGGGCCAGTGGGCAGCGCCCCACTTCCAAAGGCCCCCCTCAGCTGG

# UMIS
(base) [k1506263@login3(rosalind) HUMfeqzR]$ zcat UMI/210601_I100400190019_V350004379_L1_CHKPE85221050164-0/001MS_out2.fq.gz | head
@V350004379L1C001R0010000005:0:0:0:0 2:N:0:TCTCTTCA
TACAGGGGAG
+
)EEE>CFC9F
@V350004379L1C001R0010000006:0:0:0:0 2:N:0:TCTCTTCA
GTATAACCCT
+
F>EEFFEFFD
@V350004379L1C001R0010000028:0:0:0:0 2:N:0:TCTCTTCA

I believe that the Sureselect protocol produces reads from "duplex" fragments where you can identify pairs of reads coming from opposite DNA strands of the same DNA fragment and use that information to group the mapped reads. But I can't seem to find how that information is encoded within the UMI tags. In other words, how can which reads come from the same DNA fragment but opposite strands since they have different UMI tags?

Below is the pipeline I was using to process the data. I am using gatk (Picard) and fgbio and I was ignoring any "duplex" DNA information as I do not know how to identify or group reads that "come" from the same DNA fragment. In particular, I was using GroupReadsByUmi with the "adjacency" option, not paired as I didn't have UMI tags in the format described in fgbio docs (e.g. ACGT-CCGG) and CallMolecularConsensusReads (instead of CallDuplexConsensusReads).

Anybody can give any clues on how to approach this? Thanks in advance for you help.

##Generate read alignments andUMIs
#"Read trimming".
# e.g.
trim_galore --fastqc --paired --cores 4 --gzip --output_dir /scratch/groups/hodgkinsonlab/gcarbajo/projects/ctDNA/analysis/trimgalore/001MS /scratch/groups/hodgkinsonlab/gcarbajo/projects/ctDNA/data/ctDNA_IBIMA/HUMfeqzR/001MS/*_1.fq.gz /scratch/groups/hodgkinsonlab/gcarbajo/projects/ctDNA/data/ctDNA_IBIMA/HUMfeqzR/001MS/*_2.fq.gz
# Filter the UMI file to keep only the UMIs for the reads that passed the trimming
# Custom perl script that keeps the IDs for reads that pass the trimming
# Generate and "unmapped bam" with Fastq2Bam. Based on (https://github.com/fulcrumgenomics/f...ead-Structures ) to codify the barcode (10M)

e.g.
fgbio FastqToBam -i ../analysis/trimgalore/001MS/*val_1.fq.gz ../analysis/trimgalore/001MS/*val_2.fq.gz 001MS.trimmed.UMIs.fq.gz -o /dev/stdout --read-structures +T +T 10M --sample=001MS --library=001MS | samtools sort -n - |samtools view -Sb - > ../analysis/unmappedBamWithUMIs/001MS.unmapWUMI.srt.v3.bam
# Map the reads with BWA ( piping amToFastq, bwa y MergeBamAlignment)
# e.g.
# Map it and merge with the UMIS
gatk --java-options -Xmx32g SamToFastq -I ../analysis/unmappedBamWithUMIs/001MS.unmapWUMI.srt.v3.bam -F /dev/stdout -INTERLEAVE true | bwa mem -p -t 8 /mnt/lustre/groups/hodgkinsonlab/gcarbajo/annotation/hsapiens/hg38/GCA_000001405.15/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta.gz /dev/stdin | samtools sort -n - | samtools view -Sb - | gatk --java-options -Xmx32g MergeBamAlignment -UNMAPPED ../analysis/unmappedBamWithUMIs/001MS.unmapWUMI.srt.v3.bam -ALIGNED /dev/stdin -O ../analysis/mappedBamWithUMIs/001MS.mappedWUMI.v3.bam -R /mnt/lustre/groups/hodgkinsonlab/gcarbajo/annotation/hsapiens/hg38/GCA_000001405.15/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta.gz -SO coordinate -ALIGNER_PROPER_PAIR_FLAGS true -MAX_GAPS -1 -ORIENTATIONS FR -VALIDATION_STRINGENCY SILENT -CREATE_INDEX true

# Genrate consensus reads using "Adjacency" (details here: https://genome.cshlp.org/content/27/3/491.full)

#e.g.
fgbio GroupReadsByUmi -s Adjacency -i ../analysis/mappedBamWithUMIs/001MS.mappedWUMI.v3.bam -o ../analysis/GroupedReadsByUmi/001MS.mappedWUMI.v3.gbu.bam -f ../analysis/GroupedReadsByUmi/001MS.mappedWUMI.v3.gbu.hist.txt > logs/001MS.GroupReadsByUmi.out 2>logs/001MS.GroupReadsByUmi.err

#Generate consensus sequences:
#e.g.
fgbio CallMolecularConsensusReads --input ../analysis/GroupedReadsByUmi/001MS.mappedWUMI.v3.gbu.bam --output ../analysis/CallMolecularConsensusReads/001MS.mappedWUMI.v3.gbu.cons_unmap.bam --min-reads=1 --reject ../analysis/CallMolecularConsensusReads/001MS.mappedWUMI.v3.gbu.cons_rej.bam > logs/001MS_CallMolecularConsensusReads.out 2>logs/001MS_CallMolecularConsensusReads.error


#Map consensus.
#e.g.
gatk --java-options -Xmx32g SamToFastq -RGT A -I ../analysis/CallMolecularConsensusReads/001MS.mappedWUMI.v3.gbu.cons_unmap.bam -F /dev/stdout -INTERLEAVE true | bwa mem -p -t 8 /mnt/lustre/groups/hodgkinsonlab/gcarbajo/annotation/hsapiens/hg38/GCA_000001405.15/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta.gz /dev/stdin | gatk SortSam -I /dev/stdin -O /dev/stdout -SO queryname | samtools view -Sb - | gatk --java-options -Xmx32g MergeBamAlignment -UNMAPPED ../analysis/CallMolecularConsensusReads/001MS.mappedWUMI.v3.gbu.cons_unmap.bam -ALIGNED /dev/stdin -O ../analysis/consensusAligned/001MS.mappedWUMI.v3.bam -R /mnt/lustre/groups/hodgkinsonlab/gcarbajo/annotation/hsapiens/hg38/GCA_000001405.15/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta.gz -SO coordinate -ALIGNER_PROPER_PAIR_FLAGS true -MAX_GAPS -1 -ORIENTATIONS FR -VALIDATION_STRINGENCY SILENT -CREATE_INDEX true
Filtra las secuencias consenso (aquí me quedo con familias de 2 o más pero también he probado de 3)
Por ejemplo:
fgbio FilterConsensusReads -M 2 -N 30 -i ../analysis/consensusAligned/001MS.mappedWUMI.v3.bam -o ../analysis/FilterConsensusReads/001MS.mappedWUMI.v3.FCons.min2.bam -r /mnt/lustre/groups/hodgkinsonlab/gcarbajo/annotation/hsapiens/hg38/GCA_000001405.15/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta.gz > logs/001MS.FilterConsensusReads.min2.out 2> logs/001MS.FilterConsensusReads.min2.err

"Call variants" con Mutect2
Cortar las secuencias en las que se solapan las reads FW y reverse para solo contar una vez por fragmento
Por ejemplo:
fgbio ClipBam -Xmx32000M --clip-overlapping-reads -i ../analysis/FilterConsensusReads/001MS.mappedWUMI.v3.FCons.min2.bam -o ../analysis/ClipBam/001MS.mappedWUMI.v3.FCons.min2.clip.bam -m ../analysis/ClipBam/001MS.mappedWUMI.v3.FCons.min2.clip.metrics.txt -r /mnt/lustre/groups/hodgkinsonlab/gcarbajo/annotation/hsapiens/hg38/GCA_000001405.15/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta.gz > logs/001MS.ClipBam.out 2> logs/001MS.ClipBam.err
Detectar variantes:
Por ejemplo:
gatk Mutect2 -R /mnt/lustre/groups/hodgkinsonlab/gcarbajo/annotation/hsapiens/hg38/GCA_000001405.15/GRCh38_no_alt_analysis_set_GCA_000001405.15.fasta -I ../analysis/ClipBam/001MS.mappedWUMI.v3.FCons.min2.clip.bam -L ../data/bed_files/list.interval_list --germline-resource ../data/gatk-best-practices/somatic-hg38_af-only-gnomad.hg38.vcf.gz --panel-of-normals ../data/gatk-best-practices/somatic-hg38_1000g_pon.hg38.vcf.gz -O ../analysis/Mutecc2/001MS.min2.clip.vcf.gz > logs/001MS.mutect.min2.clip.out 2> logs/001MS.mutect.min2.clip.errh