PostgreSQL遺傳學應用-矩陣相似距離計算(歐式距離,…XX距離)

德哥發表於2017-12-27

SQL矩陣

背景

生物科學中相當重要的工作之一解開遺傳密碼?

歐式空間計算，是其中的一個需求，很有意思吧，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基因配對，產生優良下一代》

在化學分析方面，也有相似的外掛：

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