基因组变异位点注释

安装工作流程所需的biconductor包

source("http://bioconductor.org/workflows.R")
workflowInstall("variants")

背景

VariantAnnotation包能够有效的从Variant Calling Format(VCF)文件读取部分或所有内容。
这些文本文件包括元信息行(meta-information lines)，标题行(header line)和数据行(data lines)，其中数据行每一行都含有基因组位置信息。这类格式同样包含每个位置上样本的基因型信息。更多该文件相关的信息可以看VCF specs

配置

本文所介绍的工作流程需要一些Biocondutor的包，下面几节会仔细介绍每个包的具体用法。

library(VariantAnnotation)
library(cgdv17)
library(org.Hs.eg.db)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
library(BSgenome.Hsapiens.UCSC.hg19)
library(PolyPhen.Hsapiens.dbSNP131)

可以用biocLite安装那些未安装的包

source("https://bioconductor.org/biocLite.R")
biocLite("mypackage")

探索TRPV基因家族的变异位点

本工作流程着眼于17号染色体上Transient Receptor Potential Vanilloid (TRPV)基因家族的变异位点。样本数据来自于Bioconductor的cgdv17实验数据包，内部包含46个17号染色体上的完整的基因组多样性面板数据（pannel data).如果想知道这些数据是如何组织的信息，可以查看包的小品文。

browseVignettes("cgdv17")

我们所使用的包中的VCF文件，是CEU群体其中一个17号染色体的子集。

library(VariantAnnotation)
library(cgdv17)
file <- system.file("vcf", "NA06985_17.vcf.gz", package = "cgdv17")

检查VCF文件的标题数据

为了大致了解该文件有哪些数据，我们可以查看标题部分。scanVcfHeader()解析文件的标题部分，将解析的内容存入VCFHeader对象，然后就可以使用info()和geno()存取器（accessor）提取字段特定(field-specific)数据.

hdr <- scanVcfHeader(file)
info(hdr)
geno(hdr)

由下可知，VCF中的变异比对到NCBI构建的基因组GRCh37.

meta(hdr)$META
## DataFrame with 5 rows and 1 column
##                       Value
##                 <character>
## fileformat          VCFv4.1
## fileDate           20111102
## source     masterVar2VCFv40
## reference    build37.fa.bz2
## phasing             partial

将基因符号(symbol)变为基因ID

使用org.Hs.eg.db包将基因符号转为基因ID。

library(org.Hs.eg.db)
## 定义哪些基因
genesym <- c("TRPV1", "TRPV2", "TRPV3")
## 从org.Hs.eg.db挑选这些符号的ID,select如何使用见之前文章
geneid <- select(org.Hs.eg.db, keys=genesym, keytype="SYMBOL",
                 columns="ENTREZID")
geneid
##   SYMBOL ENTREZID
## 1  TRPV1     7442
## 2  TRPV2    51393
## 3  TRPV3   162514

创建基因范围(gene ranges)

我们使用USCS的hg19已知基因轨道（hg19 known gene track)识别TRPV基因范围。这些基因范围最终会根据VCF文件提取变异位点。
载入注释包

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
txdb

我们的VCF已经比对到NCBI的基因组，并且已知基因轨道来自于UCSC。这些机构对染色体有不同的命名传统。为了在匹配（match)或者重叠(overlap)操作用到这些数据,染色体命名方式（或者叫seqlevels)需要匹配。我们会修改txdb以匹配VCF文件.

txdb <- renameSeqlevels(txdb, gsub("chr", "", seqlevels(txdb)))
txdb <- keepSeqlevels(txdb, "17")

根据基因创建转录本列表

txbygene = transcriptsBy(txdb, "gene")

为TRPV基因创建基因范围

gnrng <- unlist(range(txbygene[geneid$ENTREZID]), use.names=FALSE)
names(gnrng) <- geneid$SYMBOL

提取变异位点子集

ScanVcfParam对象用于提取数据子集。该对象能够指定基因组坐标（范围）或单独的VCF元素。提取范围(vs 字段)需要一个tabix索引。使用?indexTabix查看细节。

param <- ScanVcfParam(which = gnrng, info = "DP", geno = c("GT", "cPd"))
param
## 从VCF文件提取TRPV范围
vcf <- readVcf(file, "hg19", param)
## 查看VCF对象 'fixed', 'info' and 'geno'
vcf
head(fixed(vcf))
geno(vcf)

基因模型的变异位点位置

locateVariants根据基因结构（例如exon, utr, splice site等）判断变异位点的位置。我们使用之前加载的TxDb.Hsapiens.UCSC.hg19.knownGene包内的基因模型。

## 使用'region'参数定义目标区域
## 详细内容见？locateVariants
cds <- locateVariants(vcf, txdb, CodingVariants())
five <- locateVariants(vcf, txdb, FiveUTRVariants())
splice <- locateVariants(vcf, txdb, SpliceSiteVariants())
intron <- locateVariants(vcf, txdb, IntronVariants())

all <- locateVariants(vcf, txdb, AllVariants())

CDS的每一行都代表一个变异位点-转录本匹配，因此一行变异位点对应多行也是可以的。如果我们对基因中心的问题感兴趣，数据就可以根据基因进行描述性分析，而不用考虑转录本。

## variants是一对多关系么
table(sapply(split(mcols(all)$GENEID, mcols(all)$QUERYID), 
      function(x) length(unique(x)) > 1))

## 根据基因分析变异位点数量
idx <- sapply(split(mcols(all)$QUERYID, mcols(all)$GENEID), unique)
sapply(idx, length)


## 根据基因分析变异位点位置
sapply(names(idx),
    function(nm) {
        d <- all[mcols(all)$GENEID %in% nm, c("QUERYID", "LOCATION")]
        table(mcols(d)$LOCATION[duplicated(d) == FALSE])
    })

非同义突变位点的氨基酸编码改变

可用predictCoding函数得到非同义变异的氨基酸改变。BSgenome.Hsapiens.UCSC.hg19包用作参考等位基因的源。变异的等位基因由使用者提供。

library(BSgenome.Hsapiens.UCSC.hg19)
seqlevelsStyle(vcf) <- "UCSC"
seqlevelsStyle(txdb) <- "UCSC"
aa <- predictCoding(vcf, txdb, Hsapiens)

predictCoding仅仅返回编码变异位点的结果。与locateVariants一样，每个变异位点-转录匹配项的输出都有一行，因此每个变异位点可以有多行。

## 是否有多个突变位点位于同一基因
table(sapply(split(mcols(aa)$GENEID, mcols(aa)$QUERYID), 
        function(x) length(unique(x)) > 1))

## FALSE
##    17

## 根据基因总结变异位点数量
idx <- sapply(split(mcols(aa)$QUERYID, mcols(aa)$GENEID, drop=TRUE), unique)
sapply(idx, length)

## 162514  51393   7442
##      6      3      8

## 根据基因总结变异位点后果
sapply(names(idx),
       function(nm) {
           d <- aa[mcols(aa)$GENEID %in% nm, c("QUERYID","CONSEQUENCE")]
           table(mcols(d)$CONSEQUENCE[duplicated(d) == FALSE])
       })

##                162514 51393 7442
## nonsynonymous       2     0    2
## not translated      1     0    5
## synonymous          3     3    1

当predictCoding调用时,变异位点“not translated”在抛出的警告进行说明。 在varAllele中缺少varAllele或“N”的变异位点不会被翻译。 如果varAllele替换已经导致了移位，则后果将是“frameshift”。 有关详细信息，请参阅?predictCoding

使用ensemblVEP包进行注释

ensemblVEP包能够访问在线Ensembl Variant Effect Predictor (VEP tool)。VEP工具输出的已知或者未知变异位点的功能后果预测，通过序列本体论(Sequence Ontology)或Ensembl报告。可选输出有Regulatory region consequences, HGNC, Ensembl protein identifiers, HGVS, co-located variants。 ensemblVEP()接受VCF文件名，在R工作环境中返回一个磁盘上的VCF或者GRanges.
加载ensemblVEP:

library(ensemblVEP)

ensemblVEP的file参数必须是硬盘上的VCF文件。我们之后会以TRPV变异位点写出VCF对象，并且上传到ensemblVEP.
以tempfile()写出一个VCF对象到vcf文件中

dest <- tempfile()
writeVcf(vcf, dest)

使用只含有TRPV变异位点的文件和自定义VERPParam对象调用ensemblVEPCall ensemblVEP

gr <- ensemblVEP(file = dest)
head(gr, 3)

来自VEP工具的数据作为元数据列返回给GRanges（'gr'对象）。 您可以通过设置VEPParam中的运行时选项来进一步控制返回哪些字段。

VEPParam()

每个选项，例如“input”，“basic”，“cache”等具有可以设置的多个字段。

basic(VEPParam())

更多运行选项的信息可以在VEPParam的man page或者
Ensembl VEP web site中找到.

了解更多信息

想要更多了解VariantAnnotation的使用方法，有以下方法：

help(package="VariantAnnotation")
?predictCoding
browseVignettes(package="VariantAnnotation")

顺便放一下自己的知识星球，如果你觉得我对你有帮助的话。

知识星球