生物信息文件夹

再整理一次测序数据去重流程

  • 2021年 08月22日

二代测序PCR过程中会产生duplications，为了下游分析的正确，一般需要进行去重操作。最常用的去重工具是picard MarkDuplicates。picard MarkDuplicates默认计算比对后的Reads，当存在Start与End以及序列一致的情况时，再计算这些reads的比对质量值之和，取其中最大的作为模板，其他作为duplications并在flag值中加上1024进行标记。

参考xGen Prism DNA Library Prep Kit。

以下包括常规的MarkDuplicates去重流程、有UMI下的MarkDuplicates去重流程，以及单端和双端的fgbio去重流程。

无UMI

使用组织作为样本检测时，很少会加入UMI序列，在比对后，使用MarkDuplicates进行去重。 gatk4已集成picard所有功能，所以使用gatk4的MarkDuplicates进行去重。默认是仅标记重复，不去除重复。

去重

gatk MarkDuplicates \
	-I sample.bam -O sample.marked.bam -M sample.dups.txt

也可以使用速度更快的sambamba，去重策略与MarkDuplicates一致。

sambamba markdup -t 8 -p sample.bam sample.marked.bam

有UMI

首先要提取UMI序列，然后再各种方式去重。

提取UMI

先把fastq转为ubam。

gatk FastqToSam \
	-F1 sample_1.fq.gz \
	-F2 sample_2.fq.gz \
	-O sample.ubam \
	-SM sample \
	-PL illumina \
	-PU Unit \
	-RG MarkDups

再使用fgbio提取UMI序列。如4bp UMI，一般会跳过UMI后的2bp碱基。

java -jar fgbio.jar ExtractUmisFromBam \
	-i sample.ubam \
	-o sample.umi.ubam \
	-r 4M2S+T 4M2S+T -t ZA ZB -s RX

进行比对。

samtools fastq sample.umi.ubam \
	| bwa mem -t 8 -p -R "@RG\tPL:illumina\tSM:sample\tPU:unit\tID:MarkDups" \
	ref.fa /dev/stdin \
	| samtools view -bSh - > sample.umi.bam

把提取UMI后未比对的ubam文件与比对后的bam文件合并，即把ubam中的UMI tag加入到bam中。

gatk MergeBamAlignment \
	-R ref.fa \
	-ALIGNED sample.umi.bam \
	-UNMAPPED sample.umi.ubam \
	-O sample.bam \
	--ALIGNER_PROPER_PAIR_FLAGS true \
	--ATTRIBUTES_TO_RETAIN XS \
	-MAX_GAPS -1 -ORIENTATIONS FR \
	--VALIDATION_STRINGENCY SILENT

去重

使用MarkDuplicates进行去重。记得加上–BARCODE_TAG RX。进行此步已经能去重，但是由于MarkDuplicates原理是保留最高MQ值的reads（默认参数），并不能很好的将测序错误合并为一致性reads。对于测序错误较多的数据，可选用下面的其他流程而不是使用MarkDuplicates。

gatk MarkDuplicates \
	-I sample.bam \
	-O sample.marked.bam \
	-M sample.dups.txt \
	--BARCODE_TAG RX

有UMI 续

UMI校对

同样进行了前面的提取UMI序列的步骤后，如果有确切的UMI序列信息，可以进行这一步。对于长UMI，mismatches（-m）可以设大一点，而短UMI则设小一点。把已知的UMI序列写入Umis.txt，一行一个。

java -jar fgbio.jar CorrectUmis \
	-i sample.bam \
	-o sample.fixUMI.bam \
	-m 1 -d 1 \
	-M sample.metrics.txt \
	-r reject.bam \
	-U Umis.txt -t RX

分组

对提取UMI后的sample.bam或校正后的sample.fixUMI.bam进行分组。其中-m是质量值，-e是允许修改的UMI之间允许修改的值。

对于单端UMI数据

java -jar fgbio.jar GroupReadsByUmi \
	-s Adjacency -i sample.bam -o sample.group.bam \
	-m 15 -e 0

对于双端UMI数据

java -jar fgbio.jar GroupReadsByUmi \
	-s Paired -i sample.bam -o sample.group.bam \
	-m 15 -e 0

获得一致性序列

这里把获得一致性序列中的-M设置为1，在后续的filter中还可以进行调整，比较灵活，但是在这步会牺牲计算速度。

对于单端UMI数据

java -jar fgbio.jar CallMolecularConsensusReads \
	-i sample.group.bam -o sample.con.ubam \
	-1 45 -2 30 -R MarkDups -M 1 -m 15

对于双端UMI数据

java -jar fgbio.jar CallDuplexConsensusReads \
	-i sample.group.bam -o sample.con.ubam \
	-1 45 -2 30 -R MarkDups -M 1 -m 15

再次比对

再次比对并合并，将信息加入到比对后数据中，merge步骤会默认将reads顺序调整为croodinate，所以不用再排序。

samtools fastq sample.con.ubam \
	| bwa mem -t 8 -p -R "@RG\tPL:illumina\tSM:sample\tPU:unit\tID:MarkDups" \
	ref.fa /dev/stdin \
	| gatk MergeBamAlignment \
	-R ref.fa \
	-ALIGNED /dev/stdin \
	-UNMAPPED sample.con.ubam \
	-O sample.con.merge.bam \
	--ALIGNER_PROPER_PAIR_FLAGS true \
	--ATTRIBUTES_TO_RETAIN XS \
	-MAX_GAPS -1 -ORIENTATIONS FR \
	--VALIDATION_STRINGENCY SILENT

过滤

对合并后的数据进行过滤

java -jar fgbio.jar FilterConsensusReads \
	-i sample.con.merge.bam \
	-o sample.con.merge.filter.bam \
	-r ref.fa -R true -M 3 -E 0.05 -N 4 -e 0.1 -n 0.1

Clip

将重合区域进行clip，避免多次计算

java -jar fgbio.jar ClipBam \
	-i sample.con.merge.filter.bam -o sample.con.merge.filter.clip.bam \
	-r ref.fa --clip-overlapping-reads true
搜索