05

用GEMINI来探索vcf格式的突变数据

第一次听说这个软件,是一个香港朋友推荐的:http://davetang.org/muse/2016/01/13/getting-started-with-gemini/ 他写的很棒,但是我当初以为是一个类似于SQLite的数据库浏览模式,所以没在意。实际上,我现在仍然觉得这个软件没什么用!

软件官网有详细的介绍:https://gemini.readthedocs.io/en/latest/

而且提供丰富的教程:

We recommend that you follow these tutorials in order, as they introduce concepts that build upon one another.

  • Introduction to GEMINI, basic variant querying and data exploration. html pdf
  • Identifying de novo mutations underlying Mendelian disease html pdf
  • Identifying autosomal recessive variants underlying Mendelian disease html pdf
  • Identifying autosomal dominant variants underlying Mendelian disease html pdf
  • Other GEMINI tools html pdf

软件本身并不提供注释,虽然它的功能的确包括注释,号称可以利用(ENCODE tracks, UCSC tracks, OMIM, dbSNP, KEGG, and HPRD.)对你的突变位点注释,比如你输入1       861389  .       C       T       ,它告诉你这个突变发生在哪个基因,对蛋白改变如何?是否会产生某些疾病?

虽然它本身没有注释功能,但是它会调用snpEFF或者VEP进行注释,你需要自己先学习它们。

1

软件安装:

GEMINI是用python写的,有一个小脚本可以自动完成安装过程:

7.3K May  4 14:44 gemini_install.py

下载这个脚本,然后安装即可

wget https://github.com/arq5x/gemini/raw/master/gemini/scripts/gemini_install.py

python gemini_install.py $tools $data

PATH=$tools/bin:$data/anaconda/bin:$PATH

where $tools and $data are paths writable on your system.

我把$tools用的就是当前文件夹,$data也是当前文件夹下面的gemini文件夹。

这样就会在当前文件夹下面生成两个文件夹,bin是存储程序,gemini是存储数据用的,而且注意要把bin目录的全路径添加到环境变量!

输入数据:

我们可以直接下载软件作者提供的测试数据

首先是22号染色体的所有突变位点经过WEP注释的文件

然后是一个三口直接的突变ped格式数据

数据存放在亚马逊云,所有的教程pdf也在

http://s3.amazonaws.com/gemini-tutorials

如果是你自己的vcf文件,需要自己用VEP注释一下

1

运行命令:

2

结果解读:

产生是chr22.db就是一个数据库格式的文件,但是需要用gemini 来进行查询,个人认为,并没有多大意思!

你只要熟悉mySQL等SQL语言,完全可以自己来!

05

用VEP对vcf格式的突变数据进行注释

VEP是国际三大数据库之一的ENSEMBL提供的,也是非常主流和方便,但它是基于perl语言的,所以在模块方面可能会有点烦人。跟snpEFF一样,也是对遗传变异信息提供更具体的注释,而不仅仅是基于位点区域和基因。如果你熟悉外显子联盟这个数据库EXAC(ExAC.r0.3.sites.vep.vcf.gz),你可以下载它所有的突变记录数据,看看它对每个变异位点到底注释了些什么,它就是典型的用VEP来注释的。 Continue reading

05

用snpEFF对vcf格式的突变数据进行注释

这个软件比较重要,尤其是对做遗传变异相关研究的,很多人做完了snp-calling后喜欢用ANNOVAR来进行注释,但是那个注释还是相对比较简单,只能得到该突变位点在基因的哪个区域,那个基因这样的信息,如果想了解更具体一点,就需要更加功能化的软件了,snpEFF就是其中的佼佼者,而且是java平台软件,非常容易使用!而且它的手册写的非常详细:http://snpeff.sourceforge.net/SnpEff_manual.html Continue reading

14

居然可以下载千人基因组计划的所有数据bam,vcf数据

它有两个ftp站点存储所有的数据!
直接看最新版的数据,共有NA编号开头的1182个人,HG开头的1768个人!
每个人的目录下面都有 四个数据文件夹
Oct 01 2014 00:00    Directory alignment
Oct 01 2014 00:00    Directory exome_alignment
Oct 01 2014 00:00    Directory high_coverage_alignment
Oct 01 2014 00:00    Directory sequence_read
这些数据实在是太丰富了!
也可以直接看最新版的vcf文件,记录了这两千多人的所有变异位点信息!
可以直接看到所有的位点,具体到每个人在该位点是否变异!
不过它的基因型信息是通过MVNcall+SHAPEIT这个程序call出来的,具体原理见:http://www.ncbi.nlm.nih.gov/pubmed/23093610
它有两个ftp站点存储所有的数据!
直接看最新版的数据,共有NA编号开头的1182个人,HG开头的1768个人!
每个人的目录下面都有 四个数据文件夹
Oct 01 2014 00:00    Directory alignment
Oct 01 2014 00:00    Directory exome_alignment
Oct 01 2014 00:00    Directory high_coverage_alignment
Oct 01 2014 00:00    Directory sequence_read
这些数据实在是太丰富了!
也可以直接看最新版的vcf文件,记录了这两千多人的所有变异位点信息!
可以直接看到所有的位点,具体到每个人在该位点是否变异!
不过它的基因型信息是通过MVNcall+SHAPEIT这个程序call出来的,具体原理见:http://www.ncbi.nlm.nih.gov/pubmed/23093610
clipboard

 

05

使用oncotator做突变注释

功能:vcf格式突变数据进一步注释成maf格式

做过癌症数据分析的童鞋都知道,TCGA里面用maf格式来记录突变!那么maf格式的数据是如何得来的呢,我们都知道,做完snp-calling一般是得到vcf格式的突变记录数据文件,然后再用annovar或者其它蛋白结构功能影响预测软件注释一下,还远达不到maf的近100条记录。

而大名鼎鼎的broad institute就规定了maf格式的突变注释文件,他就是利用了十几个常见的已知数据库来注释我们得到的vcf突变记录,通常是对somatic的突变才注释成maf格式的数据!
大名鼎鼎的broadinstitute出品的突变注释工具:http://www.ncbi.nlm.nih.gov/pubmed/25703262
本身也是一个在线工具:
集成了下面所有的分析资源
而且还提供了API

Genomic Annotations

  • Gene, transcript, and functional consequence annotations using GENCODE for hg19.
  • Reference sequence around a variant.
  • GC content around a variant.
  • Human DNA Repair Gene annotations from Wood et al.

Protein Annotations

  • Site-specific protein annotations from UniProt.
  • Functional impact predictions from dbNSFP.

Cancer Variant Annotations

Non-Cancer Variant Annotations

因为要下载的数据有点多,我这里就不用自己的电脑测试了,安装过程也很简单的!

 

17

转载-VCF格式详解

CHROM(chromosome):染色体

POS - position:参考基因组variant碱基位置,如果是INDEL(插入缺失),位置是INDEL的第一个碱基位置

ID - identifier: variant的ID。比如在dbSNP中有该SNP的id,则会在此行给出;若没有,则用’.'表示其为一个novel variant。

REF - reference base(s):参考碱基,染色体上面的碱基,必须是ATCGN中的一个,N表示不确定碱基

ALT - alternate base(s):与参考序列比较发生突变的碱基

QUAL - quality: Phred格式(Phred_scaled)的质量值,表 示在该位点存在variant的可能性;该值越高,则variant的可能性越大;计算方法:Phred值 = -10 * log (1-p) p为variant存在的概率; 通过计算公式可以看出值为10的表示错误概率为0.1,该位点为variant的概率为90%。

FILTER - _filter status: 使用上一个QUAL值来进行过滤的话,是不够的。GATK能使用其它的方法来进行过滤,过滤结果中通过则该值为”PASS”;若variant不可靠,则该项不为”PASS”或”.”。

INFO - additional information:  这一行是variant的详细信息,具体如下:

DP-read depth:样本在这个位置的reads覆盖度。是一些reads被过滤掉后的覆盖度。          DP4:高质量测序碱基,位于REF或者ALT前后

MQ:表示覆盖序列质量的均方值RMS Mapping Quality

FQphred值关于所有样本相似的可能性

AF1 AF(Allele Frequency) 表示Allele的频率,AF1为第一个ALT allele 发生频率的可能性评估

AC1AC表示Allele(等位基因)的数目,AC1为对第一个ALT allele count的最大可能性评估

AN:AN(Allele Number) 表示Allele的总数目

IS插入缺失或部分插入缺失的reads允许的最大数量

ACAC(Allele Count) 表示该Allele的数目

G3ML 评估基因型出现的频率

HWE:chi^2基于HWE的测试p值和G3

CLR在受到或者不受限制的情况下基因型出现可能性log值

UGT:最可能不受限制的三种基因型结构

CGT:最可能受限制三种基因型的结构

PV4四种P值得误差,分别是(strand、baseQ、mapQ、tail distance bias)

INDEL:表示该位置的变异是插入缺失

PC2非参考等位基因的phred(变异的可能性)值在两个分组中大小不同

PCHI2后加权chi^2,根据p值来测试两组样本之间的联系

QCHI2:Phred scaled PCHI2.

PR置换产生的一个较小的PCHI2

QBD:Quality by Depth,测序深度对质量的影响

RPB序列的误差位置(Read Position Bias)

MDV:样本中高质量非参考序列的最大数目

VDB:Variant Distance Bias,RNA序列中过滤人工拼接序列的变异误差范围

GT样品的基因型(genotype)。两个数字中间用’/'分 开,这两个数字表示双倍体的sample的基因型。0 表示样品中有ref的allele; 1 表示样品中variant的allele; 2表示有第二个variant的allele。因此: 0/0 表示sample中该位点为纯合的,和ref一致; 0/1 表示sample中该位点为杂合的,有ref和variant两个基因型; 1/1 表示sample中该位点为纯合的,和variant一致。

GQ基因型的质量值(Genotype Quality)。Phred格式(Phred_scaled)的质量值,表示在该位点该基因型存在的可能性;该值越高,则Genotype的可能性越 大;计算方法:Phred值 = -10 * log (1-p) p为基因型存在的概率。

GL三种基因型(RR RA AA)出现的可能性,R表示参考碱基,A表示变异碱基

DV高质量的非参考碱基

SPphred的p值误差线

PL:指定的三种基因型的质量值(provieds the likelihoods of the given genotypes)。这三种指定的基因型为(0/0,0/1,1/1),这三种基因型的概率总和为1。和之前不一致,该值越大,表明为该种基因型的可能 性越小。 Phred值 = -10 * log (p) p为基因型存在的概率。

FORMAT BC1-1-base.sorted.bam这两行合起来提供了’ BC1-1-base′这个sample的基因型的信息。’ BC1-1-base′代表这该名称的样品,是由BAM文件中的@RG下的 SM 标签决定的。