循环肿瘤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 4 
  SampleID_probes.bed > S1.vcf