Skip to content

排序和标记重复

排序和标记重复都是为了后面更好的找变异，从gatk best practice来说，还需要一部加入测序信息的步骤。排序和标记重复均可使用samtools或者picard进行。但是，gatk4已经内置了picard的全部功能，这里将会使用gatk4进行这些操作。如果是不是走gatk call snp的流程，其实更建议使用samtools。这一步虽然名字叫排序和标记重复，实际做的事情不止这些，只是排序和标记重复是最重要的。

转为二进制bam格式

在进行下面的操作之前，先把sam文件转为二进制的bam文件，可以节省空间以及提高效率。 使用samtools进行转化。 

1
samtools view -bSh B17NC.sam > B17NC.bam

建议将此步更新为使用sambamba，速度更快。 

1
sambamba view -S -h -f bam B17NC.sam -t 8 -o B17NC.bam

排序

使用gatk4进行排序。按coordinate进行排序，意思就是按照染色体编号进行排序。 

1
gatk SortSam -I B17NC.bam -O B17NC.sorted.bam -SO coordinate

嫌弃gatk4太慢的，也可以用samtools。事实上，如果对bam文件的要求不怎么严谨，那么这一步已经可以算是最后一步了，当前生产出的bam已经可以用于call snp了。 

1
samtools sort B17NC.bam -@ 8 -o B17NC.sorted.bam

建议此步使用sambamba，速度更快。 

1
sambamba sort B17NC.bam -t 8 -o B17NC.sorted.bam

标记重复序列

标记重复序列实际上修改bam文件的flag值，重复reads的flag会加上1024用于标记。标记重复是为了去除PCR时产生的大量重复，获得较准确的突变丰度。另外，部分标记重复软件会形成新的tag用于标记，可使用picard/gatk等对tag来进行去重。

这里使用gatk4进行。

1
gatk MarkDuplicates -I B17NC.sorted.bam -O B17NC.mdup.bam -M B17NC.dups.txt

此步可以使用sambamba，速度更快，回报格式与picard/gatk等同。 

1
sambamba markdup -t 8 B17NC.sorted.bam -o B17NC.mdup.bam

UMI

对游离DNA进行超高深度测序时一般会加入UMI序列，去重步骤与不加入UMI略有不同。可使用fastp 加上gencore的流程进行去重。但是gencore的去重方式是直接去掉而不是标记，并且也不是改写flags值（加上1024）的方式，因此这里使用fgbio 加上 gatk （gatk已包含picard的功能）的方式。

IDT有一篇主要使用fgbio进行去重的流程，可以参考。而以下的流程仅使用fgbio进行UMI提取，再使用gatk MarkDuplicates进行去重。

先将过滤后的fastq转为ubam。 

1
2
3
gatk FastqToSam -F1 clean.1.fq.gz -F2 clean.2.fq.gz \
    -O sample.ubam \
    -SM sample -PL illumina -PU unit

然后提取UMI，以双端3bp在Index的3‘的UMI 为例，用的结构式是3M1S+T，即3bp UMI然后跳1bp，剩余的都是Template。更多的Reads结构写法指导参考此篇。

1
2
3
java jar fgbio.jar ExtractUmisFromBam \
    -i sample.ubam -o sample.umi.ubam \
    -r 3M1S+T 3M1S+T -t ZA ZB -s RX

然后将ubam比对到参考基因组，要注意的是，由于fgbio默认将Read Group中的ID写为A，所以此时比对最好也是将ID写为A，后续可以用gatk AddOrReplaceReadGroups来修改。主要bwa mem的-p参数必不可少。

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

将ubam和bam合并，使比对后的bam中能带上未比对的ubam中的UMI tag。

1
2
3
4
5
6
gatk MergeBamAlignment \
    -R hg38.fa \
    -ALIGNED sample.umi.bam -UNMAPPED sample.umi.ubam \
    -O sample.umi.merge.bam \
    --ALIGNER_PROPER_PAIR_FLAGS true \
    --ATTRIBUTES_TO_RETAIN XS

最后还是使用MarkDuplicates来去重，记得--BARCODE_TAG RX这个参数。

1
2
3
4
5
6
gatk MarkDuplicates \
    -I sample.umi.merge.bam \
    -O sample.marked.bam \
    -M sample.dups.txt \
    --CREATE_INDEX true \
    --BARCODE_TAG RX

添加或改写头信息

使用gatk4来改变头信息，如果后续用gatk4来进行体细胞突变分析，头信息还是很重要的。当然，如果使用bwa比对时已经用来-R来加入信息，就不需要这一步了。 

1
2
gatk AddOrReplaceReadGroups -I B17NC.mdup.bam -O B17NC.addhead.bam \
    -LB lib1 -PL illumina -PU unit1 -SM B17NC

创建索引

可以使用gatk4创建索引，也可以使用samtools。 

1
gatk BuildBamIndex -I B17NC.addhead.bam

更建议使用samtools 

1
samtools index B17NC.addhead.bam

质量值校正

接下来就是使用gatk的质量值校正部分了。gatk4使用机器学习的方式，以各种现有的数据集作为训练集来校正质量值。实际上，在bwa比对后再经过samtools排序的bam文件就已经可以进行call snp了。只不过gatk多一手校正。

indel区域重比对

1
2
3
4
5
6
gatk BaseRecalibrator -I B17NC.addhead.bam \
    -R Homo_sapiens_assembly38.fasta \
    --known-sites Mills_and_1000G_gold_standard.indels.b38.vcf.gz \
    --known-sites dbsnp_138.vcf \
    --known-sites 1000g_phase.indels.vcf \
    -O B17NC.recal.table

重比对区域校正

1
2
3
4
gatk ApplyBQSR -R Homo_sapiens_assembly38.fasta \
    -I B17NC.addhead.bam \
    --bqsr-recal-file B17NC.recal.table \
    -O B17NC.final.bam