Seqanswers Leaderboard Ad

Collapse

Announcement

Collapse
No announcement yet.
X
 
  • Filter
  • Time
  • Show
Clear All
new posts

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

    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

Latest Articles

Collapse

  • seqadmin
    Exploring the Dynamics of the Tumor Microenvironment
    by seqadmin




    The complexity of cancer is clearly demonstrated in the diverse ecosystem of the tumor microenvironment (TME). The TME is made up of numerous cell types and its development begins with the changes that happen during oncogenesis. “Genomic mutations, copy number changes, epigenetic alterations, and alternative gene expression occur to varying degrees within the affected tumor cells,” explained Andrea O’Hara, Ph.D., Strategic Technical Specialist at Azenta. “As...
    07-08-2024, 03:19 PM
  • seqadmin
    Exploring Human Diversity Through Large-Scale Omics
    by seqadmin


    In 2003, researchers from the Human Genome Project (HGP) announced the most comprehensive genome to date1. Although the genome wasn’t fully completed until nearly 20 years later2, numerous large-scale projects, such as the International HapMap Project and 1000 Genomes Project, continued the HGP's work, capturing extensive variation and genomic diversity within humans. Recently, newer initiatives have significantly increased in scale and expanded beyond genomics, offering a more detailed...
    06-25-2024, 06:43 AM

ad_right_rmr

Collapse

News

Collapse

Topics Statistics Last Post
Started by seqadmin, Today, 11:09 AM
0 responses
14 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-19-2024, 07:20 AM
0 responses
146 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-16-2024, 05:49 AM
0 responses
120 views
0 likes
Last Post seqadmin  
Started by seqadmin, 07-15-2024, 06:53 AM
0 responses
111 views
0 likes
Last Post seqadmin  
Working...
X