肿瘤全外显子--记录
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 runmultiqc *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 >configcat 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.vcfdone
#!/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种,这里暂时不对此展开。
变异的注释
肿瘤全外显子--记录相关推荐
- Python命令行自动补全和记录历史命令
2019独角兽企业重金招聘Python工程师标准>>> ~$ cat .pythonstartup import os import readline import rlcomple ...
- java queue GATK_GATK 4.0 全外显子call variant
测试数据:KPGP的WES测序数据,下载地址ftp://ftp.kobic.re.kr/pub/KPGP/2017_release_candidate/WES/,分别下载了KPGP-00265,KPG ...
- unreal ue4 PixelStreaming 局域网及公有云部署全流程记录
PixelStreaming 局域网及公有云部署全流程记录 发表于 2020-04-10 | 更新于: 2020-04-20 | 分类于 Unreal Engine | 717 写在前面 本篇是 ...
- jqgrid for asp.net 单页全选记录ID
官网给的例子里单页全选得不到ID,一个一个选能得到,所以我要添加JS方法把rowid存到一个hidden里以便让后台也能收到,使全选时能存储ID. 选中状态的方法为.setSelection(rowi ...
- 利用python实现深度学习生成对抗样本模型,为任一图片加扰动并恢复原像素的全流程记录
利用python实现深度学习生成对抗样本,为任一图片加扰动并恢复原像素 一.前言 (一)什么是深度学习 (二)什么是样本模型 (三)什么是对抗样本 1.对抗的目的 2.谁来对抗? 3.对抗的敌人是谁? ...
- CentOS7+CDH5.14.0安装全流程记录,图文详解全程实测-8CDH5安装和集群配置
Cloudera Manager Server和Agent都启动以后,就可以进行CDH5的安装配置了. 准备文件 从 http://archive.cloudera.com/cdh5/parcels ...
- WeX5 3.8开发工具之蓝牙打印(全流程记录不是最全,胜似最全)
技术分享 记录踩过的坑和别的大佬没有叽歪的点 开始做蓝牙打印看到网上很多做混合开发 和安卓开发,ios开发的例子,插件等等版本大同小异,并不是像其他博客上所说的那么简单,[下载插件,无需改动,连接打印 ...
- 推荐一款高引超6000次的全基因组/全外显子组变异注释工具
SnpEff是韦恩州立大学Douglas M. Ruden团队于2012年发表的一款变异注释工具.到现在已近10年历史,持续更新,已至5.0版,总引用6044次. SnpEff的优点 基于Java环境 ...
- object怎么转list_PaddleOCR识别模型转Pytorch全流程记录
这篇文章主要负责记录自己在转PaddleOCR 模型过程中遇到的问题,以供大家参考. 重要的话说在最前面,以免大家不往下看: 本篇文章是把 "整个" ppocr 模型 转成了 py ...
- mycafe目前服务器正在维护,【图片】咖啡厅全剧情记录(修正重发)【mycafe吧】_百度贴吧...
该楼层疑似违规已被系统折叠 隐藏此楼查看此楼 Level 9 简介: Mary:她还在与Bill约会.这对年轻情侣的关系已经发展到同居了. Bill:他还在与Mary约会.他忍受着跟她一起生活,还英勇 ...
最新文章
- mysql group by using filesort优化
- romfs, cramfs和ramdisk
- 什么是Knative
- java设计模式之行为型设计模式
- C语言 后面,c语言++放在前面和后面的区别分析
- VLC支持的视频和音频文件扩展名
- 金蝶k3系统中间服务器不可用,【金蝶软件】客户端登陆时提示远程服务器不存在或不可用(金蝶K3系统)...
- Chrome浏览器误删书签恢复
- 《管理的常识》7-“什么是计划”读后感及读书笔记
- word如何插入和删除脚注,尾注
- 电驴显示没有连接服务器,电驴未连接到服务器
- Gazebo模型下载
- SyntaxError: Non-UTF-8 code starting with '\xb5' in file“问题解决办法
- 下单账户与实付账户不一致_如何保护您的不一致帐户
- 终极 Shell——ZSH
- css的`class`选择器选择前缀.
- Resin 3.0.14 和 IIS6 整合
- 硬盘分区备忘(主分区,扩展分区和逻辑分区)以及Linux硬盘分区工具parted 介绍
- 计算机网络第1章(概述)- 湖科大计算机网络课程笔记整理
- Excel如何给数字拼接加双引号或者加单引号加逗号