SeqArray概述

Xiuwen郑博士(西雅图华盛顿大学的生物统计学)

2016年6月25日

介绍

全基因组测序(WGS)数据生成正在以前所未有的速度

1000人基因工程三期(1公斤)
- 8100万变异和2504个人
- https://www.internationalgenome.org/
变体调用格式(VCF)
- 一个通用和灵活的基于文本的格式
- VCF文件大,数据检索是相对较慢

方法

CoreArray (c++库)

为全基因组变异的大规模数据管理而设计的
数据格式(GDS)来存储多个array-oriented数据集在一个单一的文件中

两个R包

gdsfmt - R接口CoreArray基因组数据结构(GDS)文件
SeqArray——专门为数据管理的全基因组序列变异从变量调用格式(VCF)文件

方法——优势

作为VCF SeqArray提供了相同的功能
将数据存储在一个二进制和array-oriented方式
- 使用R语言有效的访问
基因型存储在一个压缩的方式
- 2比特的数组来存储等位基因(网站bi-allelic 95%)
- 罕见变异:高度压缩
- 1公斤,2035亿基因型,保存在4.3 g如果一个字节存储基因型(2.26%)
并行访问
- 多个集群节点和/或核

方法——文件内容

文件:SeqArray / extdata / CEU_Exon。gds (387.3 k) | - +描述[]* | +样品。id {Str8 90 ZIP_ra (30.8%)、222 b} | - +变体。id {Int32 1348 ZIP_ra (35.7%), 1.9 k} | - +位置{Int32 1348 ZIP_ra (86.4%), 4.6 k} | - +染色体{Str8 1348 ZIP_ra (2.66%)、91 b} | - +等位基因{Str8 1348 ZIP_ra (17.2%)、928 b} | - +基因型[]* | \ +数据{Bit2 2 x90x1348 ZIP_ra (28.4%), 16.8 k} * | +阶段[]| \——+数据{Bit1 90 x1348 ZIP_ra (0.36%)、55 b} * | - +注释[]| | - + id {Str8 1348 ZIP_ra (41.0%), 5.8 k} | | - + 8 {Float32 1348 ZIP_ra (0.91%), 49 b} | | - +过滤{Int32,因素1348 ZIP_ra (0.89%)、48 b} * | | - +信息[]| | | - + AA {Str8 1348 ZIP_ra (24.2%)、653 b} * | | \——+款HM2 {Bit1 1348 ZIP_ra (117.2%)、198 b} * | \[] | - +格式\——+ DP[] * | \ +数据{Int32 90 x1348 ZIP_ra (33.8%), 160.3 k} \——+样品。注释[]\——+家庭{Str8 90 ZIP_ra (34.7%)、135 b}

方法——关键功能

表1:SeqArray包中的关键功能。

函数	描述
seqVCF2GDS	重新格式化VCF文件
seqSetFilter	定义了一个样品或变量的数据子集
seqGetData	得到数据从SeqArray文件定义过滤器
seqApply	适用于一个用户定义的函数对数组的利润率
seqParallel	应用函数在计算集群

基准

数据集
- 1000人基因工程三期,染色体1
- 6468094个变异,2504人
- 原来的VCF。gz文件:1.2 g
- 重新格式化SeqArray文件:458 (zlib压缩)
计算参考等位基因的频率
1. R代码(顺序版)
2. R代码(并行版本)
3. 通过Rcpp R和C + +集成方案

基准测试1(按顺序)

#加载R包图书馆(SeqArray)#打开文件genofile < -seqOpen(“1 kg_chr1.gds”)#应用一个用户定义函数/变量system.time(afreq < -seqApply(genofile“基因”,有趣的=函数(x) {的意思是(x= =0 l,na.rm =真正的)},人群收税=“替身”,利润=“by.variant”))

10.8分钟在Linux上,英特尔至强处理器@2GHz和128 gb RAM函数(x){意味着(x = = 0 l, na.rm = TRUE)}是一个用户定义的函数,在哪里x是一个整数矩阵:

样本等位基因(,1][,2][,3][,4][,5](1,)010NA1(2,)00010

0 -参考等位基因,1 -第一个替代等位基因

基准测试2(并行)

seqParallel ()基因型分裂成4重叠的部分根据不同的内核。

#加载R包图书馆(平行)#创建一个计算集群4核seqParallelSetup(4)#并行运行system.time(afreq < -seqParallel(gdsfile =genofile,有趣的=函数(f) {seqApply(f,“基因”,人群收税=“替身”,利润=“by.variant”,有趣的=函数(x)的意思是(x= =0 l,na.rm =真正的))},分=“by.variant”))

3.1分钟(与10.8年的测试(1)

基准测试3 (c++集成)

图书馆(Rcpp)#动态定义内联C / c++函数RcppFunction(“双RefAlleleFreq (IntegerMatrix x) {int nrow = x.nrow (), ncol = x.ncol ();int问= 0,zero_cnt = 0, g;for (int i = 0;我< nrow;我+ +){for (int j = 0;j < ncol;j + +) {如果((g = x (i, j)) ! = NA_INTEGER) {问+ +;如果(g = = 0) zero_cnt + +;}}}返回的两倍(zero_cnt) /问;}')system.time(afreq < -seqApply(genofile“基因”RefAlleleFreq,人群收税=“替身”,利润=“by.variant”))

1.5分钟(快很多!与10.8年的测试(1)

结论

SeqArray是极大的兴趣

R用户参与数据分析大规模的序列变异
尤其是并行/高性能计算的经验有限

SeqVarTools (Bioconductor)

变异分析,类似等位基因频率,HWE,孟德尔错误,等等
函数显示基因型/注释以一种可读的格式

SNPRelate (Bioconductor)

并行计算工具集亲缘和主成分分析

资源

https://gds-stat.s3.amazonaws.com/download/1000g/index.html

1000人基因工程阶段3:

常染色体(2.60 gb, 2504人和81271745个变异):1 kg_all.autosome.phase3_shapeit2_mvncall_integrated_v5a.20130502.genotypes.gds
染色体X (94.1 mb, 2504个人和3468093个变异):1 kg_all.chrx.phase3_shapeit2_mvncall_integrated_v1b.20130502.genotypes.gds
染色体Y (2.70 mb, 1233男性和62042个变异):1 kg_all.chry.phase3_integrated_v2a.20130502.genotypes.gds

确认

部门——西雅图华盛顿大学生物统计学

遗传分析中心:

斯蒂芬妮·m·Gogarten
大卫·莱文
凯茜劳里