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

大家好，欢迎来到IT知识分享网。

作者，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_cnn) case "$2" in "&#