生信小道具:Plinkの実戦演習(2)
6454 ワード
前回の内容に引き続き、Plinkとその基本的なフォーマットを紹介した後、もちろん本格的な時になって、実戦の練習をして、以下は例を通じて仲間たちと一緒にPlinkの使用を熟知します.
公式チュートリアルで使用されているテストソフトウェアをダウンロードします.
テストファイルを解凍して、中に何があるか見てみましょう.
このテストセットには、次の4つのファイルが含まれています.
これらの試験ファイルを簡単に紹介する:このhapmap太字テキストファイルには89個体からランダムに選択された遺伝子型(約80000個の常染色体SNP)が含まれている.qt.phe:数量性状は単独の代替代表型ファイルにあります.pop.pheファイルには仮想表型が含まれており、中国人の個体コードは1、日本人コードは2である.集団間の違いを調べるために使用します.
このコマンドが完了すると、次の情報も表示され、この変異ファイルには83534と89個の個体が含まれていることが簡単にわかります.
生成されたファイルをバイナリPEDファイルとしてさらに圧縮することができます.
たとえば、高周波遺伝子分型(少なくとも95%以上の完全性)を持つ個体のみを含む新しいファイルを作成する場合は、次のように実行できます.
準備が終わった後、SNPのデータを統計し始めました.まず、損失率
ここでは、ファイルをバイナリbed形式に変換したので、
2つのファイルが生成されます.
1つのファイルは
ここではSNP毎に失われた個体数(N＿MISS)と失われた個体の割合(F＿MISS)を示した.別のファイルを見てみましょう.
ここでは個々について示したが、失われたSNPの合計数(N＿MISS)と、各個体における失われたSNPの合計数(F＿MISS)の割合(F＿MISS)は、一般的に、この個体における失われたSNPの合計数(N＿MISS)の割合が高い場合、後続の分析の前に、その個体を除去することを選択する.ここでは,このテストデータの品質は悪くなく,個体のSNP損失の割合は非常に低い.
また、ある染色体の損失率にのみ興味がある場合は、対応する染色体の損失率情報を出力するように指定することもできます.
ここでfreq_が生成されますstat.frqファイルは、SNP毎の二次等位遺伝子周波数(MAF)を含む.
ここでも指定されたSNPに対して、その等位遺伝子の頻度を調べることができます
ここでMAFといえば、簡単にMAFとは何かを拡張しますが、なぜMAFでフィルタリングするのですか?
MAFは、副等位遺伝子の頻度、すなわち、頻度の低い等位遺伝子の頻度である.MAFは様々な理由で濾過することができる.そのうちの1つは,低い被覆量または深さを持つ配列では,非常に低いMAFを持つ試料が多くの騒音を発生し,偽陽性の下流分析を引き起こす可能性があることである.極端な例を挙げると、1つの個体にのみ存在する等位遺伝子(非常に低い等位遺伝子周波数)は、集団構造推定に有用で信頼できる情報を追加しない.データが全ゲノム関連研究(GWAS)に使用される場合、非常にまれな等位遺伝子を有意に記述するために非常に強い統計的能力が必要であるため、通常、比較的高いMAF閾値が用いられる.一般的によく使われるバルブ値は0.05と0.01です.しかし、最適なMAFフィルタしきい値はあなたのデータ型に依存し、計画された分析、データ品質、サンプルサイズも統計分析で計算できます.
ここでは、バルブ値0.05を使用してフィルタリングする例を示します.
濾過後、24114変異部位が濾過され、59420変異部位が保持されることがわかる.
MAFに加えて、
遺伝子型喪失率が10%より大きい個体サンプルを除去する:
95%を超える個体が持っていない変異部位を除去する.
もちろん、すべてのフィルタリングは同時に実行できます.
最後に簡単に説明すると、PlinkはLDの計算にも使用できます.LDの計算についてもっと詳しく議論します.
もちろん、いくつかのRパッケージ、例えば
ここまで私が日常的に使っているPlink機能を紹介しましたが、これらの機能は実はPlinkの氷山の一角にすぎず、その強さはmanualをよく読んでこそ深く体得することができます.
テストデータのダウンロード
公式チュートリアルで使用されているテストソフトウェアをダウンロードします.
wget http://zzz.bwh.harvard.edu/plink/hapmap1.zip
テストファイルを解凍して、中に何があるか見てみましょう.
unzip hapmap1.zip
rm hapmap1.zip
ls -thrl
このテストセットには、次の4つのファイルが含まれています.
-rw-rw-r-- 1 hhu hhu 1693668 Jun 6 2006 hapmap1.map
-rw-rw-r-- 1 hhu hhu 29739617 Jun 8 2006 hapmap1.ped
-rw-rw-r-- 1 hhu hhu 979 Jun 8 2006 pop.phe
-rw-rw-r-- 1 hhu hhu 1869 Jun 8 2006 qt.phe
これらの試験ファイルを簡単に紹介する:このhapmap太字テキストファイルには89個体からランダムに選択された遺伝子型(約80000個の常染色体SNP)が含まれている.qt.phe:数量性状は単独の代替代表型ファイルにあります.pop.pheファイルには仮想表型が含まれており、中国人の個体コードは1、日本人コードは2である.集団間の違いを調べるために使用します.
mapとped形式のファイルをフォーマット変換する
--fileオプションのmapおよびpedフォーマットファイルの接頭辞では、対応するbed,fam,binフォーマットのファイルが生成されます.plink --file hapmap1
このコマンドが完了すると、次の情報も表示され、この変異ファイルには83534と89個の個体が含まれていることが簡単にわかります.
PLINK v1.90b6.3 64-bit (17 Jul 2018)
Options in effect:
--file hapmap1
Hostname: zeus-1
Working directory: /scratch/pawsey0149/hhu/test/plink/new
Start time: Sun Jul 28 15:26:17 2019
Random number seed: 1564298777
257868 MB RAM detected; reserving 128934 MB for main workspace.
Scanning .ped file... done.
Performing single-pass .bed write (83534 variants, 89 people).
--file: plink.bed + plink.bim + plink.fam written.
End time: Sun Jul 28 15:26:18 2019
生成されたファイルをバイナリPEDファイルとしてさらに圧縮することができます.
plink --file hapmap1 --make-bed --out hapmap1
たとえば、高周波遺伝子分型(少なくとも95%以上の完全性)を持つ個体のみを含む新しいファイルを作成する場合は、次のように実行できます.
plink --file hapmap1 --make-bed --mind 0.05 --out highgeno
SNP変異のデータをまとめる
紛失率の確認
準備が終わった後、SNPのデータを統計し始めました.まず、損失率
missing rateを見てみましょう.plink --bfile hapmap1 --missing --out miss_stat
ここでは、ファイルをバイナリbed形式に変換したので、
--bfileオプションを使用します.2つのファイルが生成されます.
-rw-r----- 1 hhu hhu 3.6M Jul 28 15:41 miss_stat.lmiss
-rw-r----- 1 hhu hhu 4.6K Jul 28 15:41 miss_stat.imiss
1つのファイルは
miss_stat.lmissであり、ここでは変異missing rate情報が含まれており、もう1つのファイルはmiss_stat.imissであり、提案されたmissing rate情報が含まれている.この2つのファイルを簡単に配布して、よりよく理解させます.more miss_stat.lmiss
CHR SNP N_MISS F_MISS
1 rs6681049 0 0
1 rs4074137 0 0
1 rs7540009 0 0
1 rs1891905 0 0
1 rs9729550 0 0
1 rs3813196 0 0
1 rs6704013 2 0.0224719
1 rs9439440 2 0.0224719
ここではSNP毎に失われた個体数(N＿MISS)と失われた個体の割合(F＿MISS)を示した.別のファイルを見てみましょう.
more miss_stat.imiss
FID IID MISS_PHENO N_MISS F_MISS
HCB181 1 N 671 0.00803266
HCB182 1 N 1156 0.0138387
HCB183 1 N 498 0.00596164
HCB184 1 N 412 0.00493212
HCB185 1 N 329 0.00393852
HCB186 1 N 1233 0.0147605
HCB187 1 N 258 0.00308856
ここでは個々について示したが、失われたSNPの合計数(N＿MISS)と、各個体における失われたSNPの合計数(F＿MISS)の割合(F＿MISS)は、一般的に、この個体における失われたSNPの合計数(N＿MISS)の割合が高い場合、後続の分析の前に、その個体を除去することを選択する.ここでは,このテストデータの品質は悪くなく,個体のSNP損失の割合は非常に低い.
また、ある染色体の損失率にのみ興味がある場合は、対応する染色体の損失率情報を出力するように指定することもできます.
plink --bfile hapmap1 --chr 1 --out res1 --missing
等位遺伝子の頻度を検査する
ここでfreq_が生成されますstat.frqファイルは、SNP毎の二次等位遺伝子周波数(MAF)を含む.
plink --bfile hapmap1 --freq --out freq_stat
more freq_stat.frq
CHR SNP A1 A2 MAF NCHROBS
1 rs6681049 1 2 0.2135 178
1 rs4074137 1 2 0.07865 178
1 rs7540009 0 2 0 178
1 rs1891905 1 2 0.4045 178
1 rs9729550 1 2 0.1292 178
1 rs3813196 1 2 0.02809 178
ここでも指定されたSNPに対して、その等位遺伝子の頻度を調べることができます
plink --bfile hapmap1 --snp rs1891905 --freq --out snp1_frq_stat
ここでMAFといえば、簡単にMAFとは何かを拡張しますが、なぜMAFでフィルタリングするのですか?
MAFは、副等位遺伝子の頻度、すなわち、頻度の低い等位遺伝子の頻度である.MAFは様々な理由で濾過することができる.そのうちの1つは,低い被覆量または深さを持つ配列では,非常に低いMAFを持つ試料が多くの騒音を発生し,偽陽性の下流分析を引き起こす可能性があることである.極端な例を挙げると、1つの個体にのみ存在する等位遺伝子(非常に低い等位遺伝子周波数)は、集団構造推定に有用で信頼できる情報を追加しない.データが全ゲノム関連研究(GWAS)に使用される場合、非常にまれな等位遺伝子を有意に記述するために非常に強い統計的能力が必要であるため、通常、比較的高いMAF閾値が用いられる.一般的によく使われるバルブ値は0.05と0.01です.しかし、最適なMAFフィルタしきい値はあなたのデータ型に依存し、計画された分析、データ品質、サンプルサイズも統計分析で計算できます.
ここでは、バルブ値0.05を使用してフィルタリングする例を示します.
plink --bfile hapmap1 -maf 0.05 -out test --make-bed
83534 variants loaded from .bim file.
89 people (89 males, 0 females) loaded from .fam.
89 phenotype values loaded from .fam.
Using 1 thread (no multithreaded calculations invoked).
Before main variant filters, 89 founders and 0 nonfounders present.
Calculating allele frequencies... done.
Total genotyping rate is 0.99441.
24114 variants removed due to minor allele threshold(s)
(--maf/--max-maf/--mac/--max-mac).
59420 variants and 89 people pass filters and QC.
Among remaining phenotypes, 44 are cases and 45 are controls.
--make-bed to test.bed + test.bim + test.fam ... done.
濾過後、24114変異部位が濾過され、59420変異部位が保持されることがわかる.
MAFに加えて、
--mindおよび--genoも一般的なフィルタパラメータであると考えられる.遺伝子型喪失率が10%より大きい個体サンプルを除去する:
plink --file hapmap1 --recode --mind 0.10 --out hapmap1_complete
95%を超える個体が持っていない変異部位を除去する.
plink --file hapmap1 --make-bed --geno 0.05
もちろん、すべてのフィルタリングは同時に実行できます.
plink --file hapmap1 --make-bed --mind 0.10 --maf 0.05 --geno 0.05 --out hapmap1_clean
最後に簡単に説明すると、PlinkはLDの計算にも使用できます.LDの計算についてもっと詳しく議論します.
plink -bfile hapmap1 -r2 -ld-window-kb 1 -ld-window 1000 -ld-window-r2 0 -out test
もちろん、いくつかのRパッケージ、例えば
adegenetやstructure図を描くstructureは、特殊な入力ファイルフォーマット(出力可能なフォーマットが多く、ここではplinkに戻るmanual検索を詳細にサポートする1つだけを列挙しています)が必要な場合があります.この場合、Plinkは同時に手伝うことができます.plink -bfile hapmap1 -aec -out man -recode A
ここまで私が日常的に使っているPlink機能を紹介しましたが、これらの機能は実はPlinkの氷山の一角にすぎず、その強さはmanualをよく読んでこそ深く体得することができます.