怎么用bioconductor的annotation包,批量找基因注释
发布网友
发布时间:2022-04-22 10:23
我来回答
共1个回答
热心网友
时间:2023-10-10 15:08
安装工作流程所需的 biconctor 包
VariantAnnotation包能够有效的从Variant Calling Format(VCF)文件读取部分或所有内容。
这些文本文件包括元信息行(meta-information lines),标题行(header line)和数据行(data lines),其中数据行每一行都含有基因组位置信息。这类格式同样包含每个位置上样本的基因型信息。更多该文件相关的信息可以看 VCF specs
本文所介绍的工作流程需要一些Biocontor的包,下面几节会仔细介绍每个包的具体用法。
可以用 biocLite 安装那些未安装的包
本工作流程着眼于17号染色体上Transient Receptor Potential Vanilloid (TRPV)基因家族的变异位点。样本数据来自于Bioconctor的cgdv17实验数据包,内部包含46个17号染色体上的完整的基因组多样性面板数据(pannel data).如果想知道这些数据是如何组织的信息,可以查看包的小品文。
我们所使用的包中的VCF文件,是CEU群体其中一个17号染色体的子集。
为了大致了解该文件有哪些数据,我们可以查看标题部分。 scanVcfHeader() 解析文件的标题部分,将解析的内容存入 VCFHeader 对象,然后就可以使用 info() 和 geno() 存取器(accessor)提取字段特定(field-specific)数据.
由下可知,VCF中的变异比对到NCBI构建的基因组GRCh37.
使用 org.Hs.eg.db 包将基因符号转为基因ID。
我们使用USCS的hg19已知基因轨道(hg19 known gene track)识别TRPV基因范围。这些基因范围最终会根据VCF文件提取变异位点。
载入注释包
我们的VCF已经比对到NCBI的基因组,并且已知基因轨道来自于UCSC。这些机构对染色体有不同的命名传统。为了在匹配(match)或者重叠(overlap)操作用到这些数据,染色体命名方式(或者叫seqlevels)需要匹配。我们会修改txdb以匹配VCF文件.
根据基因创建转录本列表
为TRPV基因创建基因范围
ScanVcfParam 对象用于提取数据子集。该对象能够指定基因组坐标(范围)或单独的VCF元素。提取范围(vs 字段)需要一个tabix索引。使用 ?indexTabix 查看细节。
locateVariants 根据基因结构(例如exon, utr, splice site等)判断变异位点的位置。我们使用之前加载的 TxDb.Hsapiens.UCSC.hg19.knownGene 包内的基因模型。
CDS的每一行都代表一个变异位点-转录本匹配,因此一行变异位点对应多行也是可以的。如果我们对基因中心的问题感兴趣,数据就可以根据基因进行描述性分析,而不用考虑转录本。
可用 predictCoding 函数得到非同义变异的氨基酸改变。 BSgenome.Hsapiens.UCSC.hg19 包用作参考等位基因的源。变异的等位基因由使用者提供。
predictCoding 仅仅返回编码变异位点的结果。与 locateVariants 一样,每个变异位点-转录匹配项的输出都有一行,因此每个变异位点可以有多行。
当 predictCoding 调用时,变异位点“not translated”在抛出的警告进行说明。 在varAllele中缺少varAllele或“N”的变异位点不会被翻译。 如果varAllele替换已经导致了移位,则后果将是“frameshift”。 有关详细信息,请参阅 ?predictCoding
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:
ensemblVEP的 file 参数必须是硬盘上的VCF