# 13章　TabixとSQLite: メモリを使わないアプローチ
- メモリ上ではなく、ディスク上にデータを保存して作業を行う計算上の戦略
- 計算は遅くなるが、メモリ内処理やストリーミング処理が適切でないときにとるべき手段
- [13章Github](https://github.com/vsbuffalo/bds-files/tree/master/chapter-13-out-of-memory)

## 13.1 BGZFとTabixを使用した索引付きのタブ区切りファイルへの高速アクセス
- バイオインフォマティクスではランダムアクセスが必要
- ランダムアクセス：[例](https://wa3.i-3-i.info/word19512.html#:~:)
- ファイルシステムのインデックスを参照してデータの場所を算出し、直接そこにアクセスする方法
- 逆に端から順にアクセスしていく手法をシーケンシャルアクセスというらしい
- `BGZF(Blocked GNU Zip Format)`: ランダムアクセスをしたい大きなファイルを圧縮するときに使われる。普通の圧縮よりは少しサイズが大きくなる
- 例えば遺伝子アノテーションファイル(gtf)など
- 一定のサイズごとにファイルを切り分けて分割して圧縮
- [日本語でのちょっと詳しめなBGZFの説明](https://qiita.com/kojix2/items/d8a66703e0c63d017056#:~:)
- `Tabix`:　タブ区切りファイルにインデックスを張るソフトウェア
- bgzfに圧縮されたファイルに、Tabixを使ってインデックスをはる

### ゲノムの場所や範囲にリンクされたデータへアクセスする際の問題点
1. データが完全にメモリにのらない可能性がある
2. RDBを使っても、何百万ものエントリーがあった場合には計算が遅くなる

> これらの問題に対応するのに`BGZF`や`Tabix`は有効(高速ランダムアクセスを可能にするのに有用!)

### 13.1.1 bgzipを使ってTabixのためにファイルを圧縮する
- gzipされたマウスの`gtf`ファイルを例題として使う
1. `gzcat`でファイルのヘッダーを確認
- `gzcat ファイル名(Mus_musculus.GRCm38.75.gtf.gz)`
2. ヘッダー以外の箇所を、染色体と開始位置でソートして、`bgzip`で圧縮しなおす
- `zgrep -v`でヘッダー以外の場所を選択
- `sort`でソート
- `bgzip`で圧縮
- サブシェルを使って、ヘッダーとそれ以外をそれぞれ圧縮にパイプしている
- `(zgrep "^#" Mus_musculus.GRCm38.75.gtf.gz; zgrep -v "^#" Mus_musculus.GRCm38.75.gtf.gz | sort -k1,1 -k4,4n) | bgzip > Mus_musculus.GRCm38.75.gtf.bgz`

### 13.1.2 Tabixによるファイルの索引作成
- Tabixコマンドラインツールを使って、BGZF圧縮ファイルのインデックスを作成する
- `-p`オプションでファイル形式を選択
- カスタムのタブ区切りの場合(bgzでない形式を利用する場合)は-pオプションなしで実行することもできるが、、
- 基本的にはbgzのようなバイオインフォマティクスの形式を使うことが推奨
- `tabix -p gff Mus_musculus.GRCm38.75.gtf.bgz `

### 13.1.3 Tabixを使う
- これで高速ランダムアクセスが可能になった
- tabixコマンドラインツールでインデックスをつけたファイルの検索もできる
- tabixは、HTTPやFTPサーバーに対しても機能する
- 例えばbgzipしてTabixでインデックスをつけたファイルをサーバーにおいておけば、色んな人がそのファイルにリモートアクセス可能
- 例: 16番染色体の23146536から23158028までの特徴領域にアクセスしたいとき

tabix Mus_musculus.GRCm38.75.gtf.bgz 16:23146536-23158028 | head -n3
16      protein_coding  UTR     23146536        23146641        .       +       .       gene_id "ENSMUSG00000022878"; transcript_id "ENSMUST00000023593"; gene_name "Adipoq"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "Adipoq-001"; transcript_source "ensembl_havana"; tag "CCDS"; ccds_id "CCDS28075";
16      protein_coding  exon    23146536        23146641        .       +       .       gene_id "ENSMUSG00000022878"; transcript_id "ENSMUST00000023593"; exon_number "1"; gene_name "Adipoq"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "Adipoq-001"; transcript_source "ensembl_havana"; tag "CCDS"; ccds_id "CCDS28075"; exon_id "ENSMUSE00000131550";
16      protein_coding  gene    23146536        23158028        .       +       .       gene_id "ENSMUSG00000022878"; gene_name "Adipoq"; gene_source "ensembl_havana"; gene_biotype "protein_coding";

- その中でもexonの行だけを取りたい場合
- awkを使う

tabix Mus_musculus.GRCm38.75.gtf.bgz 16:23146536-23158028 | awk '$3 ~ /exon/ {print}'
16      protein_coding  exon    23146536        23146641        .       +       .       gene_id "ENSMUSG00000022878"; transcript_id "ENSMUST00000023593"; exon_number "1"; gene_name "Adipoq"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "Adipoq-001"; transcript_source "ensembl_havana"; tag "CCDS"; ccds_id "CCDS28075"; exon_id "ENSMUSE00000131550";
16      protein_coding  exon    23146646        23146734        .       +       .       gene_id "ENSMUSG00000022878"; transcript_id "ENSMUST00000171309"; exon_number "1"; gene_name "Adipoq"; gene_source "ensembl_havana"; gene_biotype "protein_coding"; transcript_name "Adipoq-002"; transcript_source "havana"; exon_id "ENSMUSE00000875966"; tag "cds_end_NF"; tag "mRNA_end_NF";

- 使わない？？
- 一般的にはGalaxyとか使った方がいい？
- bcltoolを使うときにもしかしたら使う？？

## 13.2 SQLiteを使ったリレーショナルデータベース(RDB)の操作

- これまで紹介されてきたフラットファイル形式のデータは、単純で柔軟性があり、ポータブルなのでバイオインフォマティクスで広く使われている

- しかしフラットファイル自体は表同士の関係を持っていないため、データ同士を適切に関係付けて格納および操作するには、**リレーショナルデータベースが役に立つ**
- `SQL`という言語を利用

- `SQLite`は、リレーショナルデータベース管理システム
- 導入が簡単
- ギガバイト単位のデータベースであれば、うまく機能する

### 13.2.1 バイオインフォでRDBが必要になるときは?

- RDBとの直接的なやりとりが発生する頻度は少なめ?
- 理由1: 大抵の大規模RDBはAPIを提供している
- 理由2: バイオインフォの作業ではRDBの使用が最良ではないことが多い
- キュレーションされた遺伝子モデルのセット > RDBの活用が役に立つ
- データの追加とマージにおいては、既存のレコードを照会しつつ重複データを防ぎながら行う必要がある。RDBはそれを可能にするデータ追加メソッドがある
- データの更新についても、過去の遺伝子モデルを安全にアーカイブできるし、データ間のリレーションの更新を簡素化
- データの操作について、特定の列にインデックスをつけれるので大幅に高速化できる
- シークエンシング実験から得られたアラインメントデータのセット > RDBは不要
- 基本的に生の実験データを使う場合は、ストリーミングアルゴリズムで事足りるし高速
- それぞれのファイルの専門のツールを使うのが計算効率が高い
- samtoolsやbcftoolsなど

### 13.2.2 SQLiteのインストール
- 自分は`conda`でインストールしました
- 本では`Homebrew`や`apt-get`でのインストール方法が紹介されています

### 13.2.3 CLIによるSQLiteの使い方
- `sqlite3 データベース名`: データベース作成
- `.table`： 全テーブルを表示
- `.schema`: テーブル作成に使用されたSQLのCREATE文を出力する
- 列には型がなく、格納されるデータに型がある
- text型(文字列型)
- integer型(整数型)
- real型(実数型)
- NULL型(データの欠落または値なしの場合)
- BLOB(Binary Large Object)型: 任意のデータをバイト列として格納する
- 列をすべて同じ型にする必要はない、が、
- 下流での処理を楽にするために、同じ列に格納されるデータは全て同じ型を保つ方がいい

### 13.2.4 データの操作
- ここからは使い方の説明
- データ検索の基本形
`select * from テーブル名;`
- SQLite内部に入らずに欲しいデータを出力する

sqlite3 gwascat.db "select * from gwascat" > results.txt
head -n 1 results.txt

- limitで表示するエントリー数を制限する
`select * from gwascat limit 2`

select trait,chrom,position,strongest_risk_snp,pvalue　from gwascat limit 5;

- 結果の表示を見やすくする
`.header on`, `.mode column`

select trait,strongest_risk_snp,pvalue　from gwascat order by pvalue desc limit 5;

trait                                       strongest_risk_snp  pvalue    
------------------------------------------  ------------------  ----------
Serum protein levels (sST2)                 rs206548            90000000.0
Periodontitis (Mean PAL)                    rs12050161          4000000.0 

- whereを使って検索条件づけをする

select chrom, position, trait, strongest_risk_snp,　pvalue　from gwascat where chrom = "22" and pvalue < 10e-15;

### 13.2.5 SQLite関数

- 既存の列から新しい列を作成するときに、`as`を使う
- traitの最初の文字を小文字にして、asで新たな列にするには

sqlite> select lower(trait) as trait,
   ...> "chr" || chrom || ":" || position as region from gwascat limit 5;
trait                 region        
--------------------  --------------
asthma and hay fever  chr6:32658824 
asthma and hay fever  chr4:38798089 
asthma and hay fever  chr5:111131801
asthma and hay fever  chr2:102350089
asthma and hay fever  chr17:39966427

- 全てのNULLをNAに置き換えるには、`ifnull`関数を使う

sqlite> select ifnull(chrom, "NA") AS chrom, ifnull(position,"NA") AS position,
   ...> strongest_risk_snp, ifnull(pvalue,"NA") AS pvalue FROM gwascat
   ...> WHERE strongest_risk_snp = "rs429358";
chrom  position  strongest_risk_snp  pvalue 
-----  --------  ------------------  -------
19     44908684  rs429358            5.0e-14
19     44908684  rs429358            1.0e-06
NA     NA        rs429358            NA  

### 13.2.6 SQLite集約関数
- 集約関数は、SQL文から返された列全体を入力として受け取り、単一の値を返すもの
- `count(*)`: 列の行数を数える
- `avg(*)`: 列の平均を数える
- その他は本を参照、またはグーグル参照

select "2007" as year, count(*) as number_entries from gwascat where date between "2007-01-01" and "2008-01-01";

select count(distinct strongest_risk_snp) as　unique_rs from gwascat;

- `group by`: 列をいくつかの列の値でグループにまとめる

select chrom, count(*) as nhits from gwascat group by chrom　order by nhits desc;

select strongest_risk_snp, count(*) as count from gwascat group by strongest_risk_snp order by count desc limit 5;

select strongest_risk_snp, strongest_risk_allele, count(*) as count from gwascat group by strongest_risk_snp, strongest_risk_allele  order by count desc limit 10;

- substr関数は、文字列の中の指定した位置から指定した長さの部分文字列を取得できる

select substr(date,1,4) as year from gwascat group by year

select substr(date,1,4) as year, round(avg(pvalue_mlog),4) as mean_log_pvalue, count(pvalue_mlog) as n from gwascat group by year

select substr(date,1,4) as year, round(avg(pvalue_mlog),4) as mean_log_pvalue, count(pvalue_mlog) as n from gwascat group by year having count(pvalue_mlog) > 10;

### 13.2.7 サブクエリー

- それぞれの論文がもつ形式関連の数を数える

sqlite> select substr(date,1,4) as year,author,pubmedid, count(*) as num_assoc from gwascat group by pubmedid limit 5;
year  author         pubmedid  num_assoc
----  -------------  --------  ---------
2005  Klein RJ       15761122  1        
2005  Maraganore DM  16252231  1        
2006  Arking DE      16648850  1        
2006  Fung HC        17052657  3        
2006  Dewan A        17053108  1 

- 先程のクエリーを括弧で囲い、その部分がFromで指定するデータベースであるかのように扱う(サブクエリー)
- これによって、年ごとの形質関連の形質関連のエントリー数平均を表示する

sqlite> select year, avg(num_assoc) from (select substr(date,1,4) as year, author, count(*) as num_assoc from gwascat group by pubmedid) group by year;
year  avg(num_assoc)  
----  ----------------
2005  1.0             
2006  1.6             
2007  5.87837837837838
2008  7.64566929133858
2009  6.90625         
2010  9.21660649819495
2011  7.4968152866242 
2012  13.4536741214058
2013  16.605504587156 

- 直近の研究の方が形質関連のエントリーが平均して増えてきていることがわかる

### 13.2.8 リレーショナルデータベースの編成と結合(Join)
- `PubMedID`が`24388013`の研究に関わるエントリーを確認する

sqlite> select date, pubmedid, author, strongest_risk_snp from gwascat where pubmedid = "24388013" limit 5;
date        pubmedid  author       strongest_risk_snp
----------  --------  -----------  ------------------
2013-12-30  24388013  Ferreira MA  rs9273373         
2013-12-30  24388013  Ferreira MA  rs4833095         
2013-12-30  24388013  Ferreira MA  rs1438673         
2013-12-30  24388013  Ferreira MA  rs10197862        
2013-12-30  24388013  Ferreira MA  rs7212938   

- いくつかの列の値が重複してしまう >　テーブルの不要な冗長性(`スプレッドシート症候群`)が生まれる　> ディスク容量が不必要に増加する
- これを避けるために、`データベースの正規化`を行う
- データベースの正規化とは、正規形のレベルを階層的に上げていく操作？？
- データの重複をなくし整合的にデータを取り扱えるようにデータベースを設計すること
- データの追加、更新、削除に伴うデータの不整合や損失をふせぎ、メンテナンスの効率を高める
- 正規形とは?
- 行の各列には、*1つ*のデータ値のみが含まれるようにデータを整理することが最善
- 1つのフィールドに複数の値を持つ行がある場合は、複数の行に分割する(`第一正規形`)
- あるフィールド値に重複がある場合には、テーブルを分割する（分割されたテーブルを`第二正規形`、`第三正規形`という)

#### 内部結合
- 結合述語: ２つのテーブルで共通する列
- `true`であれば、2つのテーブルで行数が等しい
- `false`であれば、2つのテーブルで行数が不揃い
- 内部結合は、結合述語が`true`の時に使う
- ここでは第二正規形としたテーブルを左のテーブルといい、第三正規形としたテーブルを右のテーブルと呼ぶ
- 今回の例の結合述語は`study_id`
- `select * from 左のテーブル join 右のテーブル on 左のテーブル名.結合述語 = 新たなテーブル名`

#### 外部結合
- 結合述語が`false`の時に使う
- SQLiteでは左外部結合のみサポート
- つまり、左のテーブルの全ての行に対して、**可能な限り**右のテーブルの行を結合する作業
- `select * from 左のテーブル left outer join 右のテーブル on 右のテーブル名.結合述語 = 新たなテーブル名`
- 右のテーブルで結合できなかった行は`NULL`が入る

### 13.2.9 データベースへの書き込み
#### テーブル作成
- 作成する全てのテーブルには主キー列(テーブル内の行を識別するために使用される一意の整数)が必要

例
sqlite> create table テーブル名(
   ...> id integer primary key,
   ...> 列名1 text,
   ...> 列名2 integer,
   ...> 列名3 integer,
   ...> 列名4 text,
   ...> 列名5 text);

#### 行を挿入する

sqlite> insert into テーブル名(id,chrom,start,end,strand,name)
   ...> values(NULL,"16",48224287,48224287,"+","rs17822931");
sqlite> select * from variants;
1|16|48224287|48224287|+|rs17822931

#### インデックス作成
`create index インデックス名 テーブル名(列名)`
- 注意点としてはディスク容量をくうこと
- テーブルが非常に大きい場合はインデックスも大きくなる場合があるので注意
- インデックスの削除等は本を参照

### 13.2.10 テーブルの削除とデータベースの削除
- テーブルの削除は、`drop table テーブル名:`
- 全テーブル削除する場合は、`db`ファイルを削除する

### 13.2.11 Pythonを使ってSQLiteを使う

### 13.2.12 データベースのダンプ
- 指定されたデータベースを新規に完全に再構築するために必要な一連のSQLコマンド
- データベースのバックアップや複製に役立つ
- ダンプする場合は`sqlite3 dbの名前 ".dump" > dump.sql`で対応
- dumpコマンドにテーブル名をつければ単一のテーブルだけをダンプできる
- ダンプファイルを使ってデータベースを再現する場合には、`sqlite3 新たなdbの名前 < dump.sql`
- 頻繁にダンプしといた方がいい