多位点取样肿瘤全外显子分析

首先介绍一个之前发表的肿瘤异质性数据库，希望能帮助此方向的相关研究人员。http://cailab.labshare.cn/cancertracer/index.html

文章 https://pubmed.ncbi.nlm.nih.gov/31701131/

肿瘤内异质性可能促进肿瘤的进化，是导致肿瘤治疗失败和耐药的关键因素。肿瘤内异质性的发展类似于生长中的树木，树枝和枝条的结构越来越复杂。肿瘤的普遍突变存在于每个肿瘤亚克隆和区域，可以用主干来表示。分支是在不同地理位置上分离的肿瘤区域或亚克隆，在单个活检中携带异质性突变，但不是在每个肿瘤区域都存在

CancerTracer数据库主要提供两类数据，分别是

1. 肿瘤内或转移内异质性:原发肿瘤或单个转移病灶内存在多个亚克隆

2. 转移后异质性:同一患者不同转移灶存在不同亚克隆

该数据库的主页面也提供了完整的教学导航如下：

如果对肿瘤异质性感兴趣的小伙伴就快跟着教学导航学起来吧~

1. 什么是外显子测序呢

利用序列捕获或者靶向技术将全基因组外显子区域DNA富集后再进行高通量测序的基因组分析。外显子组包含约1%的基因组（约30MB），却包含约85%致病突变，与个体表型相关的大部分功能性变异也都集中在染色体的外显子区。能够识别单核苷酸变异体(SNVs)、小插入缺失(InDels)以及能够解释复杂遗传疾病的罕见的原发性突变。

2. 外显子捕获试剂盒

Illumina、Agilent、BGI、罗氏NimbleGen四家外显子捕获试剂。NimbleGen和Illumina使用的是DNA探针；Agilent和BGI使用的是RNA探针。

3. 全外显子测序工作流程

首先将基因组DNA打断成200-300bp,然后末端修复之后加A加接头，LM-PCR的线性扩增之后使用指定的捕获试剂kit对目标片段进行杂交捕获,再经过LM-PCR的线性扩增后Q-PCR定量，如果文库检测合格后就可上机测序，接下来就是数据下机后的信息分析

4. 介绍外显子捕获效率

外显子测序过程中要用到杂交。在人的染色体上有许多与外显子有同源性的部分，在杂交过程中也被捕获。由于是随机打断，可能一条reads上面有外显子区域也有侧翼区域会被一起抓下，就会使测到的序列中有一部分不是外显子序列。把比对到外显子的序列占全部测序序列的比列称为捕获效率。reads捕获效率一般在70%-80%，碱基的捕获效率一般在50%-70%。

5. 外显子测序一般建议做多少倍的覆盖

一般做100X或者150X。较高的覆盖倍数，对于测异质性的遗传变异可以发现小比例的突变。外显子测序的覆盖不是特别均匀，这样较高的平均覆盖率有利于保证大部分的区域有足够的覆盖倍数。100X的平均深度下，至少有90%的区域覆盖度可以达到10X以上。

6. 外显子捕获可以像全基因组测序那样做CNV

外显子测序因为有一个杂交捕获的过程，这里就会有一个杂交效率的问题。各个外显子的杂交效率是不同的，其同源竞争的情况也不同，所以不同的外显子的覆盖率的差异就很大。所以一般外显子测序不能用于CNV的检测。但在癌症研究中，利用癌组织和癌旁（或者血液样品）进行对照，有方法可以检测CNV。

7. 外显子测序特有优点

外显子测序是全基因重测序的一个较为经济的替代手段，对研究基因的SNP、InDel等具有较大的优势，但无法研究基因组结构变异如染色体断裂重组等，一般在疾病研究中，可结合转录组测序一起研究。
（1）人的全基因是3G，一般要平均30x的覆盖，也就是要超过90G的数据，如果样本多的话，费用负担不起。外显子只占人全部基因序列的1%，而且是最关键的1%。
（2）肿瘤本身存在着较大异质性，且往往肿瘤样品的纯度不高。肿瘤的基因序列是不稳定的，一直在变的，也就是肿瘤的深部、浅部可能其基因序列是不一样的。为了测出各种突变(低频)，就需要有较深的测序深度，比如100x甚至200x的测序深度。这时候外显子测序就可以做到高的测序覆盖度，同时费用不会太高。

8. 什么是遗传变异

生物信息学中各种基因组研究的基础就是遗传变异的研究，比如进化和各种表型的研究。遗传变异包括单核苷酸多态性(SNP)，小片段的插入缺失(Indel)，结构变异(SV)，拷贝数变异(CNV)等等。区分这些变异简单的方法就是变了几个，其中SNP是单个核苷酸的改变，indel通常是50bp以下的变异，SV和CNV则要更长。

9. SNV 和 SNP

SNP 和 SNV 都是单碱基的突变，但是SNP 是多了一个频率属性的SNV，比如在群体中1%以上

10. biallelic and multiallelic

biallelic 表示在基因组的某个位点上有两个等位基因，存在一个和参考基因组相同的碱基和一个和参考基因组不同的碱基或者是一个deletion。 multiallelic 多等位基因表示在基因组某个位点可以观测到三个或者多个等位基因，在vcf文件中可以看到两个或者三个非参考基因组的突变。多等位基因并不常见.

11. Transition vs Transversion

转换(transition)则是嘌呤被嘌呤，或嘧啶被嘧啶替代，颠换(transversion)是指嘌呤与嘧啶的变化。

12. SNP 种类

全基因组SNP突变可以分成6类（C>A, C>G, C>T, A>C, A>G, A>T）。以A:T>C:G为例，此种类型SNP突变包括A>C和T>G。由于测序数据既可比对到参考基因组的正链，也可比对到参考基因组的负链，当T>C类型突变出现在参考基因组正链上，A>G类型突变即在参考基因组负链的相同位置，所以通常也会将T>C和A>G划分成一类。

13. SNP的可能影响

编码区SNP，根据密码子简并性 SNP 不一定会引起编码氨基酸的改变，不引起任何变化的叫做Synonymous SNP， 而引起氨基酸变化的叫做Non-Synonymous SNP。如果编码的某种氨基酸的密码子变成另一种，会导致多肽链的氨基酸种类和序列发生改变，这就是一个错义突变。当突变使一个编码氨基酸的密码子变成终止子时，则蛋白质合成进行到该突变位点时会提前终止，这时就是无义突变。

14. 同一份BAM文件samtools depth和samtools mpileup的结果不同？

Single-End的数据差异应该比较小，不会太明显，Pair-End上差异可能会比较大，原因有两点：
（1）mpileup会默认将PE reads中比对情况异常的那些reads（包括双端都比上，但是两条配对reads之间的比对距离明显偏离了插入片段的长度分布，或者一端比对上而另一端没比对上）丢弃，则在计算depth时，这些被丢弃的reads不参与统计，而depth则不会，depth默认不做任何过滤 。mpileup中提供了一个-A选项来保留异常比对的reads，如果设置了这个选项，那么在这点上mpileup的depth统计就和samtools depth相同了
（2）默认mpileup还会过滤掉测序质量值低于13的碱基，depth默认不过滤
从上面列出的两点差异可以看出，mpileup默认输出的是高质量的覆盖深度，mpileup功能被开发出来是为了与bcftools组合，将samtools mpileup的输出作为bcftools的输入用于下游的snp-calling，当然需要保证数据的质量

当然可以通过设置对应的参数使得它的属于结果与depth的一致，但是不推荐这么做

---

比对结果的质控

• flagstat + multiqc
  #flagstat.sh的脚本
#!/bin/bash
ls *bam >name
cat name |while read id
do
samtools flagstat $id > $id.stat
done
  then run multiqc *stat
  可以得到所有bam文件的统计信息
• qualimap
  qualimap使用外显子坐标的bed文件
  用qualimap的bamqc统计比对bam文件的深度和覆盖度
  #qualimap.sh的内容
#!/bin/bash
exon_bed=.../hg38.exon.bed
out_dir=./QC/qualimap
ls *.bam >1
cat 1 |while read id
do
qualimap bamqc --java-mem-size=20G -gff $exon_bed -bam $id -outdir $out_dir
mv $out_dir/genome_results.txt $out_dir/$id.txt
done
  qualimap的结果包括了mean mapping quality，这是flagstat和后面的gatk没有包含的而这个mapping quality的计算方法和碱基质量值一样，Q=-10loge
• CollectHsMetrics?

---

变异质控和过滤

(1)软过滤

使用GATK的VQSR，使用机器学习的方法，利用多个不同数据的特征训练一个模型（高斯混合模型）对变异的结果进行质控
但是软过滤是有一定的条件的：
（1）需要一个精心准备的已知变异集合，作为训练质控模型的真集，比如：Hapmap，OMINI，1000G，dbsnp等这些国际性项目的数据
（2）要求新检测的结果中有足够多的变异，不然VQSR在进行模型训练时，会因为可用的变异位点数目不足而无法进行。
很多非人的物种在完成变异检测后，没法使用VQSR进行质控。同样的，一些小panel、外显子测序，由于最后的变异位点数目不够而无法进行VQSR。
全基因组分析或者多个样本的全外显子组分析可以使用此法。

(2)硬过滤

对一些过滤指标一定的阈值，然后一刀切,主要涉及以下的过滤指标：
（1）QualByDepth（QD）：Variant confidence normalized by unfiltered depth of variant samples (QD)。QD是变异质量值除以覆盖深度得到的比值这里的变异质量值就是VCF中QUAL的值——用来衡量变异的可靠程度。这里的覆盖深度是这个位点上所有含有变异碱基的样本的覆盖深度之和，通俗一点说，就是这个值可以通过累加每个含变异的样本（GT为非0/0的样本）的覆盖深度（VCF中每个样本里面的DP）而得到。
QD这个值描述的实际上就是单位深度的变异质量值，也可以理解为是对变异质量值的一个归一化，QD越高一般来说变异的可信度也越高。在质控的时候，相比于QUAL或者DP（深度）来说，QD是一个更加合理的值。因为我们知道，原始的变异质量值实际上与覆盖的read数目是密切相关的，深度越高的位点QUAL一般都是越高的，而任何一个测序数据，都不可避免地会存在局部深度不均的情况，直接使用QUAL或者DP都会很容易因为覆盖深度的差异而带来有偏的质控结果。

（2）FisherStrand（FS）：Strand bias estimated using Fisher’s exact test (FS)。FS是一种通过Fisher检验得到的pValue转换而来的值，它描述的是在测序或者比对的时候，对只含有变异的read以及只含有参考序列碱基的read是否有明显的正负链特异性（strand bias，或者说是差异性）。这个差异反应了测序过程不够随机，或者是比对算法在基因组的某些区域存在一定的选择偏向。如果测序过程是随机的，比对是没问题的，那么不管read是否含有变异，以及是否来自基因组的正链或者负链，只要是真实的它们就都应该是比较均匀的，不会出现链特异的比对结果，FS应该接近于零。

（3）StrandOddsRatio (SOR):Strand bias estimated by the symmetric odds ratio test (SOR)。同样是对链特异性的一种描述。由于很多时候read在外显子区域末端的覆盖存在着一定的链特异（这个区域的现象其实是正常的），往往只有一个方向的read，这个时候该区域中如果有变异位点的话，那么FS通常会给出很差的分值，这时SOR就能够起到比较好的校正作用了。
StrandBiasBySample：Number of forward and reverse reads that support REF and ALT alleles (SB)。

（4）RMSMappingQuality（MQ）:Root mean square of the mapping quality of reads across all samples (MQ)。MQ这个值是所有比对到该位点上的read的比对质量值的均方根(先平方、再平均、然后开方)。它和平均值相比更能够准确地描述比对质量值的离散程度。

(5)MappingQualityRankSumTest（MQRankSum）：Rank sum test for mapping qualities of REF versus ALT reads (MQRankSum)。

(6)ReadPosRankSumTest (ReadPosRankSum)：Rank sum test for relative positioning of REF versus ALT alleles within reads (ReadPosRankSum)

上诉的这些指标一般是用来进行硬过滤常用的指标现在的问题是，这些指标设置怎样的值才算是合理的，GATK给的推荐一般是：
SNP：
QD<2.0
MQ<40.0
FS>60.0
SOR>3.0
MQRankSum <-12.5
ReadPosRankSum<-8.0
INDEL:
QD<2.0
FS>200.0
SOR>10.0
MQRankSum <-12.5
ReadPosRankSum<-8.0

这里需要说明的是包含在这些指标内的变异是被判断不好的变异，会被过滤掉

#gatk_filter.sh的内容
#!/bin/bash
mutation=./mutation
vcf_clean=./vcf_clean
cd $mutation
ls *.vcf >1
cat 1|while read id
do
gatk SelectVariants -select-type SNP -V $id -O $vcf_clean/${id}_snp.vcf
gatk SelectVariants -select-type INDEL -V $id -O $vcf_clean/${id}_indel.vcf
gatk VariantFiltration -V $vcf_clean/${id}_snp.vcf --filter-expression "QD<2.0 || MQ<40.0 || FS>60.0 || SOR>3.0 || MQRankSum <-12.5 || ReadPosRankSum<-8.0" --filter-name "Filter" -O $vcf_clean/${id}_filter.snp.vcf
gatk VariantFiltration -V $vcf_clean/${id}_indel.vcf --filter-expression "QD<2.0 || FS>200.0 || SOR>10.0 || MQRankSum <-12.5 || ReadPosRankSum<-8.0" --filter-name "Filter" -O $vcf_clean/${id}_filter.indel.vcf
done
满足–filter-expression表达式的，表示质量异常，将会被标记成Filter
剩下的表示通过，将会自动标记为PASS

变异结果的评估

(1) 可以通过 GATK 的VariantEval工具获得每个样本的 Ti/Tv ratio，以此来进一步确定结果的可靠程度。

Ti/Tv是一个被动指标，它是对最后质控结果的一个反应，我们是不能够在一开始的时候使用这个值来进行变异过滤的。
如果没有选择压力的存在，Ti/Tv将等于0.5，因为从概率上讲Tv将是Ti的两倍。但现实当然不是这样的，比如对于人来说，全基因组正常的Ti/Tv在2.1左右，而外显子区域是3.0左右，新发的变异（Novel variants）则在1.5左右。

(2) 计算杂合变异和纯和变异的比例

使用gatk的CollectVariantCallingMetrics功能
#collect.sh的内容
#!/bin/bash
resource=/home/ubuntu/WGS_new/output/bwa_mapping/GATK/resource
ls *filter*vcf >1
cat tmp.txt|while read id
do
echo $name
gatk CollectVariantCallingMetrics --DBSNP $resource/dbsnp_146.hg38.vcf.gz -I $id -O $id
done
将生成的文件合并
#count.sh的内容
#!/bin/bash
ls *indel*detail_metrics >1
ls *snp*detail_metrics >2
cat 1 |cut -d "." -f 1 >3
paste 1 2 3 >config

cat config |while read name
do
arr=($name)
indel=${arr[0]}
snp=${arr[1]}
id=${arr[2]}

grep -v "#" $indel > ${id}_indel_detail
rep -v "#" $snp > ${id}_snp_detail

```done``

ls *indel_detail >4
ls *snp_detail >5

paste 4 5 >config2

cat config2 |while read name
do
arr=($name)
indel=${arr[0]}
snp=${arr[1]}

cut -f 1,2,13,14,15 $indel > ${indel}.txt

cut -f 1,2,6,7,8,11,12 $snp >${snp}.txt
done
cat *indel*txt > final_metrics_indel_tmp.txt
cat *snp*txt > final_metrics_snp_tmp.txt

sed -i "/^$/d" final_metrics_indel_tmp.txt
sed -i "/^$/d" final_metrics_snp_tmp.txt

awk '{if (NR==1) {print $0;} else {if (NR%2==0) print $0;}}' final_metrics_indel_tmp.txt >final_metrics_indel.txt
awk '{if (NR==1) {print $0;} else {if (NR%2==0) print $0;}}' final_metrics_snp_tmp.txt >final_metrics_snp.txt

文件列合并用paste，行合并用cat

这里的dnsbp_Ti/Tv指的是包含在dbsnp数据库中的Ti/Tv比例
Novel_Ti/Tv指的是不包含在dbsnp中的新的snp的Ti/Tv的比例
总的来说似乎离3有一定偏差
杂合/纯和比例似乎是比较低，有点奇怪

制作突变频谱

一般可以分为6类
C>A、C>T、C>G 、A>T 、A>C 、A>G
在这之前，先把PASS过的vcf文件提取出来，顺便indel一起做了

#extract.sh 的内容
#!/bin/bash
ls *filter.indel.vcf >1
ls *filter.snp.vcf >2
sed 's/\./\t/g' 1 |cut -f 1 >3
paste 1 2 3 >config

cat config |while read id 
do
arr=($id)
indel=${arr[0]}
snp=${arr[1]}
name=${arr[2]}
grep "PASS" $indel > ${name}_indel_PASS.vcf
grep "PASS" $snp > ${name}_snp_PASS.vcf

done

#!/bin/bash
ls *snp*PASS* >1
cat 1 |while read id
do
cut -f 4,5 $id|sort |uniq -c |grep -v "," >${id}.spectrum.txt
done

先sort再uniq -c 是因为sort之后相同的部分会排到一起，这时uniq -c才能算出数量
以NPC10F-N和NPC10F-T为例

library(ggplot2)
install.packages("devtools")
library(devtools)
install_github("easyGgplot2", "kassambara")
library(easyGgplot2)


a = read.table("NPC10F-T_snp_PASS.vcf.spectrum.txt")
dat <- data.frame(type=c('C>A(G>T)','C>T(G>A)','C>G(G>C)','T>A(A>T)','T>G(A>C)','T>C(A>G)'))
counts1=c(a[4,1]+a[9,1],a[6,1]+a[7,1],a[5,1]+a[8,1],a[10,1]+a[3,1],a[12,1]+a[1,1],+a[11,1]+a[2,1])

p1=ggplot(dat,aes( x=type,y=counts1))+geom_bar(stat="identity",aes(fill=type))+labs(title="NPC10F-T")+theme(plot.title = element_text(hjust = 0.5),axis.text.x = element_text(size=10,vjust = 0.5, hjust = 0.5,face="bold",angle = 75)) 

b= read.table("NPC10F-N_snp_PASS.vcf.spectrum.txt")
counts2=c(b[4,1]+b[9,1],b[6,1]+b[7,1],b[5,1]+b[8,1],b[10,1]+b[3,1],b[12,1]+b[1,1],+b[11,1]+b[2,1])
p2=ggplot(dat,aes( x=type,y=counts2))+geom_bar(stat="identity",aes(fill=type))+labs(title="NPC10F-N")+theme(plot.title = element_text(hjust = 0.5),,axis.text.x = element_text(size=10,vjust = 0.5, hjust = 0.5,face="bold",angle = 75)) 

ggplot2.multiplot(p1,p2, cols=2)

似乎很明显能发现，C>T(G>A)和T>C(A>G)突变的频率最高
这种其实就是所谓的转换（transition），即嘧啶到嘧啶，嘌呤到嘌呤
这也正好验证了之前计算的Ti/Tv比例，之前计算的大致在2~3之间，正好说明这种突变频率高，前后得出结论一致。

正常情况下SNV突变类型就6种，但是如果结合突变的上下文就可以变成96种，这里暂时不对此展开。

变异的注释