分析带UMI标签的测序数据

检测癌组织的低频突变，为了提高检测低频突变的灵敏度，往往进行高深度的测序。但样本之间存在交叉污染，测序有存在一定概率的错误，这些因素会导致高深度测序过程中将假阳性的信号放到，得到假阳性的结果。解决交叉污染的方法，有公司比如IDT采用唯一配对的样本index，只有配对的index中的reads才属于特定样本。解决测序错误的方法，研究人员在建库的时候，先对分子加上UMI碱基，unique molecular identifier -> UMI，然后根据来源于同一个分子的测序数据进行测序错误修正，得到正确的分子序列。两种方法结合可以减少交叉污染提高准确性（https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5759201/）。

如图中所示，左侧一个分子被测了5次，其中第二次有一个测序错误，但该错误并没有在每个测序数据中出现，所以在后续合成一个分子的时候，测序错误被修正，只保留了真正的突变。（https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5852328/）

常规的肿瘤配对测序分析，或者遗传性突变位点的分析，并不需要UMI信息，所以包含UMI的数据分析是需要不一样的分析流程来得到准确的分析结果，其中包括提取UMI分子标签，合并来自同一个分子的测序reads，低频突变检测而非胚系突变检测等。

大致流程为：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24


-----prepare analysis ready BAM file------
'         FASTQ
'            ↓
'          uBAM
'            ↓  Extract UMI
'          uBAM
'             ↓  Align uBAM and merge 
'           BAM
-----call consensus------
'             ↓  Group Reads By Umi
'           BAM
'             ↓  Group Reads By Umi
'           BAM
'             ↓  Call Molecular Consensus Reads
'           uBAM
'             ↓  Align uBAM and merge
'           BAM
'             ↓  Filter Consensus Reads
'           BAM
'             ↓  Clip
'           BAM
-----Vardict------
'             ↓  Call
'            VCF

一，得到包含UMI分子标签信息的BAM文件

UMI信息，应该从fastq中的配对的reads中提取，但fastq不能存储更多的信息，所以需要先将fastq转成uBAM文件，提取uBAM文件中的UMI分子标签信息，将该信息通过RX标签写入uBAM文件中，通过uBAM和BAM文件合并，把RX信息合并到比对到的BAM文件中进行下一步分析。

1）生成uBAM

1
2
3
4
5
6
7
8
9


java -Xmx8G -jar picard.jar FastqToSam \
    FASTQ=$fq1    \
    FASTQ2=$fq2   \
    OUTPUT=test.uBAM \
    READ_GROUP_NAME=test \
    SAMPLE_NAME=test \
    LIBRARY_NAME=test \
    PLATFORM_UNIT=HiseqX10  PLATFORM=illumina \
    RUN_DATE=`date --iso-8601=seconds`

2）提取UMI信息

1
2
3
4
5
6


java -jar fgbio.jar ExtractUmisFromBAM \
    --input=test.uBAM  \
    --output=test.umi.uBAM \
    --read-structure=2M148T 2M148T \
    --single-tag=RX \
    --molecular-index-tags=ZA ZB

此时，uBAM文件中RX标签记录着UMI的信息，2M148T表示前两个碱基是UMI分子标签，148个template 测序序列，配对read在uBAM文件中有如下信息

1
2


R1----ZA:Z:TC ZB:Z:TA RG:Z:test       RX:Z:TC-TA
R2----ZA:Z:TC ZB:Z:TA RG:Z:test       RX:Z:TC-TA

3）比对uBAM文件中的reads

1
2


samtools fastq  test.umi.uBAM ' bwa mem -t 8 -p reference.fa /dev/stdin ' samtools view -b
 > test.umi.BAM

4）uBAM和BAM合并

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


java -Xmx8G -jar picard.jar MergeBAMAlignment R=reference.fa \
    UNMAPPED_BAM=test.umi.uBAM  \
    ALIGNED_BAM=test.umi.BAM \
    O=test.umi.merged.BAM  \
    CREATE_INDEX=true    \
    MAX_GAPS=-1 \
    ALIGNER_PROPER_PAIR_FLAGS=true \
    VALIDATION_STRINGENCY=SILENT \
    SO=coordinate \
    ATTRIBUTES_TO_RETAIN=XS

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

二，Call Consensus Reads

1）Group Reads By Umi

该步会生成一个包含MI tag的文件，存储每个read的最初分子的ID，低质量的read应该过滤掉，因为低比对质量的read上mismatch较多，容易造成假阳性。

1
2
3
4


java -jar fgbio.jar GroupReadsByUmi \
    --input=test.umi.merged.BAM \
    --output=test.umi.group.BAM  \
    --strategy=paired  --min-map-q=20  --edits=1 --raw-tag=RX

2）Call Molecular Consensus Reads

根据UMI分子标签的信息和Read1、Read2的位置，从一组read中识别出最初的 molecular 分子序列。min-reads是必填的，如果测序深度较低，单个read也可以用来call consensus，此时设置1，如果深度较高可以设置2或者更高，但1是最安全的，因为后续可以基于此参数进一步过滤。但同时该值的设置显著影响效率，此时我们设置1。此时生成的文件是uBAM格式的，需要进一步比对

1
2
3
4
5


java -jar fgbio.jar  CallMolecularConsensusReads \
    --min-reads=1 \
    --min-input-base-quality=20 \
    --input=test.umi.group.BAM \
    --output=test.consensus.uBAM

3）比对uBAM文件中的reads

1

samtools fastq $dir/$smp.callconsensus.BAM ' bwa mem -t 8 -p reference.fa  /dev/stdin ' samtools view -b - ' test.consensus.BAM

4）uBAM和BAM合并

 1
 2
 3
 4
 5
 6
 7
 8
 9
10


java -Xmx8G -jar picard.jar MergeBAMAlignment R=reference.fa \
    UNMAPPED_BAM=test.consensus.uBAM  \
    ALIGNED_BAM=test.consensus.BAM \
    O=test.consensus.merge.BAM  \
    CREATE_INDEX=true    \
    MAX_GAPS=-1 \
    ALIGNER_PROPER_PAIR_FLAGS=true \
    VALIDATION_STRINGENCY=SILENT \
    SO=coordinate \
    ATTRIBUTES_TO_RETAIN=XS

4）Filter Consensus Reads

1
2
3
4
5
6
7
8


java -jar fgbio.jar FilterConsensusReads \
    --input=test.consensus.merge.BAM  \
    --output=test.consensus.merge.filter.BAM \
    --ref=reference --min-reads=2 \
    --max-read-error-rate=0.05 \
    --max-base-error-rate=0.1 \
    --min-base-quality=30 \
    --max-no-call-fraction=0.20

5）Clip

来源于同一个template的read，若果有重叠，重叠部分的突变其实应该只计一次，所以要clip一下read。

1
2
3
4


java -jar fgbio.jar ClipBAM \
    --input=test.consensus.merge.filter.BAM   \
    --output=test.consensus.merge.filter.clip.BAM \
    --ref=$ref  --soft-clip=false --clip-overlapping-reads=true

三、Variant Call

1
2
3
4
5
6


AF_THR="0.01"
VarDict/bin/VarDict -G reference.fa -f $AF_THR -N test \
    -b test.consensus.merge.filter.BAM \
    -z -c 1 -S 2 -E 3 -g 4 -th 4 target.bed ' \
    VarDictJava/VarDict/teststrandbias.R ' \
    VarDictJava/VarDict/var2vcf_valid.pl -N test -E -f $AF_THR > test.vcf