cfDNA（带UMI标签的fq数据--WES）的生信处理call突变流程

作者，Evil Genius

我们先来看看什么是UMI，这个UMI跟单细胞的UMI有什么不同呢？

单细胞空间测序的UMI
单细胞测序的步骤中增加了 UMI（unique molecular identifiers），UMIs 是由 4-10 个随机核苷酸组成的序列，在 mRNA 反转录后，进入到文库中，每一个 mRNA，随机连上一个 UMI，因此可以计数不同的 UMI，最终计数 mRNA 的数量。
10X genomics单细胞测序通过Barcode来标记细胞，UMI 来标记转录本，这样与参考基因比对后就可以定量细胞以及基因的数量。
Cell Ranger 进一步将外显子reads与参考转录本比对，寻找兼容性。注释到单个基因信息的reads认为是一个特异的转录本，只有注释到转录本的reads才用于UMI计数。

UMI简介

UMI（unique molecular identifer）是一种用来降低二代测序实验误差（sequencer, polymerase, DNA damage error and so on）的有效技术。建文库时在adapter处加入UMI，在生信分析时通过UMI区分真阳性突变和假阳性突变，提高检测的Sensitivity和Specificity。

UMI的原理

UMI的原理就是给每一条原始DNA片段加上一段特有的标签序列，经PCR扩增后一起进行测序。这样根据不同的标签序列，生信人员就能区分不同来源的DNA模板，分辨哪些是PCR扩增及测序过程中的随机错误造成的假阳性突变，哪些是患者真正携带的突变，从而提高检测灵敏度和特异性。

一些需要较高正确度的应用，如肿瘤稀有突变分析，常会对5%、甚至1%左右的稀有变异作检测。这时一旦出现测序错误、或者PCR偏差，便会产生大量假阳性结果，因此需要添加UMI进行校正。

一句话总结：做稀有突变检测、校正测序错误与PCR偏差，UMI“劳苦功高”

生信call变异流程

生物信息学分析流程（带UMI标签）

Trimmomatic或cutadapt过滤,一般采用Fastp也可以

去接头、测序引物、低质量碱基、短序列及高N比率的碱基。

FastQc数据质控

统计clean reads中碱基质量超过Q20及Q30的占比、reads数、GC比等。

picard进行uBAM转换

picard FastqToSam将fastq转uBAM（uBAM就是未经过比对软件比对的BAM文件）。

fgbio提取UMI信息

fgbio ExtractUmisFromBam提取uBAM中的UMI信息，生成带有UMI的uBAM, UMI信息存储在标签RX中。

samtools或picard进行fastq转换

samtools fastq或者picard SamToFastq将未比对带有UMI的BAM转换为fastq。

bwa比对到参考基因组

bwa index构建参考基因组索引、bwa mem比对得到sam文件。

samtools进行bam转换

samtools view -bS将sam转换为bam格式。

picard 合并uBAM和BAM

picard MergeBamAlignment将fgbio提取UMI信息的BAM与上一步比对所得BAM合并，得到一个包含各种信息包括RX tag等的BAM文件。

fgbio GroupReadsByUmi

fgbio GroupReadsByUmi通过UMI将reads分组

fgbio CallMolecularConsens

fgbio CallMolecularConsens根据UMI和read1、read2的位置，从一组read中识别出最初的分子序列，得到consensus.uBAM

比对前一步uBAM文件中的reads

通过samtools fastq、bwa mem、samtools view最终得到Consens BAM。

合并consensus.uBAM和Consens BAM

使用picard MergeBAMAlignment合并consensus.uBAM和Consens BAM得到consensus.merge.BAM。

fgbio FilterConsensusReads

使用fgbio FilterConsensusReads 对consensus.merge.BAM过滤得consensus.merge.filter.BAM。

fgbio ClipBAM

fgbio ClipBAM对来源于同一个template的read重叠部分突变clip，得consensus.merge.filter.clip.BAM。

call变异

我们来看看代码过程

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

picard FastqToSam F1=test_read1.fq.gz F2=test_read2.fq.gz O=test.uBam SM=testsample
 
fgbio  -Xmx50G ExtractUmisFromBam  -i  test.uBam  -o  test.umi.uBam  -r  5M2S+T 5M2S+T  -s  RX -t ZA  ZB

比对去掉umi的序列(hg38参考基因组，有需要可以换成hg19)

samtools fastq  test.umi.uBam | bwa mem -t 50 -p /data/ref/hg38/hg38 /dev/stdin | samtools view -b > test.umi.Bam

合并uBam 和 Bam 得到带有UMI信息的比对文件

picard MergeBamAlignment R=/data/ref/hg38/hg38.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

Call Consensus Reads

fgbio GroupReadsByUmi \
    --input=test.umi.merged.Bam \
    --output=test.umi.group.Bam  \
    --strategy=paired  --min-map-q=20  --edits=1 --raw-tag=RX
 
fgbio CallMolecularConsensusReads \
    --min-reads=1 \
    --min-input-base-quality=20 \
    --input=test.umi.group.Bam \
    --output=test.consensus.uBam
 
samtools fastq test.consensus.uBam | bwa mem -t 50 -p /data/ref/hg38/hg38  /dev/stdin | samtools view -b - > test.consensus.Bam
 
picard MergeBamAlignment R=/data/ref/hg38/hg38.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
 
fgbio FilterConsensusReads \
    --input=test.consensus.merge.Bam  \
    --output=test.consensus.merge.filter.Bam \
    --ref=/data/ref/hg38/hg38.fa --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
 
fgbio ClipBam \
    --input=test.consensus.merge.filter.Bam   \
    --output=test.consensus.merge.filter.clip.Bam \
    --ref=/data/ref/hg38/hg38.fa  \
    --clip-overlapping-reads=true

拿到bam之后，就需要进行正常的call 突变流程了。

大家流程可以参考外显子（WES）panel分析基础篇

当然了，我们不能人工一步一步的跑这个过程，需要直接封装好直接把带UMI的数据分析出我们想要的结果，包括ANNOVAR注释，OncoKB注释，MSI、CNV、SV等，封装脚本的示例如下:

bash  gatk_UMI_sh sample.clean.R1.fastq.gz sample.clean.R2.fastq.gz sample(样本名称) 输出路径 -t bed文件  -n normal_fq1 -N normal_fq2  -c CNV基线  -y 癌症类型

完整的代码如下：

#!/bin/bash
####zhaoyunfei

###############################################################
###############################################################

###############################################################
#path of software
fastp=fastp软件路径
bwa=/home/jinsai/software/bwa/bwa-master/bwa软件路径
samtools=/home/jinsai/software/samtools/samtools-1.9/samtools软件路径
gatk1=gatk软件路径
annovar=table_annovar.pl软件路径
cnvkit=cnvkit.py软件路径
factera=factera.pl软件路径
bedtools=bedtools软件路径
sambamba=sambamba软件路径
msi=msisensor2软件路径
msi_model=models_hg19_GRCh37
msisensor_pro=msisensor_pro软件路径
picard=picard路径
fgbio=fgbio路径
###############################################################

#############################################################################################
#reference file
hg19=ucsc.hg19.fasta路径
Mills_and_1000G_gold_standard=Mills_and_1000G_gold_standard.indels.hg19.sites.clean.vcf路径
DBSNP=dbsnp_138.hg19.clean.vcf路径
INDELS=1000G_phase1.snps.high_confidence.hg19.sites.clean.vcf路径
MILLS=Mills_and_1000G_gold_standard.indels.hg19.sites.clean.vcf路径
cnvkit_reference=reference.cnn路径
af_only_gnomad=somatic-b37_af-only-gnomad.new.sites.vcf路径
Mutect2_Filter=Mutect2_Filter路径
humandb=humandb路径
ref_bit=ucsc.hg19.2bit路径
exons_bed=exons.bed路径
refFlat=refFlat.txt路径
access_hg19=access.hg19.bed路径
small_exac_common=somatic-b37_small_exac_common_3.new.vcf路径
#############################################################################################

absname(){
    basepath=$(cd `dirname $1`; pwd)
    basefile=$(basename $1)
    echo $basepath"/"$basefile
}
function usage(){
    echo -e "Usage: `basename $0` "'fq1 fq2 library out_dir'
    echo -e "    -n <normal fq1>"
    echo -e "    -N <normal fq2>"
    echo -e "    -m|--max_threads <Integer: Default value: 24. max threads >"
    echo -e "    -t|--target <target file, default:Panel.bed>"
    echo -e "    -c|--reference_cnn < reference_cnn file ,required = T > "
    echo -e "    -y|--cancer_type < cancer type > "
}

nor_fq1=''
nor_fq2=''
threads=36
ARGS=`getopt -o n:N:t:m:c:y:h --long a:,N:,target,max_threads,reference_cnn,cancer_type,help -n "$0" -- "$@"`
if [ $? != 0 ]; then
    usage;
    exit 1
fi
eval set -- "${ARGS}"

while true
do
    case "$1" in
        -h|--help)
            usage;
            exit;
            ;;
        -n|--n)
            case "$2" in
                "")
                    shift 2
                    ;;
                *)
                    n_fq1=$2
                    shift 2;
                    ;;
            esac
            ;;
        -N|--N)
            case "$2" in
                "")
                    shift 2
                    ;;
                *)
                    n_fq2=$2
                    shift 2;
                    ;;
            esac
            ;;
        -t|--target)
            case "$2" in
                "")
                    shift 2
                    ;;
                *)
                    target=$2
                    shift 2;
                    ;;
            esac
            ;;
        -m|--max_threads)
            case "$2" in
                "")
                    shift 2
                    ;;
                *)
                    threads=$2
                    shift 2;
                    ;;
            esac
            ;;
        -c|--reference_cn