[plink][genetics]遺伝子解析ソフトplink tutorial

初めに

「ピーリンク」ではなく「プリンク」と呼ぶ。Population-based linkageの略。
https://www.cog-genomics.org/plink/

遺伝学の大家Purcellが作成。
https://en.wikipedia.org/wiki/Shaun_Purcell

なぜ特殊なデータ形式が必要なのか

通常のデータセットであれば、被験者が\(1\)万人いて、変数が\(5,000\)個あったとすると、\(10,000\times 5,000\)のデータフレームができる。これはこれで十分大きいが、遺伝子データの場合、\(1\)万人の被験者に対して\(2,200\)万の遺伝子変数がある。マイクロソフトのエクセルの読み込み上限は\(1,048,576\)行、\(16,384\)列なので、到底読み込めない。遺伝子変数を含むデータはギガバイト、時にはテラバイトにいもなるので、以下に記すような特殊な形式が必要になるというわけである。

Plinkで用いられるデータの種類

.pedと.mapはテクストファイルで、残りはすべてバイナリーファイルとなっている。バイナリーファイルの方が読み込みが早く利便性が高い。また、共変量のファイルも必要となる。

具体例、例えば妊娠糖尿病の患者について考えてみる。.bedファイルにはすべての患者の遺伝子型が含まれる。この場合、ケースコントロールスタディーとすると、妊娠糖尿病の患者と、健康な人（対照群）のすべてのデータが含まれる。.famファイルには個人に関係するデータ（データ上の他の個人との家族の相互関係、性別、妊娠糖尿病の診断など）が含まれる。.bimファイルにはSNPsの実際の遺伝子上の位置が含まれ、その他必要となる共変量を含むファイルも使用される。また、ファイル上の変数の説明も以下に記載する。

.pedはPedigree（血統）の略である。各行は個人に相当し、最初の\(6\)列は個人情報を提供している。ファイルはヘッダーや変数名を含まない。最初の\(4\)列はFID, ID, F, Mとなっているが、常に含まれるわけではない。IDのみの場合も多い。\(5, 6\)列は性別と対象の表現型を含む。残りの列は遺伝子情報を表す。例えば$$\begin{eqnarray}Ch123456\ \ N123456\ \ 0\ \ 0\ \ 2\ \ 1\ \ G\ \ G\ \ \cdots \\ Ch123457\ \ N123457\ \ 0\ \ 0\ \ 1\ \ 1\ \ G\ \ G\ \ \cdots \\ Ch123458\ \ N123458\ \ 0\ \ 0\ \ 2\ \ 2\ \ C\ \ G\ \ \cdots \\ \vdots\end{eqnarray}$$というデータがあるとき、ID 123456の人は最初のSNPにGGというアレルを有し、ID 123458の人はCGというSNPを有する。なので、SNPの数が\(K\)個のとき、.pedファイルは\(6 + (2K)\)列を有する。.pedファイルはテキスト形式なので、通常のテキストアプリケーションで開くことができるが、多数のSNPがあるときにはパッと見て判読することは困難だろう。

.pedファイルは.mapファイルとペアになっている。.mapファイルは個人の遺伝子型についての完全な情報を提供する。.mapファイルはどのSNPがジェノタイピングされ、ゲノム上でどのように位置づけられるかの情報が記載されている。

.mapファイルの先頭列は染色体の番号を表す。\(2\)列目はSNP IDで、各々のSNPにはrs\(\cdots\)という番号が付与されている。\(3, 4\)列目にはSNPの位置が記載される。\(3\)列目はセンチモルガンで、\(1\)センチモルガンは\(100\)回の減数分裂につき\(1\)回の遺伝子組み換えが起こることを意味する。\(4\)列目は塩基対座標あるいは塩基対の遺伝的距離が記載される。例えば、バリアントの間の分子の数などがこれに相当する。\(1\)センチモルガンというのは人間では\(100\)万塩基対に相当する。染色体ごとの総センチモルガンはHuman Reference Genomeが基準となっている。何をreferenceするかによってSNPの位置が変わることは注意すべきである。結局、.mapファイルはSNPの数が\(K\)個のとき\(4\)列\(K\)行のデータセットとなる。

plink binary file

.mapと.pedファイルは通常のテキストエディタで開けるが、遺伝子解析ではデータ容量が膨大になるので、データ容量が不足する。こういった場合はバイナリーファイルでデータを保存するのが良く、PLINKでは.bed, .fam, そして.bimが用意されている。

.bedファイルはテキストエディタで読み込めず、遺伝子情報を圧縮して含んでいる。

.famファイルは.pedファイルと同様に個人の情報を含んでいる。最初の\(6\)列は.pedファイルと同様である。

.bimファイルはSNPの情報を含んでいる。これは.mapファイルに相当する。PLINKでは.map/.pedファイルをバイナリーファイルに変更したり戻したりするのは容易である。

Imputeされた遺伝子情報データ

Imputationとは欠損値の推定・代入である。読み取られる遺伝子情報とreference panelのLD (Linkage Disequilibrium)の情報から欠損している情報を挿入する。reference panelは例えば1000 Genome ProjectやHaplotype Reference Consortiumといった既存のSNPデータベースである。Imputateされた遺伝子情報はImputation probabilityと関連している（しばしばgenotype callsと呼ばれる）。これはある遺伝子情報がどのくらいreference panelに基づいているかを示す。例えば、rs2777888というSNPがCCである確率が28%、TCである確率が52%、そしてTTである確率が21%、という具合である。Imputationによる不確実性を排除するには、最も確率が高い配列のみを選択するというものである。これは妥当な選択であり、しばしばポリジェネックスコアを計算するときに採用される手法である。PLINKではこのImputation probabilityを扱うことができる。

PLINK 2.0 (2018年にリリースされた)ではImputationに対応している。以下の\(3\)つのファイルが追加されている。

• .pvarファイル：.bimファイルと同様に遺伝子マーカーについての情報を記載している。
• .psamファイル：.famファイルと同様にサンプル中の個人情報について記載している。
• .pgenファイル：これはテキストエディタでは読み込めず、圧縮されたgenotype probabilityを収納している。

PLINK 2.0の.psamファイルには遺伝子変異についてのより多くの情報が含まれている。配列情報だけでなく、Imputationの質についても記載されている。これはQUALあるいはINFOという列に記載されている。.pgenファイルには遺伝的変異がどの程度の確率で特定の遺伝子型を持つかに関する情報が記載されている。この情報はQuality control作業を行うときに使用される。

PLINK 2.0のファイル形式について

.pvarファイルは以下のようになる。

.psamファイルは以下のようになる。

Oxfordファイル形式

.bgenとかのフォーマット。日本語の情報が少ないのでここに記載する価値はあるかと。Oxford大学とWellcome Trust Centreのグループが開発。GTOOLやSNPTESTなどで使われる。.pedや.mapと同じように、\(2\)つのファイルに遺伝子情報がストアされる。genotypeファイルはSNPごとに行に遺伝子情報が掲載され、各列は個人を表す（これはPLINKの.bedファイルと逆になっている）。Oxford形式のgenotypeファイルは以下のようになる。

plinkのbinrary形式に変換するときは、PLINK2を用い、

plink2 --bgen bgenfile.bgen "ref-first" --sample samplename.sample --out outfile

とする。この”ref-first”と”ref-last”に関しては、色々調べたが、とりあえず”ref-first”としとけば問題がないようである。

plink2 bgen to vcf ukbiobank

www.biostars.org

実行

以下ではLink先のTutorialを実行する。
https://zzz.bwh.harvard.edu/plink/tutorial.shtml

インストールは済んでいるものとする。MacのTerminalで作業ディレクトリに移動しておく。

$ ./plink

で起動。以下$は省略する。

plink --file hapmap1

とすることで、オプションなしでhapmap1ファイルを開く。

--file

オプションで開くファイルを指定できる。この指定を行うことで、.pedファイルと.mapファイルを探す。つまり、hapmap1.pedとhapmap1.mapが同一の作業ディレクトリに存在する必要がある。別々のフォルダにある場合は

--ped
--map

オプションを指定すれば良い。pedファイルは”pedigree”の略である。1行につき一人の遺伝情報を含む。mapファイルはpedファイル中のマーカーの名前及び位置情報を記す。pedファイルもmapファイルもプレーンテキストなので、付属のエディタなどで開くことができる。pedファイルの構成は

• #1 Family ID
• #2 Individual ID
• #3 Paternal ID
• #4 Maternal ID
• #5 Sex (1 = male, 2 = female, other = unknown)
• #6 Phenotype

の6列となる。重要なのは、一つのpedファイルは一つのphenotypeのみを含むことである。それが量的データなのか、質的データなのかをplinkは自動で判断する。プレーンテキストとして開くと、下のようになっている。

それぞれのFamily IDに対して、各個人のphenotypeを表す。mapファイルの構成は

• #1 chromosome (1-22, X, Y, or 0 if unplaced)
• #2 rs# or snp identifier
• #3 Genetic distance (morgans)
• #4 Base-pair position (bp units)

の4列となる。

Binary形式

最初にファイルをbinary形式のbedファイルに変換する。こうすることでメモリー空間を節約し解析をスピードアップできる。

plink --file hapmap1 --make-bed --out hapmap1

と実行することで、bed, bim, famファイルが作成される。bim, famファイルはテキストエディタで開くことができるが、bedファイルは編集できない。例えば、少なくとも95%以上のgenotypingを持つ個人のみをファイルに含めたい時は、

plink --file hapmap1 --make-bed --mind 0.05 --out highgeno

とする。ここで、mindはMaximum per-person missingで、個体において遺伝子型が決定されなかった割合の最大値であるこの例では5%以上遺伝子型が不明なサンプルは除外される。

最初に読み込むファイルがped, map形式ではなく、binarayファイルである時、オプションは

plink --bfile hapmap1

とする。

データのクオリティコントロール

plink --bfile hapmap1 --missing --out miss_stat

-missingオプションを付けることで、被験者一人ずつのgenotypingの結果とsnpごとのgenotypingの欠損値を集計し、データのクオリティを評価することができる。上記を実行することで、imissファイルとlmissファイルが生成される。imissのiはindividual, lmissのlはlocusの略である。imissのフォーマットは

• #1 FID: Family ID
• #2 IID: Individual ID
• #3 MISS_PHENO: Missing phenotype (Y/N)
• #4 N_MISS: Number of missing SNPs
• #5 N_GENO: Number of non-obligatory missing genotypes
• #6 F_MISS: Proportion of missing SNPs

となっている。ここで、obligatory missingとはお決まりの欠損値、個人ごと、SNPごとに90%のgenotypeが欠損しているということである。なお、ここではデフォルトフィルタリングによって、個人の最大欠損数(Maximum per-person missing, MIND)が0.1(10%)以上の個人は自動的に解析から除外されている。また、 Proportion of missing SNPs = Number of missing SNPs/Number of non-obligatory missing genotypesである。この値が3-7%の個人はDNAのクオリティに問題があるとして除外する。lmissのファイルの構成は

• #1 CHR: Chromosome Number
• #2 SNP: SNP identifier
• #3 N_MISS: Number of Individual missing this SNP
• #4 N_GENO: Number of non-obligatory missing genotypes
• #5 F_MISS: Proportion of sample missing for this SNP

となっている。ここでもデフォルトのフィルタリングにより、SNPの最大欠損数(Maximum per-SNP missing(GENO))が0.1以上のSNPは解析から除外されている。

統計量の記載

plink --bfile hapmap1 --freq --out freq_stat

とすることで、それぞれのSNPについてminor allele frequencyとallele codeを含んだファイルを生成できる。withinオプションを指定することで、phenotypeごとの解析も可能である。

plink --bfile hapmap1 --freq --within pop.phe --out freq_stat

.pheファイルは3列からなり、順にFamily ID, Individual ID, phenotypeを表す。この場合phenotypeは1が中国人、2が日本人ということである。このコマンドを実行すると、freq_stat.frqの代わりに、freq_stat.frq.stratが生成される。特定のSNPについてpopulation内での頻度を知りたいのであれば、以下のように

plink --bfile hapmap1 --snp rs# --freq --within pop.phe --out snp1_frq_stat

とする。#にrs番号を記載する。

基本的な解析

plink --bfile hapmap1 --assoc --out as1

とすることで、それぞれの各SNPごとに疾患との基本的な関連解析を行う事ができる。出力されるas1ファイルの中身は以下である。

• #1 CHR: Chromosome
• #2 SNP: SNP identifier
• #3 A1: Code for allele1 (the minor, rare allele based on the entire sample frequencies)
• #4 F_A: The frequency of this variant in cases
• #5 F_U: The frequency of this variant in controls
• #6 A2: Code for the other allele
• #7 CHISQ: The chi-squared statistic for this test (1 df)
• #8 p: The asymptotic significance value for this test
• #9 OR: The odds ratio for this test

多重補正された値をソートするには

plink --bfile hapmap1 --assoc --adjust --out as2

とする。このファイルの構成は以下のようになる。

• #1 CHR: Chromosome
• #2 SNP: SNP identifier
• #3 UNADJ: Unadjusted, asymptotic significance value
• #4 GC: Genomic control adjusted significance value.
• #5 BONF: ボンフェローニ補正されたsignificance value
• #6 HOLM: Holm step-down補正されたsignificance value
• #7 SIDAK_SS: Sidak single-step補正されたsignificance value
• #8 FDR_BH: Benjamini & Hochberg step-up補正されたFDRコントロール
• #9 FDR_BY: Benjamini & Yekutieli step-up補正されたFDRコントロール

adjustコマンドを用いると、logファイルにinflation factorとmean chi-square factorが記載される。inflation factorは帰無仮説下では1となることが期待されるが、1を上回る場合、population stratificationの補正が無いと、偽陽性となる可能性がある。つまり、人種によるsnpの頻度の違いが疾患と関連があるようにみえてしまうということである。

特定のSNPの抽出

PLINKだけでなく他のソフトウェアで解析するため、あるSNPの情報をbinaryファイルでなくテキスト形式ファイルで読み込みたいときは、

plink --bfile filename --snp rsxxxxx --recodeAD --out filename2

とすると良い。.RAWファイルが出力されるので、Rなどで解析できる。

データQuality Control(QC)について（その2）

個人のSNPについてのQCと、集団全体のについてのQCが必要である。個人の場合、

plink --bfile merge --mind 0.05 --make-bed --out merge2

とすることで、\(5\%\)以上SNPが決定できなかった個人はサンプルから除外される。

集団の場合、集団の中で何\(\%\)以上遺伝子型が決定されなかったSNPを除くという操作が可能であり、それが以下である。

plink --bfile merge --geno 0.05 --make-bed --out merge2

とすることで、call ratesが\(95\%\)未満のSNPはサンプルから除外される。低いcall ratesとは、遺伝子型の欠失率が高い状況、言い換えれば多くの個体で欠失したSNPがある状況を指す。この文脈での “calling “は、1つのユニークなSNPの推定を表すために使用される。通常、call ratesが\(95\%\)未満のマーカーは解析から除外される。