Mutect2-肿瘤不同转移时期的免疫微环境异质性研究

本来只是单纯想用Mutect2来重复下文献的分析过程，结果变成了Mutect2的使用笔记。。GATK的Mutect2相关文档非常详细，思路及理念也讲的很清楚，FAQ也很完善 https://gatkforums.broadinstitute.org/gatk/discussions/tagged/mutect2 ，是个很好的学习过程

Mutect2 resources guide

文献： Heterogeneous Tumor-Immune Microenvironments among Differentially Growing Metastases in an Ovarian Cancer Patient

别人的解读： 文章解读 | 肿瘤不同转移时期的免疫微环境异质性研究

文献的相关数据：SRA的study id是SRP095988，进入NCBI的网站： https://trace.ncbi.nlm.nih.gov/Traces/study/?acc=SRP095988 ，可知其样本总共有6个，分别是VagCuff-4-1、VagCuff-4-1、RUQ-1-4、LiverSeg5-7-1、 Normal_blood和Spleen-9-5

进入NCBI的ftp： ftp://ftp.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByStudy/sra/SRP/SRP095/SRP095988/ 下载每个样本对应的sra文件，由于网速不太快，因此首选aspera来下载数据，命令如下：

ascp -T -i /home/anlan/.aspera/connect/etc/asperaweb_id_dsa.openssh anonftp@ftp-private.ncbi.nlm.nih.gov:sra/sra-instant/reads/ByStudy/sra/SRP/SRP095/SRP095988/SRR5141040/SRR5141040.sra ./

接着用NCBI的sratoolkit工具将sra文件转化为fastq文件：

fastq-dump --split-3 ./SRR5141040.sra

然后将SRR编号全都改为样本名字（个人喜欢哈），再简单做个质控，使用较为好用fastp质控工具（使用默认参数，可以看下默认参数都是哪些过滤指标）：

fastp -i Normal_blood_1.fastq -I Normal_blood_2.fastq -o Normal_blood_1.clean.fq -O Normal_blood_2.clean.fq

从文献中可看出，作者是用Broad Institute’s Agilent exome design (Agilent SureSelect All Exon V2)，理论上后续分析最好也是用其一致的bed文件，但我这暂时用CCDS代替下了

测序长度是76 bp paired，测序深度是：

The sequencing depths of the samples were: normal blood sample (90% at 20X), primary (82% at 50X), spleen (78% at 50X), RUQ (60% at 50X), liver (89% at 50X), and vaginal cuff (77% at 50X) tumors.

文章中未提使用了哪个比对软件，所以我就按照GATK 4.0流程，从Raw reads -> Map to

-> Raw Mapped Reads -> Mark Duplicated -> Recalibrate Base Quality Scores -> Analysis-Rready Bam

参照之前的文章（用GATK 4.0处理WES数据来找variant） GATK 4.0 全外显子call variant

shell脚本如下（如Normal_blood样本，则： nohup bash bwa-BQSR.sh Normal_blood_1.clean.fq Normal_blood_2.clean.fq Normal_blood & ）：

#!/bin/bash

fq1=$1
fq2=$2
sample=$3

index=/home/anlan/reference/index/bwa/gatk_hg38/gatk_hg38
genome=/home/anlan/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta
GATK=/home/anlan/biosoft/GATK4.0/gatk-4.0.5.1/gatk
exon_interval=/home/anlan/exon.intervals 
dbsnp=/home/anlan/annotation/GATK/resources/bundle/hg38/dbsnp_146.hg38.vcf.gz
dbsnp1000G=/home/anlan/annotation/GATK/resources/bundle/hg38/1000G_phase1.snps.high_confidence.hg38.vcf.gz
dbindel100G=/home/anlan/annotation/GATK/resources/bundle/hg38/Mills_and_1000G_gold_standard.indels.hg38.vcf.gz

## BWA Alignment
if [ ! -e $sample.sam ]; then
    bwa mem -t 4 -M -R "@RG/tID:$sample/tPL:illumina/tLB:WES/tSM:$sample" $index $fq1 $fq2 1>$sample.sam 2>/dev/null
    echo bwa `date`
fi

## Sam to bam and sort
if [ ! -e $sample.sorted.bam ]; then
    $GATK --java-options "-Xmx4G -Djava.io.tmpdir=/home/anlan/tmp" SortSam /
    -I $sample.sam -O $sample.sorted.bam /
    -SO coordinate --CREATE_INDEX true
    echo gatk-SortSam `date`
fi

## MarkDuplicate
if [ ! -e $sample.sorted.marked.bam ]; then
    $GATK --java-options "-Xmx4G -Djava.io.tmpdir=/home/anlan/tmp" MarkDuplicates /
    -I $sample.sorted.bam -O $sample.sorted.marked.bam /
    -M ${sample}.metrics
    echo gatk-MarkDuplicates `date`
fi

## BQSR
if [ ! -e $sample.sorted.marked.BQSR.bam ]; then
    $GATK --java-options "-Xmx4G -Djava.io.tmpdir=/home/anlan/tmp" BaseRecalibrator /
    -R $genome -I $sample.sorted.marked.bam -O ${sample}.recal_data.table /
    --known-sites $dbsnp --known-sites $dbsnp1000G --known-sites $dbindel100G

    $GATK --java-options "-Xmx4G -Djava.io.tmpdir=/home/anlan/tmp" ApplyBQSR /
    -R $genome -I $sample.sorted.marked.bam --bqsr-recal-file ${sample}.recal_data.table /
    -O $sample.sorted.marked.BQSR.bam
    echo gatk-BQSR `date`
fi

对排序后的bam文件进行测序深度分析，看看样本的测序深度是多少，使用bedtools的genomecov工具和intersect工具，外加脚本来计算，以Normal_blood为例，如下（其中exon.bed是根据CCDS文件制造的exon的bed文件）：

bedtools genomecov -ibam Normal_blood.sorted.bam -bg >Normal_blood.bedGraph
bedtools intersect -a Normal_blood.bedGraph -b ../exon.bed >Normal_blood.intersect
python3.4 bedtools2depth.py

还试了samtoools的depth工具：

samtools depth -a -b ../exon.bed Normal_blood.sorted.bam >Normal_blood.samtoool.depth

还用了一个Github开发者写的计算测序深度/覆盖度工具 mosdepth （速度比samtools快很多），我下载了二进制版本，如果软件使用时报错：无法找到htslib.so文件，则你需要在脚本前面再加上htslib.so文件所在路径（第一次见到是这么玩的。。。）

LD_LIBRARY_PATH=~/biosoft/HTSlib/htslib-1.5 ./mosdepth --by ../exon.bed Normal_blood.depth Normal_blood.sorted.bam

但用上述统计的测序深度均和文章中的不一致。。k可能我是只统计外显子上的平均深度，未考虑两端测序长度（主要也不知道作者是拿质控后的bam文件还是原始比对bam文件做的测序深度以及捕获区域是怎么设计的）

文中作者对bam文件还做了一步过滤：

Reads with mapping quality below 30 in the BAM files were filtered out before mutation calling.

因此我也需要对每个样本 XXX.sorted.marked.BQSR.bam ，做下过滤（不确定应该是否在BQSR之前过滤还是之后），如下:

samtools view -b -q 30 Normal_blood.sorted.marked.BQSR.bam >Normal_blood.allready.bam
samtools index Normal_blood.allready.bam

接下来则是mutect2来call somatic variant，按照GATK团队建议的，最好是将正常样本建一个建一个Panel of Normals (PoN)，如 Panel of Normals (PON) 所说：

What all PONs have in common is that (1) they are made from normal samples (in this context, “normal” means derived from healthy tissue that is believed to not have any somatic alterations) and (2) their main purpose is to capture recurrent technical artifacts in order to improve the results of the variant calling analysis.

按照官方的说明，normals samples是越多越好（至少40个），但是其实一些小样本（4-10个）也是可以的，总比没有好；就算是一些unmatched PON（不是同一批样本，但是在同一测序仪上进行测序的？），也可以作为PON来过滤掉一些测序仪以及实验人为误差 Create a sites-only PoN with CreateSomaticPanelOfNormals ；如果没有组织样本，但是如果是有多个人的Normal blood对照样本，也可以用来构建PON Using Germline Blood for PON creation

但是按照这篇文献的数据来看，是无法构建PON的，其样本是来自于一个人的不同部位的组织样本以及血液样本，因此我就直接拿肿瘤样本和正常血液样本用mutect2来call somatic variant

这里有一点需要注意下，在Mutect2之前（Mutect版本），是在call variant时将dbSNP作为blacklist，并且将Cosmic作为whitelist，但是在Mutect2中取消了这种做法 GATK4 beta Mutect2 misses –cosmic option ；因为GATK团队选择用PON以及germline resource来帮助过滤possible germline variants

在官网文档的文档中是使用含有population allele frequencies的gnomAD数据集，我也看到过用一些其它germline resource，如1000 Genomes

、hapmap等；前提是这些数据集的INFO中必需要有AF

gnomAD数据集（af-only-gnomad.hg38.vcf.gz）在GATK的官方ftp的resource里有，在路径里 ftp://ftp.broadinstitute.org/bundle/ 找一找

因此对tumor和matched normal进行calls somatic variants命令如下（直接使用了默认参数，有其他需求时再调整，-A可对生成的vcf文件增加注释参数）：

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk --java-options "-Xmx2g" Mutect2 /
-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
-I Primary-IIIG.allready.bam /
-tumor Primary-IIIG /
-I Normal_blood.allready.bam /
-normal Normal_blood /
--germline-resource ~/annotation/GATK/resources/bundle/hg38/af-only-gnomad.hg38.vcf.gz /
--af-of-alleles-not-in-resource 0.0000025 / 
--disable-read-filter MateOnSameContigOrNoMappedMateReadFilter /
-L ../exon_150bp.list /
--bam-output Primary-IIIG2Normal.bam /
-O Primary-IIIG.mutect2.vcf

上面需要注意的是 -tumor 和 -normal 参数所对应的样本名称必需跟bam文件中的SM一致，也就是在bwa比对时所加的SM； --af-of-alleles-not-in-resource 对应的0.0000025是根据你所加的germline resource来定的，比如gnomAD resource有~200k exomes and ~16k genomes，因此值对应就是1/(2*exome samples)； --disable-read-filter MateOnSameContigOrNoMappedMateReadFilter 主要是将那些比对到不同contig上的paired reads也参与分析；而加 -L 参数则为了减少运行时间，对于目标区域做call somatic variant；最后加的 --bam-output 则是为了生成重新组装的比对bam文件（毕竟Mutect2的算法是基于贝叶斯分类和单倍体重组装），从而用于最后的人工检查（如IGV） Call somatic short variants and generate a bamout with Mutect2

如果你考虑cross-sample contamination对somatic突变假阳性的影响，那么可以选择 GetPileupSummaries 工具来计算肿瘤样本的污染情况；对于现在的Mutect2的流程，GATK只考虑了normal contamination of tumor这种情况，至于tumor contamination of normal是不考虑的；对于评估cross-sample contamination，其基本思路是来自GATK之前开发的ContEst工具，其中需要的一个主要文件是population germline resource，并且这个resource文件只包含common biallelic variants，什么是biallelic variants可以看 Biallelic vs Multiallelic sites

以gnomAD resource为例，先用SelectVariants工具处理下，生成一个 af-only-gnomad.hg38.SNP_biallelic.vcf.gz 用于后续分析

/biosoft/GATK4.0/gatk-4.0.5.1/gatk SelectVariants /
-R ~/reference/genome/gatk_hg38/Homo_sapiens_assembly38.fasta /
-V af-only-gnomad.hg38.vcf.gz /
--select-type-to-include SNP /
--restrict-alleles-to BIALLELIC /
-O af-only-gnomad.hg38.SNP_biallelic.vcf.gz

接着用GetPileupSummaries计算resource site位点上的count数，这样并不是计算所有的位点，而是用AF参数过滤后的（大于0.01并小于0.2），同时也可以用 -L 参数指定区域，分别对Normal和Tumor样本计算

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk GetPileupSummaries /
-I Normal_blood.allready.bam /
-L ../exon_150bp.list /
-V ~/annotation/GATK/resources/bundle/hg38/af-only-gnomad.hg38.SNP_biallelic.vcf.gz /
-O Normal_blood.pileups.table

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk GetPileupSummaries /
-I Primary-IIIG.allready.bam /
-L ../exon_150bp.list /
-V ~/annotation/GATK/resources/bundle/hg38/af-only-gnomad.hg38.SNP_biallelic.vcf.gz /
-O Primary-IIIG.pileups.table

然后结合上述Normal和Tumor数据评估污染比例，基于的原理简单的说是：

The basic idea, which comes from ContEst5 by Kristian Cibulskis and others in the Broad Institute Cancer Genome Analysis group, is simply to count ref reads at hom alt sites and subtract the number of ref reads expected from sequencing error to obtain the number of ref reads contaminating these hom alt sites.  Finally, we use the allele frequencies to account for the fact that some contaminating reads have the alt allele. The only subtlety is in distinguishing hom alt sites from loss of heterozygosity events, which we describe below.

具体原理可以看了 mutect.pdf 中的CALCULATING CONTAMINATION部分和文献 ContEst: Estimating cross-contamination of human samples in next-generation sequencing data ，我看了后也不太懂。。。总之最后通过计算出的污染比例（XXX.calculatecontamination.table文件中）来过滤掉somatic variant中的一些可能是污染导致的假阳性突变

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk CalculateContamination /
-I Primary-IIIG.pileups.table /
-matched Normal_blood.pileups.table /
-O Primary-IIIG.calculatecontamination.table

2016年Bioinformatics有一篇文章也介绍一个评估污染的工具 Conpair

最后一步则是过滤了，结合上述的cotamination评估（比如我Primary-IIIG样本的污染评估是0.00131826，那么在过滤后的VCF里，AF低于0.00131826的位点则会被标注上contamination）,下述命令使用默认过滤参数，具体情况可以再修改

~/biosoft/GATK4.0/gatk-4.0.5.1/gatk FilterMutectCalls /
-V Primary-IIIG.mutect2.vcf /
--contamination-table Primary-IIIG.calculatecontamination.table /
-O Primary-IIIG.somatic.vcf

最终vcf的FILTER列中有很多过滤信息

grep '##FILTER' Primary-IIIG.somatic.vcf

获得过滤后的somatic vcf文件后，然后用突变注释工具（ANNOVAR、snpEff/VEP等工具）注释下，就可以了，最后再人工校验下

对于每个过滤指标的解释，可以看 mutect.pdf 中的MUTECT FILTERS部分

再把文章好好看下，把自己分析的结果跟文献对照下，虽然我大部分重复文献最终的结局都是对不上的。。。可能打算再结合下其他somatic variant软件的结果，比如Strelka？

本文出自于 http://www.bioinfo-scrounger.com 转载请注明出处