• 周五. 12月 9th, 2022

5G编程聚合网

5G时代下一个聚合的编程学习网

热门标签

UMI 的原理分析带有 UMI 的数据

admin

11月 28, 2021

循环肿瘤DNA,即 ctDNA(circulating tumor DNA),是指肿瘤细胞经外泌体、细胞坏死或凋亡释放进入循环系统中的DNA。ctDNA 可作为肿瘤生物标记物,我们通过检测血液中的 ctDNA,就能够追踪到肿瘤在血液中的踪迹。

目前ctDNA的分析方法主要有三种:以 ARMS 为代表的定量 PCR 技术;以 ddPCR 为代表的数字 PCR 技术;基于下一代测序 (NGS) 的高通量检测技术。NGS 能同时检测多个基因的多种不同变异形式,是应用最广泛的的基因检测技术。然而,基于 NGS 的 ctDNA 分析方法还存在两个技术障碍:

  • 血浆中的ctDNA含量极低;

  • 二代测序背景噪声很高,用常规建库测序方法肿瘤信号将会完全淹没在背景噪声中。

所以为了在保证检测准确性和测序深度的基础上降低价格,目前越来越多的新技术先对ctDNA目的片段进行富集,再将富集后的目的序列进行深度测序。此外,分子标签技术(Unique Molecular indentifier,UMI)也被广泛应用以提高检测的灵敏度和特异性。它的原理就是给每一条原始DNA片段加上一段特有的标签序列,经文库构建及PCR扩增后一起进行测序。这样,根据不同的标签序列我们就可以区分不同来源的DNA模板,分辨哪些是PCR扩增及测序过程中的随机错误造成的假阳性突变,哪些是患者真正的携带的突变。

UMI 的原理

如图中所示,α 和 β 均是连接到 DNA 模板上的随机序列,不同的模板连接的序列不同,即为 ctDNA 加上了独特的标记。经 PCR 扩增后通过 α 和 β 的序列就可以识别出同一模板来源的片段,再根据 α 和 β 互补识别同一模板的正反链。

  • i 图中一个错误仅影响了扩增出的少数序列,统一分析同一模板来源的片段可以将这个错误过滤;

  • ii 图中出现的错误影响了所有正链序列,但反链正常,根据碱基互补原则,可将此种错误校正;

  • iii 图中则显示了排除了一些引入错误和系统错误后得到真正的突变。

分析带有 UMI 的数据

所需软件

  • Picard:https://github.com/broadinstitute/picard

  • bwa:https://github.com/lh3/bwa

  • fgbio:https://github.com/fulcrumgenomics/fgbio

  • VarDict:https://github.com/AstraZeneca-NGS/VarDictJava

流程

1. 得到带有 UMI 标签信息的 BAM 文件

先用picard FastqToSam将质控后的fastq数据转化为bam格式:

java -Xmx8G -jar /usr/local/bin/picard.jar FastqToSam 
F1=/data1/zwbao/ctDNA/sample/s1-1.fq.gz 
F2=/data1/zwbao/ctDNA/sample/s1-2.fq.gz 
O=S1.unmapped.bam 
SAMPLE_NAME=S1pwd

bam文件中提取 UMI 标签信息的bam文件,并生成带有 UMI 标签信息的未比对的bam文件:

java -Xmx8G -jar   /data/zwbao/biosoft/ctDNA/fgbio/target/scala-2.12/fgbio-0.7.0-69c80dd-SNAPSHOT.jar ExtractUmisFromBam 
--input=S1.unmapped.bam --output=S1.unmapped.withUMI.bam 
--read-structure=3M2S146T 3M2S146T --molecular-index-tags=ZA ZB 
--single-tag=RX

bam文件中 RX 标签记录着 UMI 的信息,3M2S146T 表示前三个碱基是UMI分子标签,接着两个碱基需要被跳过(因为我的数据还插入了两个随机碱基),146个碱基是测序序列。

将未比对的bam文件转为fastq格式,并进行比对,得到比对后的bam文件:

java -Xmx8G -jar   /usr/local/bin/picard.jar SamToFastq I=S1.unmapped.withUMI.bam F=S1.SamToFastq INTERLEAVE=true

bwa mem -p -t 35 /data1/zwbao/shared/hg19/hg19.fa S1.SamToFastq > S1.mapped.bam

通过合并,得到一个包含各种信息包括 RX 标签等的bam文件,该文件用于下一步 call consensus read,因为会将来源于同一个分子的read合并成一个consensus read,所以在得到bam文件之后,没有进行mark duplication

java -Xmx8G -jar   /usr/local/bin/picard.jar  MergeBamAlignment UNMAPPED=S1.unmapped.withUMI.bam ALIGNED=S1.mapped.bam O=S1.mapped.withUMI.bam R=/data1/zwbao/shared/hg19/hg19.fa SO=coordinate ALIGNER_PROPER_PAIR_FLAGS=true MAX_GAPS=-1 ORIENTATIONS=FR VALIDATION_STRINGENCY=SILENT CREATE_INDEX=true

2. Map BAM to consensus reads

该步会生成一个包含 UMI 的文件,存储每个 read 的 ID,低质量的 read 应被过滤:

java -Xmx8G -jar   /data/zwbao/softwares/ctDNA/fgbio/target/scala-2.12/fgbio-0.7.0-69c80dd-SNAPSHOT.jar GroupReadsByUmi 
  --input=S1.mapped.withUMI.bam --output=S1.GroupedReads.bam 
  --strategy=paired --edits=1 --min-map-q=20

此步操作通过识别 UMI,找到来自同一片段的数据
根据 UMI 分子标签的信息和 Read1、Read2 的位置,从一组read中识别出最初的 molecular 分子序列。此时生成的文件是未比对的bam文件,需要进一步比对。

java -Djava.io.tmpdir=tmpDir -Xmx8G -jar /data/zwbao/softwares/ctDNA/fgbio/target/scala-2.12/fgbio-0.7.0-69c80dd-SNAPSHOT.jar CallDuplexConsensusReads 
  --input=S1.GroupedReads.bam --output=S1.consensus.unmapped.bam 
  --error-rate-pre-umi=45 --error-rate-post-umi=30 
  --min-input-base-quality=30

对得到的consensus reads进行过滤:

java -Djava.io.tmpdir=tmpDir -Xmx8G -jar /data/zwbao/softwares/ctDNA/fgbio/target/scala-2.12/fgbio-0.7.0-69c80dd-SNAPSHOT.jar FilterConsensusReads 
  --input=S1.consensus.unmapped.bam 
  --output=S1.consensus.filtered.unmapped.bam 
  --ref=/data1/zwbao/shared/hg19/hg19.fa 
  --min-reads=2 1 1 
  --max-read-error-rate=0.05 
  --max-base-error-rate=0.1 
  --min-base-quality=50 
  --max-no-call-fraction=0.05

进行比对,合并bam文件:

java -Xmx8G -jar   /usr/local/bin/picard.jar SamToFastq I=S1.consensus.filtered.unmapped.bam F=S1.consensus.filtered.unmapped.Fastq INTERLEAVE=true

bwa mem -p -t 30   /data1/zwbao/shared/hg19/hg19.fa  S1.consensus.filtered.unmapped.Fastq  >  S1.consensus.filtered.mapped.bam

java -Xmx8G -jar   /usr/local/bin/SortSam.jar  INPUT=S1.consensus.filtered.mapped.bam   OUTPUT=S1.consensus.filtered.mapped.sorted.bam  SORT_ORDER=queryname

java -Xmx8G -jar   /usr/local/bin/SortSam.jar  INPUT=S1.consensus.filtered.unmapped.bam   OUTPUT=S1.consensus.filtered.unmapped.sorted.bam  SORT_ORDER=queryname

java -Xmx8G -jar   /usr/local/bin/picard.jar  MergeBamAlignment UNMAPPED=S1.consensus.filtered.unmapped.sorted.bam  ALIGNED=S1.consensus.filtered.mapped.sorted.bam O=S1.consensus.filtered.mapped.withUMI.bam R=/data1/zwbao/shared/hg19/hg19.fa SO=coordinate ALIGNER_PROPER_PAIR_FLAGS=true MAX_GAPS=-1 ORIENTATIONS=FR VALIDATION_STRINGENCY=SILENT CREATE_INDEX=true

来源于同一段序列的read,若果有重叠,重叠部分的突变其实应该只计一次,所以要clip一下read

java -Djava.io.tmpdir=tmpDir -Xmx8G -jar /data/zwbao/softwares/ctDNA/fgbio/target/scala-2.12/fgbio-0.7.0-69c80dd-SNAPSHOT.jar ClipBam 
  --input=S1.consensus.filtered.mapped.withUMI.bam --output=S1.consensus.filtered.mapped.withUMI.clipped.bam 
  --ref=/data1/zwbao/shared/hg19/hg19.fa --soft-clip=false --clip-overlapping-reads=true

3. Produce variant calls from consensus reads

/data/zwbao/softwares/ctDNA/VarDictJava/build/install/VarDict/bin/VarDict 
  -G /data1/zwbao/shared/hg19/hg19.fa 
  -N S1 
  -f 0.01  # 0.02
  -b S1.consensus.filtered.mapped.withUMI.clipped.bam 
  -z -c 1 -S 2 -E 3 -g 4 -th 
  SampleID_probes.bed > S1.vcf

发表回复

您的电子邮箱地址不会被公开。 必填项已用*标注