TCGA计划的4个找somatic mutation的软件使用体验

体细胞突变（somatic mutation）是指患者某些组织或者器官后天性地发生了体细胞变异，虽然它不会遗传给后代个体，却可以通过细胞分裂，遗传给子代细胞。体细胞突变对肿瘤的发生发展有关键性的作用，并且它也是制定肿瘤癌症靶向治疗措施的关键所在。NGS使体细胞变异的检测更加全面，成本更低，在检测多种体细胞变异上具有很大的优势，但在使用过程中还存在着挑战：如样品降解、覆盖度不足、遗传异质性和组织污染（杂质)等问题。 为应对以上挑战，降低错误率，科学家采取了不同的算法和统计模型用于检测体细胞突变。目前最受欢迎的有Varscan、SomaticSniper、 Strelka 和MuTect2 。

这些软件大都是直接对肿瘤-正常样本的每个位点进行比较，对肿瘤样本中明显高于正常样本的次等位基因进行标记，作为体细胞变异，同时排除种系突变和杂合性丢失（LOH）情况。虽然这些软件具有较高的引用率，并在不断地更新，但仍存在不足：

a 缺乏完整可靠的实验来评估检测结果；

b 缺乏金标准，不能保证检测到的灵敏度和特异性最高；

c 在实际应用中，各软件的相对优缺点在很大程度上是未知的。

下面是TCGA计划采取的软件：

• MuSE
• varscan

• MuTect
• SomaticSniper

当然，能找somatic mutation的软件还有很多，比如Strelka等，就不一一讲解啦。其实最基础的原理都是应该是 除去在normal样本里面出现过的germline变异位点，可以很简单的GATK UnifiedGenotyper followed by simple subtraction即可。

首先用mutect2

这个软件已经被整合到GATK里面啦，所以下载GATK即可使用它。

java软件，下载即可使用，GATK软件下载以前需要自行注册，目前是开放下载了，使用代码是：

GENOME=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta 
GATK=/home/jianmingzeng/biosoft/GATK/GenomeAnalysisTK.jar
DBSNP=/home/jianmingzeng/annotation/variation/human/dbSNP/All_20160601.vcf.gz
normal_bam=NPC10F-N_recal.bam
tumor_bam=NPC10F-T_recal.bam
sample=NPC10F
​
​
java -Djava.io.tmpdir=$TMPDIR   -Xmx40g -jar $GATK  -T MuTect2 \
-R $GENOME -I:tumor $tumor_bam-I:normal $normal_bam \
--dbsnp  $DBSNP   -o ${sample}-mutect2.vcf

出来的结果里面有些比较陌生的tags，需要仔细理解，这样才能看懂vcf文件，并进行进一步的过滤。

##FILTER=<ID=alt_allele_in_normal,Description="Evidence seen in the normal sample">
##FILTER=<ID=clustered_events,Description="Clustered events observed in the tumor">
##FILTER=<ID=germline_risk,Description="Evidence indicates this site is germline, not somatic">
##FILTER=<ID=homologous_mapping_event,Description="More than three events were observed in the tumor">
##FILTER=<ID=multi_event_alt_allele_in_normal,Description="Multiple events observed in tumor and normal">
##FILTER=<ID=panel_of_normals,Description="Seen in at least 2 samples in the panel of normals">
##FILTER=<ID=str_contraction,Description="Site filtered due to contraction of short tandem repeat region">
##FILTER=<ID=t_lod_fstar,Description="Tumor does not meet likelihood threshold">
##FILTER=<ID=triallelic_site,Description="Site filtered because more than two alt alleles pass tumor LOD">

旧版的mutect只是对一个位点REJECT或者PASS，但是新版增加了多种情况来解释为什么REJECT，就是上面的那些tag的组合。

即使我把N-T反过来用mutect2来call somatic mutation，仍然会有125个位点PASS，只需要是在tumor里面纯粹的野生型，在normal里面是AF非常低的杂合即可。

其次是varscan

java软件都是下载即可使用，官网可以下载，我把它下载到了~/biosoft/VarScan/VarScan.v2.3.9.jar目录，使用代码是：

​
​
TMPDIR=/home/jianmingzeng/tmp/software
normal_bam=NPC10F-N_recal.bam
tumor_bam=NPC10F-T_recal.bam
sample=NPC10F
​
normal_pileup="samtools mpileup -q 1 -f $reference $normal_bam";
tumor_pileup="samtools mpileup -q 1 -f $reference $tumor_bam";
​
# Next, issue a system call that pipes input from these commands into VarScan :
java -Djava.io.tmpdir=$TMPDIR   -Xmx40g  -jar ~/biosoft/VarScan/VarScan.v2.3.9.jar \
somatic <($normal_pileup) <($tumor_pileup) $sample
java -jar ~/biosoft/VarScan/VarScan.v2.3.9.jar processSomatic $sample.snp

结果文件如下：

  101353 NPC10F.snp
   95145 NPC10F.snp.Germline
    4831 NPC10F.snp.Germline.hc
    3503 NPC10F.snp.LOH
    3485 NPC10F.snp.LOH.hc
    2343 NPC10F.snp.Somatic
     557 NPC10F.snp.Somatic.hc

其中值得关注的就是NPC10F.snp.Somatic.hc文件里面的位点，打开浏览可以看出其实就是normal样本里面的野生型变化到tumor里面的杂合型，就是allel frequency的增加。

接着是muse

首先需要下载安装

cd ~/biosoft
mkdir muse &&  cd muse
wget http://bioinformatics.mdanderson.org/Software/MuSE/MuSEv1.0rc_submission_b391201 muse
mv MuSEv1.0rc_submission_b391201 muse
chmod 777 muse

该软件需要的input跟mutect2一样，都是：

• (1) the indexed reference genome FASTA file,
• (2) the binary sequence alignment/map formatted (BAM) sequence data from the pair of tumor and normal DNA samples, and
• (3) the dbSNP variant call format (VCF) file that should be bgzip compressed, tabix indexed and based on the same reference genome as (1).

它给的测试代码是：

./MuSE call –O Output.Prefix –f Reference.Genome Tumor.bam Matched.Normal.bam
./MuSE sump -I Output.Prefix.MuSE.txt -G –O Output.Prefix.vcf –D dbsnp.vcf.gz

需要分成两步来运行，修改一下就可以处理自己的样本啦。

reference=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta
DBSNP=/home/jianmingzeng/annotation/variation/human/dbSNP/All_20160601.vcf.gz
normal_bam=NPC10F-N_recal.bam
tumor_bam=NPC10F-T_recal.bam
sample=NPC10F
~/biosoft/muse/muse call -O $sample -f $reference $tumor_bam $normal_bam
# ~/biosoft/muse/muse sump -I ${sample}.MuSE.txt -E –O ${sample}_muse.vcf –D $DBSNP

看看结果文件，就是典型的VCF格式，而且tags不多值得注意的只有ID=SS,Number=1,Type=Integer,Description="Variant status relative to non-adjacent Normal,0=wildtype,1=germline,2=somatic,3=LOH,4=post-transcriptional modification,5=unknown

前5列CHROM POS ID REF ALT很正常，第6列QUAL全部是点，第7列FILTER 把位点分级了，统计如下：

     81 Tier1
     26 Tier2
     12 Tier3
     26 Tier4
     39 Tier5

第8列是 INFO 信息，全部是SOMATIC

第9,10,11列是GT:DP:AD:BQ:SS 格式的tumor和normal，可以看到normal都是野生型0/0, tumor全部是杂合突变1/0，只是allel frequency不同而已，介于0~1之间。

最后是SomaticSniper

第一次用这个软件，也是需要先安装

cd ~/biosoft
mkdir SomaticSniper &&  cd SomaticSniper
wget https://github.com/genome/somatic-sniper/archive/v1.0.5.0.tar.gz
tar zxvf v1.0.5.0.tar.gz
cd somatic-sniper-1.0.5.0
mkdir build && cd build 
cmake ../
make deps
make -j
make test

服务器的make好像有点问题，没有安装成功，先略过吧。

还看看strelka吧

这款软件专门设计给临床使用的，说明书非常详细 https://github.com/Illumina/strelka/tree/master/docs/userGuide 在github上面可以找到最新版地址，安装如下：

cd ~/biosoft
mkdir strelka &&  cd strelka
wget https://github.com/Illumina/strelka/releases/download/v2.8.2/strelka-2.8.2.centos5_x86_64.tar.bz2
tar xvfj strelka-2.8.2.centos5_x86_64.tar.bz2

也是需要tumor和normal的bam文件，还有参考基因组文件，就可以进行找somatic mutation啦。

reference=/home/jianmingzeng/reference/genome/human_g1k_v37/human_g1k_v37.fasta
DBSNP=/home/jianmingzeng/annotation/variation/human/dbSNP/All_20160601.vcf.gz
normal_bam=NPC10F-N_recal.bam
tumor_bam=NPC10F-T_recal.bam
sample=NPC10F 
cd ~/data/public/npc/
mkdir -p somatic/strelka
~/biosoft/strelka/strelka-2.8.2.centos5_x86_64/bin/configureStrelkaSomaticWorkflow.py \
--normalBam $normal_bam \
--tumorBam $tumor_bam \
--referenceFasta $reference \
--runDir  somatic/strelka

上面的代码只是生成了一个python脚本，还需要自行提交该脚本才能得到结果。

Successfully created workflow run script.To execute the workflow, run the following script and set appropriate options:

strelka的注意事项

对于WGS、WES和ctDNA超深度测序，需要编辑strrelka安装目录下配置文件strelka_config_bwa_default.ini

WGS：默认参数不变

WES：设置第一个参数isSkipDepthFilters=1（默认0），表示不过滤测序深度（skip depth filtration）

ctDNA：设置第一个参数isSkipDepthFilters=1，第二个参数maxInputDepth = 30000（默认10000）

其实就是需要多读软件的FAQ：https://sites.google.com/site/strelkasomaticvariantcaller/home/faq