PostgreSQL 遗传学应用 - 矩阵相似距离计算 (欧式距离,…XX距离)

4 minute read

背景

生物科学中相当重要的工作之一解开遗传密码?

欧式空间计算,是其中的一个需求,很有意思吧,PostgreSQL可以用来解开遗传密码。

https://en.wikipedia.org/wiki/Euclidean_distance

https://www.math.uci.edu/~gpatrick/source/205b06/chapviii.pdf

实际上PostgreSQL是一个扩展性非常强大的数据库,比如在文本相似计算方面,就有诸多扩展插件。

《17种相似算法与GIN索引 - pg_similarity》

https://github.com/eulerto/pg_similarity

https://baike.baidu.com/item/%E6%AC%A7%E5%87%A0%E9%87%8C%E5%BE%97%E5%BA%A6%E9%87%8F/1274107?fromtitle=%E6%AC%A7%E6%B0%8F%E8%B7%9D%E7%A6%BB&fromid=1798948

《PostgreSQL结合余弦、线性相关算法 在文本、图片、数组相似 等领域的应用 - 3 rum, smlar应用场景分析》

《PostgreSQL结合余弦、线性相关算法 在文本、图片、数组相似 等领域的应用 - 2 smlar插件详解》

《PostgreSQL结合余弦、线性相关算法 在文本、图片、数组相似 等领域的应用 - 1 文本(关键词)分析理论基础 - TF(Term Frequency 词频)/IDF(Inverse Document Frequency 逆向文本频率)》

在基因科学方面,也有扩展插件应用:

《为了部落 - 如何通过PostgreSQL基因配对,产生优良下一代》

在化学分析方面,也有相似的插件:

http://www.rdkit.org/

某个生物科技公司,有这样的一种需求:

每张表有几十万行,几万列,全部浮点类型,任意列勾选,计算欧氏距离等需求。

设计

因为数据库设计限制,不能支持一张表几万列,不过PostgreSQL可以将多列存成数组。

1、DNA结构如下:

create table dna (  
  id serial primary key,   -- 主键  
  arr float8[]             -- 浮点数组  
);  

比如每行代表一个物种的测序数据。

2、生成随机浮点数组的函数,可以方便的生成测试数据。

create or replace function gen_randarr(int) returns float8[] as $$  
  select array_agg(random()*1000) from generate_series(1, $1);  
$$ language sql strict;  
postgres=# select gen_randarr(10);  
                                                                                 gen_randarr                                                                                   
-----------------------------------------------------------------------------------------------------------------------------------------------------------------------------  
 {830.968368332833,283.642665948719,64.4483459182084,24.3995497003198,654.509209562093,762.801019474864,109.366949647665,849.462529178709,111.898560542613,650.523159187287}  
(1 row)  
  
Time: 0.758 ms  

3、生成50万条测试数据,每组2万浮点数。

vi test.sql  
  
insert into dna (arr) values (gen_randarr(20000));  
  
pgbench -M prepared -n -r -P 1 -f ./test.sql -c 50 -j 50 -t 10000  

数据大概占用86GB空间。

postgres=# \dt+ dna  
                   List of relations  
 Schema | Name | Type  |  Owner   | Size  | Description   
--------+------+-------+----------+-------+-------------  
 public | dna  | table | postgres | 86 GB |   
(1 row)  

计算欧式距离的函数

可以使用plpgsql创建计算两个浮点数组的欧式距离的函数,长度可以不一样,因为可能不同物种的遗传数据不一样,有的多,有的少。

如果使用C函数,性能会更好。

CREATE OR REPLACE FUNCTION euc_distance(l float8[], r float8[]) RETURNS float8 AS $$  
DECLARE  
  s float8 := 0;  -- 中间结果   
  x float8;       -- LOOP中的数组元素值   
  i int := 1;     -- 数组下标   
  r_len int := array_length(r,1);    -- 右边数组的长度   
  l_len int := array_length(l,1);    -- 左边数组的长度   
BEGIN  
  if l_len >= r_len then  
    foreach x in array l LOOP  
      s := s + ( (x - case when i<=r_len then r[i] else 0 end) ^ 2 );  
      i := i+1;  
    END LOOP;  
  else  
    foreach x in array r LOOP  
      s := s + ( (x - case when i<=l_len then l[i] else 0 end) ^ 2 );  
      i := i+1;  
    END LOOP;  
  end if;  
  RETURN |/ s;  
END;  
$$ LANGUAGE plpgsql;   

例子

postgres=# select euc_distance(array[1,2,3], array[1,2,3]);  
 euc_distance   
--------------  
            0  
(1 row)  
  
Time: 0.386 ms  
postgres=# select euc_distance(array[1,2,3], array[1,2,3,4,5]);  
   euc_distance     
------------------  
 6.40312423743285  
(1 row)  
  
Time: 0.470 ms  

通过这个函数,传入要计算的数组即可计算欧式距离。

计算部分指定位置的欧式距离

这个主要用于部分计算,例如人类和猴子,在某一段的相似性,那么需要从这两条记录中,分别取出要计算的部分,重新组成两个数组,然后计算它们两的欧氏距离。

例子:

select t1.id, t2.id, euc_distance(t1.arr, t2.arr) from   
  (select * from dna where id=1) t1,  
  (select * from dna where id=2) t2;  
  
  
 id | id |   euc_distance     
----+----+------------------  
  1 |  2 | 57768.4024741692  
(1 row)  
  
Time: 12.027 ms  

select t1.id, t2.id,   
euc_distance(  
  array[t1.arr[1], t1.arr[7], t1.arr[8], t1.arr[9], t1.arr[10]], -- 指定位置  
  array[t2.arr[1], t2.arr[7], t2.arr[8], t2.arr[9], t2.arr[10]]  -- 指定位置  
) from   
  (select * from dna where id=1) t1,  
  (select * from dna where id=2) t2;  
  
  
 id | id |   euc_distance     
----+----+------------------  
  1 |  2 | 679.897967241517  
(1 row)  
  
Time: 1.887 ms  

计算被勾选物种的排列组合欧式距离

比如选中了100个物种,计算它们的任意组合的欧氏距离。

需要一些辅助函数:

1、组合去重函数,只去掉重复行。

CREATE or replace FUNCTION has_dupli_val(VARIADIC arr int[]) RETURNS boolean AS $$    
  select count(distinct val)<>count(*) dist_val from unnest($1) t(val) where val is not null;    
$$ language sql strict;   

2、组合去重函数,去掉按列值排序后的重复行。

CREATE or replace FUNCTION arr_sort(arr int[]) RETURNS int[] AS $$    
  select array_agg(id order by id) from unnest(arr) t(id);    
$$ language sql strict;   

3、比如选中了1,2,3,4这四种物种,如何得到他们的排列组合呢?

select distinct on (arr_sort(array[t1.id, t2.id])) t1.id, t2.id from   
  (select unnest(array[1,2,3,4]) id) t1,  
  (select unnest(array[1,2,3,4]) id) t2  
where not has_dupli_val(t1.id, t2.id);    
  
  
 id | id   
----+----  
  1 |  2  
  3 |  1  
  1 |  4  
  2 |  3  
  4 |  2  
  4 |  3  
(6 rows)  
Time: 1.066 ms  

4、创建一个函数,用于计算输入组合物种的排列组合欧式距离。

create or replace function compute_eu_dist(  
  arr_kind int[],   -- 输入物种IDs  
  out kind1 int,    -- 物种1  
  out kind2 int,    -- 物种2  
  out euc_dist float8   -- 物种1,2的欧氏距离  
) returns setof record as $$  
declare  
  l float8[];  -- 左数组  
  r float8[];  -- 右数组  
begin  
  for kind1,kind2 in   
    select distinct on (arr_sort(array[t1.id, t2.id])) t1.id, t2.id from   
      (select unnest(arr_kind) id) t1,  
      (select unnest(arr_kind) id) t2  
      where not has_dupli_val(t1.id, t2.id)    -- 排列组合  
  loop  
    select arr into l from dna where id=kind1;   -- 获取物种1的遗传信息  
    select arr into r from dna where id=kind2;   -- 获取物种2的遗传信息  
    euc_dist := euc_distance(l,r);               -- 计算物种1,2的欧式距离  
    return next;                                 -- 返回  
  end loop;  
  return;  
end;  
$$ language plpgsql strict;  

计算例子:

输入5个物种的ID,返回这5个物种的排列组合欧式距离。

postgres=# select * from compute_eu_dist(array[1,2,3,4,5]);  
 kind1 | kind2 |     euc_dist       
-------+-------+------------------  
     2 |     1 | 57768.4024741692  
     1 |     3 | 57866.2845528097  
     1 |     4 | 57632.9837382263  
     5 |     1 |   57779.36595061  
     3 |     2 | 58004.3926579964  
     4 |     2 | 57593.0783041254  
     5 |     2 | 57802.9690538283  
     3 |     4 | 57837.6707750057  
     3 |     5 | 57921.5524014271  
     4 |     5 | 57818.9181109456  
(10 rows)  
  
Time: 100.582 ms  

并行计算

如果组合确实很多,那么可以使用PG的并行计算,让速度飞起来。

例子如下

alter function euc_distance(float8[], float8[]) PARALLEL safe;


create or replace function get_dna(int) returns float8[] as $$
  select arr from dna where id=$1;
$$ language sql strict parallel safe;

-- 生成1000个物种的排列组合
create table tmp1 as 
select distinct on (arr_sort(array[t1.id, t2.id])) t1.id id1, t2.id id2 from   
  (select unnest(array(select generate_series(1,1000))) id) t1,  
  (select unnest(array(select generate_series(1,1000))) id) t2  
where not has_dupli_val(t1.id, t2.id);   

SELECT 499500
Time: 18591.942 ms (00:18.592)

postgres=# alter table tmp1 set (parallel_workers =48);
ALTER TABLE
Time: 0.602 ms
postgres=# set max_parallel_workers_per_gather =48;
SET
Time: 0.265 ms
postgres=# set parallel_setup_cost =0;
SET
Time: 0.213 ms
postgres=# set parallel_tuple_cost =0;
SET
Time: 0.209 ms

开启了48个并行计算
postgres=# explain select id1,id2,euc_distance(get_dna(id1), get_dna(id2)) from tmp1 limit 10;
                                    QUERY PLAN                                    
----------------------------------------------------------------------------------
 Limit  (cost=0.00..7.55 rows=10 width=16)
   ->  Gather  (cost=0.00..376940.06 rows=499500 width=16)
         Workers Planned: 48
         ->  Parallel Seq Scan on tmp1  (cost=0.00..10119.56 rows=10406 width=16)
(4 rows)

小结

PostgreSQL是一个扩展性很好的数据库,内置了丰富的数据类型。

本例,使用函数编程、数组类型两个特性,解决了生物科学中的遗传计算的场景的疑难问题(上万列,任意组合计算排列组合的欧式距离)。

同时PostgreSQL还能支持并行计算,在重计算的场景,可以提高计算响应速度。如果到达海量数据级别,可以使用HybridDB for PostgreSQL(Greenplum系),实现多机并行。

参考

《17种相似算法与GIN索引 - pg_similarity》

https://github.com/eulerto/pg_similarity

https://baike.baidu.com/item/%E6%AC%A7%E5%87%A0%E9%87%8C%E5%BE%97%E5%BA%A6%E9%87%8F/1274107?fromtitle=%E6%AC%A7%E6%B0%8F%E8%B7%9D%E7%A6%BB&fromid=1798948

《为了部落 - 如何通过PostgreSQL基因配对,产生优良下一代》

《PostgreSQL结合余弦、线性相关算法 在文本、图片、数组相似 等领域的应用 - 3 rum, smlar应用场景分析》

《PostgreSQL结合余弦、线性相关算法 在文本、图片、数组相似 等领域的应用 - 2 smlar插件详解》

《PostgreSQL结合余弦、线性相关算法 在文本、图片、数组相似 等领域的应用 - 1 文本(关键词)分析理论基础 - TF(Term Frequency 词频)/IDF(Inverse Document Frequency 逆向文本频率)》

Flag Counter

digoal’s 大量PostgreSQL文章入口