基因测序文件解析及Doris数据存储

1. 基因测序与分析

1.1 基因测序技术

基因测序技术也称作DNA测序技术，即获得目的DNA片段碱基排列顺序的技术，获得目的DNA片段的序列是进一步进行分子生物学研究和基因改造的基础。

1977年，Walter Gilbert和Frederick Sanger发明了第一台测序仪，并应用其测定了第一个基因组序列，噬菌体X174，全长5375个碱基。由此开始，人类获得了探索生命遗传本质的能力，生命科学的研究进入了基因组学的时代，到至今为止的四十年时间内，测序技术已取得了相当大的发展，从第一代发展到了第三代测序技术。Sanger所发明的测序方法被称为第一代测序技术，该技术直到现在依然被广泛使用，但是其一次只能获得一条长度在700～1000个碱基的序列，无法满足现代科学发展对生物基因序列获取的迫切需求。

高通量测序 (High-Throughput Sequencing, HTS) 是对传统Sanger测序的革命性变革，其解决了一代测序一次只能测定一条序列的限制，一次运行即可同时得到几十万到几百万条核酸分子的序列，因此也被称为第二代测序技术。

第二代测序技术虽然测序的通量大大增加，但是其获得单条序列长度很短，想要得到准确的基因序列信息依赖于较高的测序覆盖度和准确的序列拼接技术，因此最终得到的结果中会存在一定的错误信息。因此，科研人员又发明了第三代测序技术也称为单分子测序技术，该技术在保证测序通量的基础上，对单条长序列进行从头测序，能够直接得到长度在数万个碱基的核酸序列信息。

1.2 全基因组关联分析

GWAS(Genome-wide association study)，即全基因组关联分析，是指在全基因组范围内找出存在的序列变异，即单核苷酸多态性(SNP)，从中筛选出与疾病相关的SNPs。全基因组关联研究是一种检测特定物种中不同个体间的全部或大部分基因，从而了解不同个体间的基因变化有多大的一种方法。不同的变化带来不同的性状，如各种疾病的不同。

1.3 基因分析工具

1.3.1 Plink工具

Plink是一个免费的开源全基因组关联分析工具集，旨在以高效率计算的方式执行一系列基本的、大规模的分析。Plink的重点分析对象是基因型或者表型数据，可为后续的可视化、注释和结果存储提供一些支持。

• Plink官网：http://zzz.bwh.harvard.edu/plink/ (opens new window)

• Plink功能：数据管理、质控、人群分层检测、关联分析-SNP、多marker预测、haplotype分析、拷贝数变异分析、异位显性分析、关联分析-基因、基因环境互作分析、meta分析等。

1.3.2 Admixture工具

Admixture是一款快速分析群体遗传结构的工具，一般进行群体分析所使用的软件为Structure，但Structure的运行速度较慢，如今Admixture凭借其高速的运算速度逐渐成为群体遗传结构分析的主流软件。

• Admixture官网：https://bioinformaticshome.com/tools/descriptions/ADMIXTURE.html (opens new window)
• 群体遗传结构（Structure）指遗传变异在物种或群体中的一种非随机分布。按照地理分布或其他标准可将一个群体分为若干亚群，处于同一亚群内的不同个体亲缘关系较高，而亚群与亚群之间则亲缘关系稍远。群体结构分析有助于理解进化过程，并且可以通过基因型和表型的关联研究确定个体所属的亚群。

2. 基因测序VCF文件处理

2.1 解析基因变异记录VCF文件

2.2.1 文件格式说明

VCF是用于描述SNP（单个碱基上的变异），INDEL（插入缺失标记）和SV（结构变异位点）结果的生物基因数据文件，通过基因测序技术获得。

文件的开头是一堆以“##”开始的注释行，包含了文件的基本信息。然后是以“#”开头的一行，共9+n个部分，前九部分标注的是后面行每部分代表的信息，相当于表头。后面部分是样本名称，可以有多个。注释行结束后是具体的突变信息，每一行分为9+n个部分，每部分之间用制表符（‘\t’）分隔。

通常会对样本压缩成.vcf.gz格式，这个我们下面使用PyVCF库可以直接解析，无需解压。如果需要解压取样的话，可以使用如下命令：

$ gunzip test.vcf.gz                             // 解压.vcf.gz得到.vcf文件，并删除原始.vcf.gz文件
$ cat test.vcf | head -n 10000 > test_1w.vcf     // 对.vcf文件进行取样

注：使用gunzip解压成功后，会自动把原始的.vcf.gz文件删掉，如果需要保留，可以提前拷贝一份。

2.2.2 使用PyVCF库解析文件

PyVCF：通过正则表达式把vcf文件信息转换成结构化的信息，简化了vcf文件的处理过程，方便后续提取相关参数及处理。

• 项目地址：https://github.com/jamescasbon/PyVCF (opens new window)

• 安装依赖：直接安装PyVCF可能会遇到error in setup command: use_2to3 is invalid.报错，因为在setuptools 58之后的版本已经废弃了use_2to3。

  $ pip3 install setuptools==57.5.0 
  $ pip3 install PyVCF==0.6.8
  

PyVCF库使用介绍

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')
for record in vcf_reader:
    print(record)
    
>>> 输出内容   
Record(CHROM=NC_006583.3, POS=548048, REF=G, ALT=[<NON_REF>])
Record(CHROM=NC_006583.3, POS=548057, REF=T, ALT=[<NON_REF>])
Record(CHROM=NC_006583.3, POS=548071, REF=C, ALT=[<NON_REF>])
Record(CHROM=NC_006583.3, POS=548075, REF=G, ALT=[<NON_REF>])
Record(CHROM=NC_006583.3, POS=548076, REF=T, ALT=[<NON_REF>])

调用vcf.Reader类处理vcf文件，vcf文件信息就被保存到vcf_reader中了。它是一个可迭代对象，它的迭代元素都是一个_Record对象的实例，保存着非注释行的一行信息，即变异位点的具体信息。通过它，我们可以很轻易地得到位点的详细信息。

vcf.gz文件可以指定读哪个染色体的哪个位置，但需要目录处有对应的.vcf.gz.tbi文件并且安装pysam库。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')

vcf_reader.fetch('chr_1', 2, 5)
for record in vcf_reader:
    print(record)

详见：https://pyvcf.readthedocs.io/en/latest/API.html#vcf-reader (opens new window)

[1] _Record对象-位点信息的储存形式

class vcf.model._Record(CHROM, POS, ID, REF, ALT, QUAL, FILTER, INFO, FORMAT, sample_indexes, samples=None)

下面分别解释一下 CHROM、POS、ID、REF、ALT、QUAL、FILTRE、INFO、FORMAT 等9部分的含义。

• CHROM：染色体名称，类型为str。
• POS：位点在染色体上的位置，类型为int。
• ID：一般是突变的rs号，类型为str。如果是‘.’，则为None。
• REF：参考基因组在该位点上的碱基，类型为str。
• ALT：在该位点的测序结果，类型为list，代表了突变的几种类型。
• QUAL：该位点的测序质量，类型为int或float。
• FILTER：过滤信息。将FILTER列按分号分隔形成的字符串列表，类型为list。如果未给出参数则为None。
• INFO：该位点的一些测试指标。将‘=’前的参数作为键，后面的参数作为值，构建成的字典，类型为dict。
• FORMAT：基因型信息。保存vcf的FORMAT列的原始形式，类型为str。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')
for record in vcf_reader:
    print(record)
    print("=====")
    print(type(record.CHROM), record.CHROM)
    print(type(record.POS), record.POS)
    print(type(record.ID), record.ID)
    print(type(record.REF), record.REF)
    print(type(record.ALT), record.ALT)
    print(type(record.QUAL), record.QUAL)
    print(type(record.FILTER), record.FILTER)
    print(type(record.INFO), record.INFO)
    print(type(record.FORMAT), record.FORMAT)
    break
    
>>> 输出内容
Record(CHROM=NC_006583.3, POS=1, REF=N, ALT=[<NON_REF>])
=====
<class 'str'> NC_006583.3
<class 'int'> 1467
<class 'NoneType'> None
<class 'str'> G
<class 'list'> [<NON_REF>]
<class 'NoneType'> None
<class 'NoneType'> None
<class 'dict'> {'END': 1471}
<class 'str'> GT:DP:GQ:MIN_DP:PL

除了这些基本属性之外，_Record对象还有一些其他属性：

• samples：把FORMAT信息作为键，后面对应的信息做为值，构建成的字典（CallData对象），以及sample名称，这两个值组成一个Call对象，共同构成samples的一个元素。这样就把sample和基因型信息给关联起来，按下标访问每一个Call对象。samples类型为list。
• start：突变开始的位置。
• end：突变结束的位置。
• alleles：该位点所有的可能情况，由REF和ALT参数组成的列表（REF类型是str，ALT参数是_AltRecord对象的子类实例），类型是list。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')
for record in vcf_reader:
    # 按下标访问Call，按.sample访问sample，按键访问FORMAT对应信息
    print(record.samples, '\n', record.samples[0].sample, '\n', record.samples[0]['GT'])  
    print(record.start, record.POS, record.end)
    print(record.REF, record.ALT, record.alleles)
    break
    
>>> 输出内容
[Call(sample=A1, CallData(GT=0/0, DP=0, GQ=0, MIN_DP=0, PL=[0, 0, 0]))] 
A1 
0/0
0 1 1
N [<NON_REF>] ['N', <NON_REF>]

_Record对象方法：

• 对象之间比较大小方法：根据染色体名称和位置信息比较。
• 迭代方法：对samples里的元素进行迭代。
• 字符串方法：只返回CHROM，POS，REF，ALT四列信息。
• genotype(name)方法，和samples按下标访问不同，这个方法提供按sample名称进行访问的功能。
• add_format(fmt)， add_filter(flt)， add_info(info, value=True)：给相应的属性添加元素。
• get_hom_refs()：拿到samples中该位点未突变的所有sample，返回列表。
• get_hom_alts()：拿到samples中该位点100%突变的所有sample，返回列表。
• get_hets()：拿到samples中该位点基因型为杂合的所有sample，返回列表。
• get_unknown()：拿到samples中该位点基因型未知的所有sample，返回列表。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')

# 按染色体名称和位置进行比较
record = next(vcf_reader)
record2 = next(vcf_reader)
print(record < record2)

# 按samples列表进行迭代
for i in record:
    print(i)
    
# 按sample名字进行访问
print(record.genotype('A1'))

>>> 输出内容
True
Call(sample=A1, CallData(GT=0/0, DP=0, GQ=0, MIN_DP=0, PL=[0, 0, 0]))
Call(sample=A1, CallData(GT=0/0, DP=0, GQ=0, MIN_DP=0, PL=[0, 0, 0]))

_Record对象还有很多有用的方法属性：

• num_called：该位点已识别的sample数目。
• call_rate：已识别的sample数目占sample总数的比例。
• num_hom_ref，num_hom_alt，num_het，num_unknown：四种基因型的数量
• aaf：所有sample等位基因的频率（即除开REF），返回列表。
• heterozygosity：该位点的杂合度，0.5为杂合突变，0为纯合突变。
• var_type：突变类型，包括‘snp’，‘indel’，‘sv’（structural variant），‘unknown’。
• var_subtype：更加细化的突变类型，如‘indel’包括‘del’，‘ins’，‘unknown’。
• is_snp，is_indel，is_sv，is_transition，is_deletion：判断突变是不是snp，indel，sv，transition，deletion等。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')

record = next(vcf_reader)
print(record)
print(record.samples)
print(record.num_called)
print(record.call_rate)
print(record.num_hom_ref)
print(record.aaf)
print(record.num_het)
print(record.heterozygosity)
print(record.var_type)
print(record.var_subtype)
print(record.is_snp)
print(record.is_indel)

>>> 输出内容
[Call(sample=A1, CallData(GT=0/0, DP=0, GQ=0, MIN_DP=0, PL=[0, 0, 0]))]
1
1.0
1
[0.0]
0
0.0
unknown
unknown
False
False

[2] Reader对象-处理vcf文件，构建结构化信息

class Reader(fsock=None, filename=None, compressed=None, prepend_chr=False, strict_whitespace=False, encoding='ascii')

在读取vcf文件时，总共有六个参数可供选择。

• fsock：目标文件的文件对象，可以用open(文件名)得到这个文件对象。
• filename：文件名，当fsock和filename同时存在时，优先考虑fsock。
• compressed：是否要解压，不提供参数时由程序自行判断（以文件名是否以.gz结尾判断是否需要解压）。
• prepend_chr：在保存染色体名称时，是否加前缀‘chr’，默认不加，如果vcf文件的染色体名称本来没有前缀‘chr’，可设置为True，自动加上。
• strict_whitespace：是否严格以制表符‘\t’分隔数据。True则表示严格按制表符分，False表示可以夹杂空格。
• encoding：文件编码。

头文件信息主要保存在Reader对象的属性中，包括alts，contigs，filters，formats，infos，metadata。

import vcf

vcf_reader = vcf.Reader(filename=r'./data/A1.raw.g.vcf.gz')
print(vcf_reader.alts)
print(vcf_reader.alts['NON_REF'].id)
print(vcf_reader.alts['NON_REF'].desc)

>>> 输出内容
OrderedDict([('NON_REF', Alt(id='NON_REF', desc='Represents any possible alternative allele at this location'))])
NON_REF
Represents any possible alternative allele at this location

2.2 处理VCF文件并生成索引

2.2.1 安装htslib及bcftools工具

htslib是bcftools使用的核心库，单独安装bcftools会缺失很多环境，建议先安装htslib，然后再安装bcftools。

以下在Debian 11 x86_64系统安装htslib是bcftools工具（不同系统的安装方式不同），基础编译环境如下：

$ sudo apt-get update  
$ sudo apt-get install autoconf
$ sudo apt-get install -y make

[1] htslib工具

HTSlib 是一个统一的 C 库的实现，用于访问SAM、CRAM 和 VCF等常见文件格式，用于高通量测序数据，是samtools和bcftools使用的核心库。

• 项目地址：https://github.com/samtools/htslib (opens new window)
• 安装说明：https://github.com/samtools/htslib/blob/develop/INSTALL (opens new window)

$ git clone http://github.com/samtools/htslib.git
$ cd htslib
$ sudo apt-get install autoconf automake make gcc perl zlib1g-dev libbz2-dev liblzma-dev libcurl4-gnutls-dev libssl-dev
$ autoreconf -i
$ git submodule update --init --recursive
$ ./configure    
$ make
$ make install

参考：https://github.com/samtools/htslib (opens new window)、https://github.com/samtools/htslib/blob/develop/INSTALL (opens new window)

[2] bcftools工具

bcftools是一组用于变异获取和操作VCF和BCF的实用工具集合。参数众多，功能也十分强大。

• 项目地址：https://github.com/samtools/bcftools (opens new window)
• 官方文档：http://samtools.github.io/bcftools/bcftools.html (opens new window)

$ git clone http://github.com/samtools/bcftools.git
$ cd bcftools
$ autoheader && autoconf && ./configure --enable-libgsl --enable-perl-filters
$ make
$ make install

2.2.2 生成索引随机读取VCF位点

建立索引可以加快数据读取速度，而且建立索引后可以实现随机读取VCF位点。bcftools工具index命令用于对VCF文件建立索引，要求输入的VCF文件必须是使用bgzip压缩之后的文件，支持.csi和.tbi两种索引，默认情况下建立的索引是.csi格式。

[1] 将vcf文件压缩转换为vcf.gz

$ bcftools view your_file.vcf -Oz -o your_file.vcf.gz

[2] 为vcf.gz构建索引

$ bcftools index your_file.vcf.gz

注：必须要先压缩成gz格式，直接用vcf否则会报错index: the file is not BGZF compressed, cannot index: your_file.vcf。在我的实验环境下，123s对一个具有23亿（由于测序VCF里只有突变位置，所以实际远小于23亿）基因位点的生物样本建立了索引。

[3] 读取任意范围的位点信息

$ tabix your_file.vcf.gz chr_1:1-100000000

注：chr_1是染色体名，1-100000000是位点范围。随机读取的性能挺强的，在我的实验环境下，13秒取出了3600w基因位点信息。

3. 数据存储技术选型

3.1 数据存储选型分析

3.1.1 数据特点分析

需求情景：每个基因测序样本的VCF文件解析完之后大约是23亿数据，目前有100个样本的数据需要存储，数据规模达到千亿级别。

数据源（VCF文件）经过各种数据集成和加工处理后，需要入库到实时数仓系统中以备将来查询使用，VCF的基因数据有以下几个特点：

• 不需要更新，测序数据生成后，因为不再变化 ，所以不需要对其进行更新操作；
• 数据量大，单个VCF文件的数据条数已经是十亿级别，一百个VCF文件就是千亿级别；
• 单条数据很短，且各字段取值范围小，和其他数据比较相似，字符串短且无语义，对其按列进行压缩会有很好的压缩比；
• 存储的后期压力大。虽然目前只有一百个VCF文件，但后期每进行一次测序都会形成一个VCF文件，随着业务的增长，数据量会爆炸性增长。
• 数据使用较不频繁，使用时对吞吐量的要求较高。

3.1.2 数据存储选型理由

针对基因数据的特点，结合硬件的现实情况和技术的可行性，我们选型了Apache Doris作为主数仓，技术上有以下几点针对性的考虑：

• 基于数据量大和不再更新的特性，传统的DBMS的数据体量太小，对不再更新的数据没有针对性的优化，无法胜任基因数据的存储任务，因而倾向于以OLAP类产品作为主数据仓库，这类产品的内部数据会按某种规则排序存放，更新性能不好，导入和查询性能则十分优秀，Doris就是近期出现的性能最好的OLAP产品，Doris的收费衍生版本SelectDB在分析型数据库性能测试排行榜ClickBench的测试中表现亮眼，领先幅度明显。
• 由于数据尺寸较为庞大，且项目前期的硬件资源较为紧张，其中也包括存储硬件资源的紧张，再加上未来对数据爆炸性增长的预期，存储资源的占用是一个很大的问题，因而势必要对数据进行压缩，列存储是不二选择，而Doris就是使用了列存储引擎。通过我们的测试，Doris列存储引擎的压缩表现也非常不错，一百个样本的数据在保留一个备份副本的情况下，仍然可以比VCF文件节省至少60%的磁盘存储空间，且数据量越大，压缩效果越好。
• 由于数据爆炸性增长的可能性大，数仓系统的可扩展性就十分重要，其需要在硬件资源告急的情况下，可以平滑扩展数仓系统的集群，而可扩展性正是Doris的主打特性，不仅操作简单，还可以根据并发量、吞吐量、磁盘空间短缺的具体情况采取针对性的扩展措施，避免硬件资源的浪费，这一点我们也已经通过前期的测试予以确认。另外，Doris的可扩展上限非常高，可驾驭数百台服务器的集群。
• 自带一些数据导入和数据流动的工具，自带高可用，兼容MySQL语法，在使用上起手容易；优化措施较为丰富，支持多种索引结构，支持物化视图，优化上限较高但是有一些学习成本；社区状况良好，后续支持有力；较为独立，不依赖于第三方系统。这些有利因素便于项目的前期推动和后期的技术保障。

3.2 Doris数据存储

3.2.1 Doris介绍

Apache Doris 是一个基于 MPP 架构的高性能、实时的分析型数据库，以极速易用的特点被人们所熟知，仅需亚秒级响应时间即可返回海量数据下的查询结果，不仅可以支持高并发的点查询场景，也能支持高吞吐的复杂分析场景。基于此，Apache Doris 能够较好的满足报表分析、即席查询、统一数仓构建、数据湖联邦查询加速等使用场景，用户可以在此之上构建用户行为分析、AB 实验平台、日志检索分析、用户画像分析、订单分析等应用。

Apache Doris 最早是诞生于百度广告报表业务的 Palo 项目，2017 年正式对外开源，2018 年 7 月由百度捐赠给 Apache 基金会进行孵化，之后在 Apache 导师的指导下由孵化器项目管理委员会成员进行孵化和运营。目前 Apache Doris 社区已经聚集了来自不同行业近百家企业的 300 余位贡献者，并且每月活跃贡献者人数也接近 100 位。 2022 年 6 月，Apache Doris 成功从 Apache 孵化器毕业，正式成为 Apache 顶级项目（Top-Level Project，TLP）

Apache Doris 如今在中国乃至全球范围内都拥有着广泛的用户群体，截止目前， Apache Doris 已经在全球超过 500 家企业的生产环境中得到应用，在中国市值或估值排行前 50 的互联网公司中，有超过 80% 长期使用 Apache Doris，包括百度、美团、小米、京东、字节跳动、腾讯、网易、快手、微博、贝壳等。同时在一些传统行业如金融、能源、制造、电信等领域也有着丰富的应用。

项目地址：https://github.com/apache/doris (opens new window) 官方文档：https://doris.apache.org/zh-CN/ (opens new window)（官方文档非常详细）

3.2.2 Doris使用场景

数据源经过各种数据集成和加工处理后，通常会入库到实时数仓 Doris 和离线湖仓（Hive, Iceberg, Hudi 中），Apache Doris 被广泛应用在以下场景中。

1）报表分析

• 实时看板 （Dashboards）
• 面向企业内部分析师和管理者的报表
• 面向用户的高并发报表分析（Customer Facing Analytics）。比如面向网站主的站点分析、面向广告主的广告报表，并发通常要求成千上万的 QPS ，查询延时要求毫秒级响应。著名的电商公司京东在广告报表中使用 Apache Doris ，每天写入 100 亿行数据，查询并发 QPS 上万，99 分位的查询延时 150ms。

2）即席查询（Ad-hoc Query）：面向分析师的自助分析，查询模式不固定，要求较高的吞吐。小米公司基于 Doris 构建了增长分析平台（Growing Analytics，GA），利用用户行为数据对业务进行增长分析，平均查询延时 10s，95 分位的查询延时 30s 以内，每天的 SQL 查询量为数万条。

3）统一数仓构建 ：一个平台满足统一的数据仓库建设需求，简化繁琐的大数据软件栈。海底捞基于 Doris 构建的统一数仓，替换了原来由 Spark、Hive、Kudu、Hbase、Phoenix 组成的旧架构，架构大大简化。

4）数据湖联邦查询：通过外表的方式联邦分析位于 Hive、Iceberg、Hudi 中的数据，在避免数据拷贝的前提下，查询性能大幅提升。

3.2.3 Doris技术架构

Doris整体架构如下图所示，Doris 架构非常简单，只有两类进程：

• Frontend（FE），主要负责用户请求的接入、查询解析规划、元数据的管理、节点管理相关工作。
• Backend（BE），主要负责数据存储、查询计划的执行。

这两类进程都是可以横向扩展的，单集群可以支持到数百台机器，数十 PB 的存储容量。并且这两类进程通过一致性协议来保证服务的高可用和数据的高可靠。这种高度集成的架构设计极大的降低了一款分布式系统的运维成本。

在使用接口方面，Doris 采用 MySQL 协议，高度兼容 MySQL 语法，支持标准 SQL，用户可以通过各类客户端工具来访问 Doris，并支持与 BI 工具的无缝对接。

在存储引擎方面，Doris 采用列式存储，按列进行数据的编码压缩和读取，能够实现极高的压缩比，同时减少大量非相关数据的扫描，从而更加有效利用 IO 和 CPU 资源。

Doris 也支持比较丰富的索引结构，来减少数据的扫描：

• Sorted Compound Key Index，可以最多指定三个列组成复合排序键，通过该索引，能够有效进行数据裁剪，从而能够更好支持高并发的报表场景
• Z-order Index ：使用 Z-order 索引，可以高效对数据模型中的任意字段组合进行范围查询
• Min/Max ：有效过滤数值类型的等值和范围查询
• Bloom Filter ：对高基数列的等值过滤裁剪非常有效
• Invert Index ：能够对任意字段实现快速检索

在存储模型方面，Doris 支持多种存储模型，针对不同的场景做了针对性的优化：

• Aggregate Key 模型：相同 Key 的 Value 列合并，通过提前聚合大幅提升性能
• Unique Key 模型：Key 唯一，相同 Key 的数据覆盖，实现行级别数据更新
• Duplicate Key 模型：明细数据模型，满足事实表的明细存储

Doris 也支持强一致的物化视图，物化视图的更新和选择都在系统内自动进行，不需要用户手动选择，从而大幅减少了物化视图维护的代价。

在查询引擎方面，Doris 采用 MPP 的模型，节点间和节点内都并行执行，也支持多个大表的分布式 Shuffle Join，从而能够更好应对复杂查询。

Doris 查询引擎是向量化的查询引擎，所有的内存结构能够按照列式布局，能够达到大幅减少虚函数调用、提升 Cache 命中率，高效利用 SIMD 指令的效果。在宽表聚合场景下性能是非向量化引擎的 5-10 倍。

Doris 采用了 Adaptive Query Execution 技术，可以根据 Runtime Statistics 来动态调整执行计划，比如通过 Runtime Filter 技术能够在运行时生成生成 Filter 推到 Probe 侧，并且能够将 Filter 自动穿透到 Probe 侧最底层的 Scan 节点，从而大幅减少 Probe 的数据量，加速 Join 性能。Doris 的 Runtime Filter 支持 In/Min/Max/Bloom Filter。

在优化器方面 Doris 使用 CBO 和 RBO 结合的优化策略，RBO 支持常量折叠、子查询改写、谓词下推等，CBO 支持 Join Reorder。目前 CBO 还在持续优化中，主要集中在更加精准的统计信息收集和推导，更加精准的代价模型预估等方面。

4. 搭建Doris分布式集群

准备了4台满足上述部署基本要求的服务器，所使用的系统为 CentOS Linux release 7.9.2009 (Core)

4.1 Doris部署基本要求

至少需要准备四台服务器，一台运行FE，三台运行BE，服务器有如下的要求：

• 操作系统版本要求Centos7.1及以上或者Ubuntu16.04及以上。

• JDK版本要求1.8及以上。

• 系统最大打开文件句柄数修改为65536。

  vim /etc/security/limits.conf
  * soft nofile 65536
  * hard nofile 65536
  

• 时钟同步，要求各服务器之间的时间差不能超过5秒。

• 要求禁用各服务器的交换分区。

• 要求各服务器的文件系统为ext4。

• 生产环境要求各服务器有不低于16核64GB内存200GB SSD的配置。

• 下面的表格列出了各服务器的通讯端口要求：

  实例名称	端口名称	默认端口	通讯方向	说明
  BE	be_port	9060	FE --> BE	BE 上 thrift server 的端口，用于接收来自 FE 的请求
  BE	webserver_port	8040	BE <--> BE	BE 上的 http server 的端口
  BE	heartbeat_service_port	9050	FE --> BE	BE 上心跳服务端口（thrift），用于接收来自 FE 的心跳
  BE	brpc_port	8060	FE <--> BE, BE <--> BE	BE 上的 brpc 端口，用于 BE 之间通讯
  FE	http_port	8030	FE <--> FE，用户 <--> FE	FE 上的 http server 端口
  FE	rpc_port	9020	BE --> FE, FE <--> FE	FE 上的 thrift server 端口，每个fe的配置需要保持一致
  FE	query_port	9030	用户 <--> FE	FE 上的 mysql server 端口
  FE	edit_log_port	9010	FE <--> FE	FE 上的 bdbje 之间通信用的端口
  Broker	broker_ipc_port	8000	FE --> Broker, BE --> Broker	Broker 上的 thrift server，用于接收请求

• 官方推荐使用带avx2的CPU，对性能有提升。

详见官方部署文档：https://doris.apache.org/zh-CN/docs/dev/install/install-deploy/ (opens new window)

4.2 安装部署Doris集群

4.2.1 安装文件下载

官方下载地址：https://doris.apache.org/zh-CN/download/ (opens new window)

以1.1.4版本为例，先使用 cat /proc/cpuinfo 命令查看 CPU 是否支持 avx2 指令集，我这里是支持的。

FE下载地址

• x64(avx2)：https://mirrors.tuna.tsinghua.edu.cn/apache/doris/1.1/1.1.4-rc01/apache-doris-fe-1.1.4-bin.tar.gz (opens new window)

BE下载地址

• 如果CPU支持avx2，需要选择avx2版本的BE：https://mirrors.tuna.tsinghua.edu.cn/apache/doris/1.1/1.1.4-rc01/apache-doris-be-1.1.4-bin-x86_64.tar.gz (opens new window)
• 如果CPU不支持avx2，需要选择noavx2版本的BE：https://mirrors.tuna.tsinghua.edu.cn/apache/doris/1.1/1.1.4-rc01/apache-doris-be-1.1.4-bin-x86_64-noavx2.tar.gz (opens new window)

4.2.2 安装部署过程

Step1：使用tar命令解压Doris的FE和BE的文件，并确保解压好的文件位于SSD提供的存储空间内。

Step2：修改配置文件，修改priority_networks项的配置为24位掩码的真实IP地址。FE和BE都要配置，配置文件位于./conf目录下。

Step3：先启动FE，再启动BE，启动脚本均位于./bin目录下。

$ ./bin/start_fe.sh --daemon    // 启动FE
$ ./bin/start_be.sh --daemon    // 启动BE

Step4：FE可通过访问管理页面判断是否启动成功，地址为http://FE_IP:8030，默认用户名为root，密码为空。

BE可查看其日志或在FE的管理页面上判断其是否启动成功，建议第一次启动BE时，先不用带--daemon参数使其前台运行，通过日志查看启动是否有问题，确认启动没问题之后再改为后台运行。

4.3 连接测试Doris集群

4.3.1 连接Doris集群

可以使用Navicat连接和操作Doris集群，数据库类型选择MySQL即可，提供FE的IP地址和端口（默认9030)，默认用户名为root，默认密码为空。

连接成功后需要将所有BE全部添加到FE，添加BE使用SQL语句，示例如下：

ALTER SYSTEM ADD BACKEND "BE_IP:9050";

4.3.2 测试Doris集群

添加全部BE后可新建测试数据库。

create database demo;

可执行官方提供的测试建表语句：

use demo;

CREATE TABLE IF NOT EXISTS demo.example_tbl
(
    `user_id` LARGEINT NOT NULL COMMENT "用户id",
    `date` DATE NOT NULL COMMENT "数据灌入日期时间",
    `city` VARCHAR(20) COMMENT "用户所在城市",
    `age` SMALLINT COMMENT "用户年龄",
    `sex` TINYINT COMMENT "用户性别",
    `last_visit_date` DATETIME REPLACE DEFAULT "1970-01-01 00:00:00" COMMENT "用户最后一次访问时间",
    `cost` BIGINT SUM DEFAULT "0" COMMENT "用户总消费",
    `max_dwell_time` INT MAX DEFAULT "0" COMMENT "用户最大停留时间",
    `min_dwell_time` INT MIN DEFAULT "99999" COMMENT "用户最小停留时间"
)
AGGREGATE KEY(`user_id`, `date`, `city`, `age`, `sex`)
DISTRIBUTED BY HASH(`user_id`) BUCKETS 1
PROPERTIES (
    "replication_allocation" = "tag.location.default: 1"
);

注：如果不能正常建表，并提示1105错误，有可能是前面的要求或步骤没有正确完成，例如部分端口没有正确开放等，可查看FE和BE的日志以确定错误原因。

4.4 修改Doris集群配置

本节都是不必要项，自身有默认值，如果需要修改则根据实际业务需求进行配置即可。

4.4.1 查看与修改内存

$ SHOW VARIABLES LIKE "%mem_limit%";          // 查看内存限制，默认为2GB
$ set global exec_mem_limit = 8589934592;     // 修改内存限制为8GB

4.4.2 查看与修改超时时间

修改超时时间：Doris修改超时时间-官方文档 (opens new window)

$ SHOW VARIABLES LIKE "%query_timeout%";      // 查看超时时间，默认为300s
$ set global query_timeout = 6000;            // 修改超时时间为6000s

5. Doris集群的基本使用

Doris由于兼容MySQL协议，基本使用与MySQL差不多，可以使用SQL语句进行查询，但也有一些地方不同。

• 虽然查询操作与SQL语句类似，但操作数据表的语句并不兼容，一些内容也不支持修改，因此使用Navicat进行操作的时候会看到一些报错。
• Doris的数据模型主要分为3类：Aggregate、Unique、Duplicate，详见：数据模型-Doris官方文档 (opens new window)

5.1 数据导入Doris集群

5.1.1 使用Insert方式导入数据

Insert into语句的使用方式和 MySQL 等数据库中的 Insert Into 语句使用方式类似。但在 Doris 中，所有的数据写入都是一个独立的导入作业。

主要的 Insert Into 命令包含以下两种：

• INSERT INTO tbl SELECT ...
• INSERT INTO tbl (col1, col2, ...) VALUES (1, 2, ...), (1,3, ...);

批量插入的效率要远高于单条插入，使用 Insert into 的第2种方式，仅用于简单测试或低频少量的操作，不要在生产环境中使用。见：插入-Doris官方文档 (opens new window)

INSERT INTO example_tbl VALUES
(1000, "baidu1", 3.25),
(2000, "baidu2", 4.25),
(3000, "baidu3", 5.25);

5.1.2 使用订阅Kafka方式导入数据

使用订阅kafka消息的方式，将数据存入Doris。

// 查询当前有哪些任务
show routine Load; 

// 停止指定任务
STOP ROUTINE LOAD FOR your_db.your_job; 

// 暂停指定任务
PAUSE ROUTINE LOAD FOR your_db.your_job;

// 暂停所有任务
PAUSE ALL ROUTINE LOAD;

// 重启指定任务
RESUME ROUTINE LOAD FOR your_db.your_job;

// 重启所有任务
RESUME ALL ROUTINE LOAD;

// 创建任务
CREATE ROUTINE LOAD your_db.your_job ON your_table COLUMNS (
	FIELD1,
	FIELD2,
	FIELD3
) PROPERTIES (
	"max_batch_interval" = "10",
	"max_batch_rows" = "100000000",
	"max_batch_size" = "524288000",
	"strict_mode" = "false",
	"format" = "json"
)
FROM
	KAFKA (
		"kafka_broker_list" = "kafka_ip:port",
		"kafka_topic" = "your_topic",
		"kafka_partitions" = "0",
		"kafka_offsets" = "0"
	)

注意事项：

1）如果删除了数据表，配置的任务也会随之删除。支持对一张表配置多个任务，以实现从多个kafka topic里共同获取数据。

2）任务的配置参数有范围限制，所以单个topic写入doris的任务是有瓶颈的，如果达不到需求，可以从生产端改成写入多个Topic，然后再配置多个任务，写入到一张表里。

5.1.3 使用DataX工具导入数据

DataX (opens new window) 是阿里云 DataWorks数据集成 的开源版本，在阿里巴巴集团内被广泛使用的离线数据同步工具/平台。DataX 实现了包括 MySQL、Oracle、SqlServer、Postgre、HDFS、Hive、ADS、HBase、TableStore(OTS)、MaxCompute(ODPS)、Hologres、DRDS 等各种异构数据源之间高效的数据同步功能。

• DataX 官方文档上虽然未显示支持 Doris，但是 Doris官方文档上有 DataX doriswriter 插件，应该是 Doris 官方自己开发的。它是利用 Doris 的 Stream Load 功能进行数据导入的，需要配合 DataX 服务一起使用，用于通过 DataX 同步其他数据源的数据到 Doris 中。

• DataX doriswriter 插件官方文档：https://doris.apache.org/zh-CN/docs/ecosystem/datax/ (opens new window)

5.2 使用Python操作Doris

使用Python操作Doris时，可以使用操作MySQL的pymysql库来实现。

# -*- coding: utf-8 -*-

import pymysql


# 查询Doris数据库
def query_data_by_doris():
    conn = pymysql.connect(host="xxx.xxx.xxx.xxx", port=9030, user="root", password="", database="demo", charset='utf8')
    mycursor = conn.cursor()
    mycursor.execute("select * from example_table")
    result = mycursor.fetchall()
    for data in result:
        print(data)


if __name__ == '__main__':
    query_data_by_doris()

5.3 Doris数据查询调优

5.3.1 查看Doris集群的数据量

由于Doris里的数据量太过于庞大，使用count(1)的方式查询会超时，Doris官方提供了查询数据量的方法。详见：SHOW-DATA (opens new window)

• 功能：该语句用于展示数据量、副本数量以及统计行数。

• 语法：SHOW DATA [FROM db_name[.table_name]] [ORDER BY ...];

• 实例：SHOW DATA FROM example_db.my_table;

  这里是生产环境真实数据表查询返回的数据，只不过这里为了脱敏，将数据库名和表名改了。
  
  +-----------+-----------+-----------+--------------+-----------------+
  | TableName | IndexName | Size      | ReplicaCount | RowCount        |
  +-----------+-----------+-----------+--------------+-----------------+
  | my_table  | my_table  | 231.382GB | 4            | 233246375487    |
  |           | Total     | 231.382GB | 4            | 233246375487    |
  +-----------+-----------+-----------+--------------+-----------------+
  
  可以看到Doris的数据压缩比是非常恐怖的，2332亿数据仅用了231.382GB存储空间。
  

5.3.2 数据查询耗时分析

查询耗时分析：Doris耗时分析-官方文档 (opens new window)

$ set enable_profile=true; 

然后登录FE的管理界面，选择到QueryProfile子菜单，可以看到该设置之后的SQL语句查询耗时。

解决方案：

• Doris库表优化：通过添加索引、物化视图等方式去优化查询速度。
• 查询方式优化：分批次去查询，将查询结果拼接起来，类似于ES的滚动查询。

5.3.3 添加Bitmap索引

bitmap index即位图索引，是一种快速数据结构，能够加快查询速度。

基本语法：

// 创建Bitmap索引
CREATE INDEX index_name ON table_name (siteid) USING BITMAP;
// 查看Bitmap索引（创建索引的过程是异步的，需要等待一段时间才可以查询到）
SHOW INDEX FROM example_db.table_name;
// 删除Bitmap索引
DROP INDEX 、index_name ON example_db.table_name;

注意事项：

• bitmap 索引仅在单列上创建。
• bitmap 索引能够应用在 Duplicate、Uniq 数据模型的所有列和 Aggregate模型的key列上。
• bitmap 索引支持的数据类型如下：TINYINT、SMALLINT、INT、BIGINT、CHAR、VARCHAR、DATE、DATETIME、LARGEINT、DECIMAL、BOOL。
• bitmap索引仅在 Segment V2 下生效。当创建 index 时，表的存储格式将默认转换为 V2 格式。

5.3.4 添加BloomFilter索引

BloomFilter是由Bloom在1970年提出的一种多哈希函数映射的快速查找算法。通常应用在一些需要快速判断某个元素是否属于集合，但是并不严格要求100%正确的场合，BloomFilter有以下特点：

• 空间效率高的概率型数据结构，用来检查一个元素是否在一个集合中。
• 对于一个元素检测是否存在的调用，BloomFilter会告诉调用者两个结果之一：可能存在或者一定不存在。
• 缺点是存在误判，告诉你可能存在，不一定真实存在。

布隆过滤器实际上是由一个超长的二进制位数组和一系列的哈希函数组成。二进制位数组初始全部为0，当给定一个待查询的元素时，这个元素会被一系列哈希函数计算映射出一系列的值，所有的值在位数组的偏移量处置为1。

下图所示出一个 m=18, k=3 （m是该Bit数组的大小，k是Hash函数的个数）的Bloom Filter示例。集合中的 x、y、z 三个元素通过 3 个不同的哈希函数散列到位数组中。当查询元素w时，通过Hash函数计算之后因为有一个比特为0，因此w不在该集合中。

那么怎么判断某个元素是否在集合中呢？同样是这个元素经过哈希函数计算后得到所有的偏移位置，若这些位置全都为1，则判断这个元素在这个集合中，若有一个不为1，则判断这个元素不在这个集合中。

基本语法：

// 创建及修改BloomFilter索引
ALTER TABLE example_db.table_name SET ("bloom_filter_columns" = "column1,column2");

// 查看BloomFilter索引（创建索引的过程是异步的，需要等待一段时间才可以查询到）
SHOW CREATE TABLE table_name;

// 删除BloomFilter索引
ALTER TABLE example_db.table_name SET ("bloom_filter_columns" = "");

使用场景：

• 首先BloomFilter适用于非前缀过滤。
• 查询会根据该列高频过滤，而且查询条件大多是 in 和 = 过滤。
• 不同于Bitmap, BloomFilter适用于高基数列。比如UserID。因为如果创建在低基数的列上，比如 “性别” 列，则每个Block几乎都会包含所有取值，导致BloomFilter索引失去意义。

注意事项：

• 不支持对Tinyint、Float、Double 类型的列建Bloom Filter索引。
• Bloom Filter索引只对 in 和 = 过滤查询有加速效果。
• 如果要查看某个查询是否命中了Bloom Filter索引，可以通过查询的Profile信息查看。

5.3.5 Doris数据查询存在的坑

[1] limit数据查询需要指定order by才能使用

• 问题描述：使用limit offset,size需手动指定order by才能使用，不兼容mysql语法。加上order by之后由于先在大数据上执行order by语句，会导致超时问题。

• 报错信息：1105, 'errCode = 2, detailMessage = OFFSET requires an ORDER BY clause: LIMIT 100000, 100000'

• 解决办法：不使用 limit offset,size 方式限制结果数量，而是根据id这种自增主键列去筛选。

[2] order by默认会自带隐藏的limit 65535

• 问题描述：添加order by后，默认被加上了一个隐藏的limit 65535，导致没有查询出全部结果。

• 解决办法：在order by语句后面手动指定limit，就会覆盖掉这个65535的默认限制。

[3] 一次性查询大量数据把doris服务给查崩了

• 报错信息：have no queryable replicas. err: 88525's backend 10015 does not exist or not alive

• 解决办法：重启doris服务节点即可。

6. 参考资料

[1] 全基因组关联分析(GWAS)神器——PLINK from 知乎 (opens new window)

[2] 种群结构分析--ADMIXTURE和fastStructure from 知乎 (opens new window)

[3] 基因测序技术发展历史及一、二、三代测序技术原理和应用 from 知乎 (opens new window)

[4] 全基因组关联分析学习资料（GWAS tutorial）from 知乎 (opens new window)

[5] 重演一篇傻瓜式GWAS分析|培养兴趣专用 from 知乎 (opens new window)

[6] 浅谈全基因组关联研究 （Genome-wide association study，GWAS）from 知乎 (opens new window)

[7] 用plink做一套GWAS分析 from 知乎 (opens new window)

[8] VCF文件格式说明 from 简书 (opens new window)

[9] 生物基因数据文件——vcf格式详解 from CSDN (opens new window)

[10] 基因组与Python --PyVCF 好用的vcf文件处理器 from CSDN (opens new window)

[11] tabix | 操作VCF文件 from 简书 (opens new window)

[12] bcftools安装及应用实例 from xinzipanghuang (opens new window)

[13] Doris基本介绍 from Doris官方文档 (opens new window)

[14] DataX doriswriter from Doris官方文档 (opens new window)

[15] Apache doris架构及组件介绍 from 知乎 (opens new window)

[16] Python连接Doirs查询实战 from CSDN (opens new window)

[17] Apache Doris (Incubating) 原理与实践 from InfoQ (opens new window)

[18] Doris学习笔记之查询 from CSDN (opens new window)

[19] Bitmap 索引 from Doris官方文档 (opens new window)

[20] BloomFilter索引 from Doris官方文档 (opens new window)

[21] support limit offset, size use default order from Github issues (opens new window)

[22] OFFSET requires an ORDER BY clause from Github issues (opens new window)

[23] 添加order by后，莫名被加上了一个默认的limit 65535 from Github issues (opens new window)

[24] have no queryable replicas. err: 15606's backend 10002 does not exist or not alive from Github issues (opens new window)

[25] 负载均衡 from Doris官方文档 (opens new window)