PLINK
以下のコード中では、”–“（ハイフン\(2\)つ）と”-“（ハイフン\(1\)つ）の違いには十分注意。
データのQuality controlのプロトコル
こちらの記事を参考に。
https://www.huble.org/wp-content/uploads/2020/06/HubLE-I-Methods-GWAS-11.pdf
必要なもの
- Genome Wide Association Study (GWAS)に必要なデータ
- コンピュータ
- Unix/Linux
- PLINK
- Human reference panel (1000 Genomes project)
1. ファイルの準備
PLINKの.bed, .bim, .famなどのbinaryファイルを準備する。最初のデータ形式が.pedや.gen/.bgenの場合は以下のコマンドを用いる。
plink --file mydata --make-bed --out mydata2こうすることで”mydata.ped/.map”という\(2\)つのファイルから”mydata2.bed/.bim/.fam”という\(3\)つのファイルが生成される。.bed/.bim/.famというbinraryファイルを読み込むときには、
plink --bfile mydata2...とする。
2. SNPのフィルタリング
2.1 Call rate
集団の中でのCall rateが低い遺伝子を以下のコマンドで除く。Call rateというのは、研究の対象となる集団で、あるSNPについて正確な配列が判明する人の率。なので、あるSNPについて100人中99人のSNPが正確に判明している場合、そのcall rateは99%となる。
plink --bfile mydata --geno 0.01 --make-bed --out mydata_geno2.2 Ninor allele frequency (MAF)
SNPの中でも稀なalleleものはどんどん集団から失われていく。そのようなalleleを解析の対象とすると、時間ばかりかかって研究成果が出ないので、除いておく。コマンドは以下。
plink --bfile mydata --maf 0.01 --make-bed -out mydata_maf2.3 Hardy-Weinberg equilibrium (HWE)
HWE検定を通過しないマーカーは除外される。この値であるが、色々な流儀があるらしいが、一般的に\(1.0\times 10^{-6}\)が用いられることが多いみたい。
plink --bfile mydata --hwe 0.000001 --make-bed --out mydata_hwe3. 個人やサンプルレベルでのQC
3.1 Call rate
同じ個人の中で、あまりにも多くのSNPが決定されていない場合、取り除かなければいけない。
plink --bfile myata --mind 0.01 --make-bed --out mydata_missingどの個人が失われたのかを要約させるには、
plink --bfile mydata --missing --out mydata_missingとする。
3.2 Heterozygosity
ヘテロ接合の観察数(O)と期待数(E)を評価することで、ヘテロ接合の多い個体や少ない個体を選別する。
plink --bfile mydata --het --out mydata_het3.3 性情報の不一致
X染色体があるときは、性染色体と性別の情報から、それが一致しない個人を除く。
plink --bfile mydata --check-sex --out mydata_sex3.4 集団の層別化
ある特定の祖先や、（遺伝的に）親しい集団というのが研究参加者に含まれていると、それは当然交絡因子になる。なぜなら、遺伝子にもアウトカム（疾患の発症など）にも当然上流から影響を与えうるからだ（同じような生活をしていて、同じような疾患発症リスクを有している可能性などを考慮できるから）。こういう場合は主成分分析を行う。また、pairwise IBS (identity by state)を行うのは、以下のコマンドとなる。
plink --bfile mydata --genome --out mydata_IBS1000 Genomesが提供するリファレンスpanelを用いた、集団間の差異を調べるコマンドは以下である。
plink --bfile mydata --read-genome plink.genome --ibs-test与えられたサンプル（参加者）のグループ層別化は、以下でプロットできる。
plink --file mydata --read-genome plink.genome --cluster --mds-plot43.5 近親度や隠れた家族関係
あるペアの個々人について、どの程度遺伝的な関係が強いのかを調査できる。IBD (identity by descent)を調べるのは、上記の”–genome”コマンドを用いる。近親交配によるSNP間の強い相関を避けるため、LD (Linkage disequilibrium)刈り取りを行い、\(r^2\)の閾値を\(0.2\)として、互いに相関がない・あるいは弱いSNPのサブセットを以下のように選択する。
plink --bfile mydata --indep-pairwise 50 5 0.2 --out mydata_IBSいくつかのTIPs
マイナーアレル頻度(MAF)の閾値は、研究の規模や変化の影響の大きさによって異なる。 大規模な研究で、\(N = 100,000\)などの場合はMAFの閾値を\(0.01\)とする傾向があるが、小規模の研究ではMAFの閾値を\(0.05\)とする。 小規模な集団に存在する稀なSNPsは、ジェノタイピングエラーを起こしやすい。PLINKでは、MAF=0.01をデフォルト値としている。
HWE平衡からの逸脱は、一般にGWASにおけるジェノタイピングエラーを示すが、集団の下部構造から生じることもある。または進化的選択に起因する場合もある。 症例対照研究において、HWEフィルタリングは対照群でのみ行われる。なぜなら、症例群ではHWEルールに反することが多いからである。
入力された血統ファイルの性別とX染色体解析で得られた性別の間に矛盾がある場合、サンプルの取り違えの可能性がある。 X染色体ホモ接合性推定値(F)の閾値は、男性で\(>0.8\) 、女性で\(<0.2\)とするのが一般的である。また、染色体異常（ターナー症候群やクラインフェルター症候群など）により、性別の不一致が生じることもある。性別の不一致の原因を視覚的に調べるツールとしては、GenomeStudioソフトウェア（イルミナ社）がある。
対立遺伝子頻度が集団間で異なることがあるため、層別化は偽の関連性をもたらし、真の関連性を覆い隠してしまう可能性がある。 Zスコアは標本の平均と分散に基づいて計算され、極端な値を持つ標本は外れ値とみなされ、削除されるべきである。
GWASでは、すべての研究参加者が非血縁であることが前提となっている。 そのため、遺伝的関係の強い者をGWASに含めると、SNPの効果推定に偏りが生じる可能性がある。IBD親族係数(PI_HAT)\(>0.1\)は、ジェノタイピング前にサンプルDNAのコンタミネーションがない限り、テスト対象のサンプル間の強い関連性または重複を示唆する。この場合、ペアごとに1つのサンプルが解析から除外される（家族構成を考慮・調整できるソフトウェアが使用されていない場合）。除外するサンプルは、通常、遺伝子型のCall rateが低いもの、またはサンプル数の損失や表現型の利用可能性を最小限にするものである。
上記をまとめたフローチャートはこちら。
コメント