SEQanswers

Go Back   SEQanswers > Bioinformatics > Bioinformatics



Similar Threads
Thread Thread Starter Forum Replies Last Post
Sureselect disappearing DNA JPC General 3 10-24-2016 03:25 AM
DNA concentration methods, SureSelect Squall Sample Prep / Library Generation 2 04-06-2016 12:27 AM
Agilent SureSelect DNA and Haloplex woodydon RNA Sequencing 2 03-16-2014 04:01 AM
low sureselect on-target reads hawaii454-0 Sample Prep / Library Generation 3 04-10-2013 10:18 AM
question on DNA fragment size after shear the DNA (in Chip-seq) kaixinsjtu Sample Prep / Library Generation 4 04-05-2012 04:36 AM

Reply
 
Thread Tools
Old 10-04-2021, 05:12 AM   #1
gcarbajosa
Junior Member
 
Location: London

Join Date: May 2009
Posts: 5
Default 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
gcarbajosa is offline   Reply With Quote
Reply

Thread Tools

Posting Rules
You may not post new threads
You may not post replies
You may not post attachments
You may not edit your posts

BB code is On
Smilies are On
[IMG] code is On
HTML code is Off




All times are GMT -8. The time now is 02:14 PM.


Powered by vBulletin® Version 3.8.9
Copyright ©2000 - 2021, vBulletin Solutions, Inc.
Single Sign On provided by vBSSO