ryamadaの遺伝学・遺伝統計学メモ このページをアンテナに追加 RSSフィード

数学・コンピュータ関連の姉妹ブログ『ryamadaのコンピュータ・数学メモ』
京都大学大学院医学研究科ゲノム医学センター統計遺伝学分野のWiki
講義・スライド
医学生物学と数学とプログラミングの三重学習を狙う学習ツール
駆け足で読む○○シリーズ
ぱらぱらめくるシリーズ
カシオの計算機
オンライン整数列大辞典

2015-07-25 Graph measures and indices

[]メモ

2014-12-24 電子書籍で振り返る1年

ryamada222014-12-24

[][][]電子書籍で振り返る1年

  • 今年は、資料を電子ファイル化して保全することに取り組みました。
  • 基本的にはRmd形式のテキストファイルで作り、そこからhtml化し、epub化して、kindleに保管をしてもらう、という形です。
  • Rmdからhtml/epub化することを覚えれば無料で入手できます。
  • kindle形式はTeXの数式との相性が悪いので数式がうまくいっていない場合があり、Rが描く図も場合によってはうまく収まっていませんが、Rmd->html/epubができない方はkindleもご利用ください。
  • 基本的には、1米ドル。少し中身に自信があるもの・分量が多いものは2米ドルくらいです。
  • キンドルでの一覧

f:id:ryamada22:20141223181348p:image

2010-10-14 量的疾患の解析手法メモ

ryamada222010-10-14

[][][][][]量的疾患の解析手法のメモ

  • 複合遺伝性(疾患)
    • 1ローカスによるメンデルの法則によって説明できず、複数の(通常は多数の)ローカスの影響によって決まる遺伝形質(疾患)
  • 遺伝形質
    • 量的形質
    • 質的形質(疾病など)
  • 多遺伝子性(Polygenic)
  • 身長の遺伝子論文
    • Hundreds of variants clustered in genomic loci and biological pathways affect human height
    • たくさんのローカスが見つかった
    • ローカスは相互に独立ではなく、パスウェイでまとまる傾向があった
    • マッピングに用いたマーカーの近傍に「らしい遺伝子」があった
    • 1つのローカスに複数の関連起源(allelic heterogeneity)があった
    • 遺伝子機能に影響を与えると想定される多型(アミノ酸置換・発現への影響)が認められることが多い
    • 個々のローカスが説明する遺伝因子の量(phenotypic variation全体のうちのどれくらいを説明するか、遺伝率のうちのどれくらいを説明するか)は小さい
    • 今回、認めた多数のローカスのすべてを合わせても、phenotypic varianceの1割未満しか説明しない
    • 今回、認めたローカスと同程度のローカスが見つけ損ねられていると考えられるが、それがほとんど見つかったとしても、15%程度しか説明しない(→もっと弱いローカスがさらにたくさん集まって説明するかもしれない 同程度のローカスのallelic heterogeneityのすべてがみつかっていないので、その分がみつかると、説明量が増える可能性はある)
  • その解析手法をぱらぱら
  • Primary genome-wide association meta-analysis (stage 1)
  • In silico follow-up (stage 2)
  • Meta-analysis of stages 1 and 2 (stage 1 と stage 2をメタで合わせる)
  • Validation by genotyping and population stratification analyses
  • Percentage variation explained and number of loci
  • Gene by Gene interaction, dominant and recessive analyses
  • Conditional analyses
  • Functional variant analyses
  • Biological enrichment analyses
  • Primary genome-wide association meta-analysis (stage 1)
    • メタで統合された個々のスタディに関して
      • サンプルQC
        • Call rate
        • その他スタディ別のクライテリア(男女の不一致、メンデリアンエラー、フェノタイプ情報の欠失、PCA outliers など)
      • マーカーQC
        • Call rate
        • HWE-p
        • MAF
      • Imputation
        • 方法:BIMBAM,MACH,IMPUTE
        • MAF cutoff
        • Quality control score(r2-hat,proper-info)
    • Association
      • sex-standardized height
      • additive inheritance model
      • Adjustments
        • age
        • その他、スタディごとの因子(PCA-based componentsなど)
      • ジェノタイプとフェノタイプの関係づけ
        • Unrelated individuals: linear regression
        • Related individuals: Variance component, linear mixed 他(こちらも参照)
      • Linear regression もVariance componentも、線形解析と言ってしまえばそれまでで、同じことですが・・・
    • メタアナリシス
      • MAFxN>3(アレルの本数が3以上)フィルタ(回帰係数が出せないから)
      • Inverse-variance fixed effects meta-analysis(一番シンプルなやり方とも言える)
      • Sample-size weighted Z-score-based fixed effect meta-analysisは補助
      • ソフトウェア:METAL
      • メタ後にGCラムダが大きくて難渋している・・・結局、メタ後に再度GC補正している
  • In silico follow-up (stage 2)
    • Stage1で上がったマーカーとその近傍でLDにないマーカー
  • Meta-analysis of stages 1 and 2 (stage 1 と stage 2をメタで合わせる)
    • 男女別で解析
    • 年齢別に2群に分けて解析
    • ソフトウェア:METAL
  • Validation by genotyping and population stratification analyses
    • 関連陽性マーカーから抽出の上、別方法でリタイピング
    • 異民族からの低身長・高身長のextreme-samplesでバリデーションテスト(Cochran-Mantel-Haenszel)
    • パワー計算:Genetic Power Calculator
    • Familial sampleでバリデーション
      • PLINKのQFAM
      • 核家族にばらしてFBAT
      • Imputationについてはbest-guessのジェノタイプを使用
      • メタするときは、sample-size weighted (ただし、effective sample sizeを用いる。Effective sample sizeとは、素のサンプルサイズをGCラムダで除したもの)
    • Ancestry Informative Marker(AIM)は、民族分離力の強いマーカー。欧州系集団を亜集団に分ける力の強い基地マーカー(HLAとlactase lociを含む)と、身長関連マーカーとの関連の強さをLDインデックスで評価
    • EIGENSTRATでも評価(南北軸は身長と相関あり)
    • 身長関連マーカーと南北軸との関係は、それ以外のマーカーと南北軸との関係と比べて、有意に強いということはなかった
    • Fstも身長関連マーカーがその他のマーカーと比べて大きいということはなかった
    • ここで、「その他のマーカー」は、アレル頻度でそろえた、マーカーセットをランダムに発生させて分布をとった
  • Percentage variation explained and number of loci
    • Phenotypic varianceのうちのどれくらいを説明するかの推定
      • オリジナル論文はNature. 2009 Aug 6;460(7256):748-52. Epub 2009 Jul 1.Common polygenic variation contributes to risk of schizophrenia and bipolar disorder.こちらのサプルの"16)"
      • Simulation study to estimate the variance explained by common variants
        • Phenotypic varianceのどの程度を説明するローカスがいくつあるかのモデルを立てる
        • ローカスが関連検定でみつかるかどうかは、ローカスの強さ(phenotypic varianceのどの程度を説明するか)と、検定閾値で決まる
        • 検定閾値を動かすと、「検出されたローカスのセット」によって説明されるphenotypic varianceの割合が変化するが、検定閾値がどのあたりで割合の変化が大きいかは、「説明ローカスの強さと数の分布」によって変わってくる
        • これを利用して、「当てはまりのよいモデル」を推定する
        • 実際のスタディと同様のシミュレーションデータを用いて行う推定なので、それを作る必要がある
          • だいたい、HWEとLEと考えられるマーカーを選んで、そこからHWE・LEを仮定してサンプリングする
          • LEの仮定は「だいたい」だけれど、実際にはある程度のLDが残っているので、それは、関連検定統計量の分散インフレーションの程度をシミュレーションデータから算出して、「effective number of markers in LE(LEだとしたときの実効マーカー数)」を出してそれを使う(このあたりは、『方便』)
        • ローカスの強さの分布をどうとるかのモデル立て(これも『方便』)
          • GRR (Genotypic relative risk)で考えれば:
            • GRRとアレル頻度に関連を持たせない
            • GRRとアレル頻度に関連を持たせる(説明分散を等しくする)
            • GRRの分布を仮定する
        • マーカーと「真の関連多型」との間のLDの強さを考慮するかしないか
      • と、このように、ごちゃごちゃしているので、推定値(たとえば、身長論文では、小数点以下1桁のパーセンテージで示されている)は、目安
  • SNPの分散と、身長の分散
    • ある2アレル多型があって、そのアレル頻度がp,1-pとする。HWE仮定のもとで、3ジェノタイプの頻度はp^2,2p(1-p),(1-p)^2。この頻度分布だと、3ジェノタイプに0,1,2という値を与えれば、その平均は2(1-p)、分散は2p(1-p)、-G*kの分散は2p(1-p)xk^2
N<-100000
p<-0.3
G<-sample(c(0,1,2),N,prob=c(p^2,2*p*(1-p),(1-p)^2),replace=TRUE)
mean(G)
var(G)
2*p*(1-p)
k<-4
var(G*k)
var(G*k)/(2*p*(1-p))
    • 今、身長がこの多型の3ジェノタイプに対応して、分散10の力があるとすれば、2p(1-p)*k^2=10なのでk=¥sqrt{¥frac{10}{2p(1-p)}}
k<-sqrt(10/(2*p*(1-p)))
var(G*k)
    • また、この多型以外はすべて環境因子できまり、それも分散10の力があるとすれば、その和Hの分散はGの分散(10)と環境要因の分散(10)の和になっていて、GはHの分散(phenotypic variance)の半分を説明している
E<-rnorm(N,mean=170-mean(G*k),sd=sqrt(10))
# 身長の平均が170になるようにEの平均を170-mean(G*k)にしている
H<-G*k+E
hist(H)
var(H)
var(E)
var(G*k)
  • さて、実際には、フェノタイプ(身長)の個人データと、ジェノタイプの個人データが得られる。そこから、どうやって、ジェノタイプがフェノタイプの分散のどれほどを説明するかを計算するか、である
    • 複数のマーカーのジェノタイプがある
    • それぞれのマーカーに重みづけをする(重みづけの方法は上述したように、GRRとして得る)
    • 重みづけが出れば、¥sum_{i=1}^{N_m} w_i G_iで、ある個人の値が決まる。ただし、N_mはマーカー数、w_iはマーカーのアレル1本あたりの重みで、G_iはこの個人が持つアレルコピー数。式からわかるように、個々のマーカーはアレル本数について線形、マーカーの合算も線形としている(線形でないモデルにしたいなら、この計算式を変えればよい)
    • 今、身長をH、ジェノタイプの個人の値をGとして、その関係を線形回帰モデルに当てはめ、adjusuted R^2を出し、それが、ジェノタイプで説明された分散
Gk<-G*k
lmout<-lm(H~Gk)
summary(lmout)
summary(lmout)$adj
    • Prediction of number of susceptibility loci
      • 元論文は Nat Genet. 2010 Jul;42(7):570-5. Epub 2010 Jun 20. Estimation of effect size distribution from genome-wide association studies and implications for future discoveries. こちら
      • Effect size
        •  e = ¥beta^2 ¥times 2f(1-f)
          • 2f(1-f)はあるローカスがHWEにあるとして、3ジェノタイプの寄与が0,1,2とadditiveなモデルにあり、アレル1コピー分の寄与を1単位としたときの分散
          • ¥betaは当該アレル1コピー分の寄与につけた重みで、Effect size eがこのローカスのHWE、additiveモデル仮定で寄与する分散量
      • 実施スタディからは、「検出されたEffect size」の分布が出る
      • 真のeffect sizeが検出されるパワーは、スタディデザインによる
      • 真のeffect size別のローカス数を推定する
      • 推定するにあたり、effect sizeについて分布を仮定して(指数分布やワイブル分布)、その分布を定めるパラメタを推定するのがパラメトリック手法。観察データの分布からカーネル平滑化して分布をとることもできル(ノンパラメトリック手法)
  • Gene by Gene interaction, dominant and recessive analyses
    • 単一ローカスのモデル解析
      • additive, dominant, recessiveそれぞれのモデルについて係数計算をしたうえで、additiveからの逸脱についても検定している
    • Genome-wide joint effect analysis
      • メソッドの記載がちょっとわかりにくいけれど、個別マーカー解析で閾値未満のp値をとったマーカーと、すべてのマーカーとの間のペアワイズの効果判定をした
      • 効果判定は、2通りのモデルで実施
        • 2マーカーのそれぞれにadditiveなモデルを入れ、2マーカーのアレルの関係の項を入れ、2マーカーアレルの関係項について検定
        • 2マーカーが作る9ジェノタイプの分布についての検定(自由度8)
      • 個別ジェノタイプのデータが得られた一部のサンプルについての実施となっている
      • ソフトウェアはPLINK
  • Conditional analyses
    • 1次・2次で出てきたマーカー180+について調整(回帰の式に180+の項を入れたということと思います)したうえで、再度、個々のスタディで関連を調べ、それをメタでまとめた
  • Functional variant analyses
    • Transformed lymphoblastoid cell linesのeQTL
    • cisとtrans
      • cisは遺伝子までの距離が1Mbのもの
    • Genotypeと発現量との関係を検定(回帰)
    • Non-genetic effectsで調整・・・と書いてあるけれど、具体的には何で調整したのか不明
    • パーミュテーション-basd FDR(サンプルの家族関係、マーカーの関係(LD)を保ったシャッフリング(単純なFDRによる期待関連数より高件数で関連同定)
    • Primary OsteoblastsのeQTL
      • 少し工夫した手法でAllelic Expression(AE)法
      • Nat Genet. 2009 Nov;41(11):1216-22. Epub 2009 Oct 18. Global patterns of cis variation in human cells revealed by high-density allelic expression analysis. (こちら)
      • トランスクリプト上にあるSNPを用いる
      • genomic DNAとmRNA由来のcDNAの両方をサンプルとして用いる
      • genomic DNAは2アレルの比が1:1であること、cDNAの方はアレルごとにトランスクリプト量の比が1:1からずれる(こともある)ことを利用する
      • 今、興味があるのは、個々のSNPではなくて、遺伝子=トランスクリプト全体なので、1遺伝子=1トランスクリプトあたり、複数(3個以上)のSNPでもって、「遺伝子=トランスクリプト」を単位としてアレル別cDNA量非を推定する手法
    • Liver, omentum, and subcutaneous fat eQTL
      • 調整因子:age race gender and surgery year
      • Genotype~Expression量
    • Non-synonymous enrichment analysis
      • 検出された180マーカーについて、そのLD的近傍のSNPのセットを定め、そのセットの中に「機能性SNP(nonsense SNP,missense SNP)」が富んでいることを検定する
      • 比較対象は、調整したうえで「ランダムに選んだ」同数のSNPのセットを1000通り作ったもの
      • 調整因子は、アレル頻度、近傍遺伝子数、遺伝子への物理的近さ
    • 他の遺伝形質との関連もテスト
  • Biological enrichment analyses
    • OMIM
      • OMIMを使って、「身長」と関係ありと目される遺伝子をリスト化
      • 検出SNPの帰属LD領域を確定
      • OMIM-basedな「身長遺伝子」がSNP帰属領域に集中気味かどうかを検定
    • Text-mining GRAIL
      • GRAILの論文はこちら
      • 検出SNPから、対象領域をLD的に定める
      • そこにある遺伝子をリストにする
      • そのようにして挙がった遺伝子をPubMedのアブストラクトについて探索したうえで、遺伝子間の距離に変える
      • アブストラクトに表れる単語の共通性が遺伝子間距離の評価基準となる
      • 検出SNP領域が「身長」に関してGLAIL的に(PubMedアブストラクト的に)有意に集中していることを示す
      • ソフトウェア:GRAIL
    • Pathway
      • Gene set enrichment analysis(GSEA)
      • 遺伝子ごとの関連の強さを、次の要素を考慮して評価
        • 遺伝子サイズ、その遺伝子に配されたマーカー数、そのLDの具合
        • 実際には、パーミュテーションをまわして、遺伝子ごとの「トップ統計量」(トップP)の分布をとり、それに応じて、補正する。こうすることで、「配したマーカー」についての補正がなされ、遺伝子に関する統計量(p値)となる
        • 遺伝子の範囲は上下流、それぞれ110kb,40kb
      • 遺伝子ごとの関連の強さで見たときの上位5%の遺伝子が、集中しているパスウェイが何かを調べる
        • 次のデータベースを使ってパスウェイの帰属遺伝子を決める
        • 「ランダムな遺伝子セット」に比べて、「パスウェイでくくられた遺伝子セット」が、上位5%遺伝子をenrichに持つことの確認
        • 評価はパーミュテーションによるP、FDRによるP

2009-12-21 ぱらぱらめくるnature reviews genetics

[][][][]Fst再考

2009-08-02 Fst,F,heterozygosity,IBD

[]集団構造化・Inbreeding・IBDの指標の意味するところ

集団の構造化をFstで表したり、Inbreeding coefficiency Fがヘテロ接合体の比率に影響を与えたりする。FstはHeterozygosityを階層的に定め、それらの関係を示す一連のF値の一つである。

これらは、どういう指標であって、どういう関係にあるのだろう。

構造化とinbreedingとは、ともにrandom matingからずれている状態を表している。

random matingからのずれは、アレルの頻度に基づいて、アレルが相互に独立であると仮定したときのdiplotypeの頻度と、実際の頻度の違いとを生じるので、構造化もinbreedingもそれを利用して定量化することができる。

ホモが増えることに着目してもよし、ヘテロが減ることに着目してもよいが、ヘテロが減ることに着目することが多い。

期待されるヘテロの程度と実際のヘテロの程度が分かったら、それの比をとるもよし、期待されるヘテロと実際のヘテロが一致する状態と、ヘテロが少ない極限の状態との2状態を基準にして、その2基準点に照らして量表現することも可能である。基準2状態に照らして量表現するときに、線形式を使って0から1までの量で表すことは、単純でわかりやすいのという理由から指標化の際に用いられる。

この期待状態と極端状態についての線形数値化の指標がFstを含む、Fの指標たちである。いろいろなFがあるのは、対象を何にして、期待状態の定義を何にするかに依存するためである。

この流れのなかで階層的構造化指標を考える。

今、何でもよいけれど、相互に何かしら関係のある状態がS=¥{s_1,s_2,...,s_n¥}あったとする。それぞれの状態について、ある0でない量V=¥{v_1,v_2,...,v_n¥}が定義出来るとする。

すると、¥prod_{i=1}^{n-1} ¥frac{v_i}{v_{i+1}}=¥frac{v_1}{v_n}

このような状態セットと対応する量の値のセットがあるときに、状態に順序があると、状態の階層構造がVによって、『どこの階層を切り出しても、切り出す階層の厚みを増減させても』関係式が維持されるので、便利です。

Heritabilityについてこれを用いて、しかもF_{i,j} = ¥frac{H_j-H_i}{H_j}という定義を用いて、(1-F_{IS})(1-F_{ST})=1-F_{IT}なる関係であることが、いわゆるF_{ST}に関係する3つのFの背景である。つまり、F_{ST}についてはI,S,Tの3状態の階層構造を特に取り出していることになるが、n=3に限ル必要は本日的にはない。

F_{IS}はinbreeding coefficient Fであるが、これは、同祖でないアレルのみからなる祖先集団から近親交配を含む交配を経て得られた集団があったときに、Fの割合で同祖のとなり1-Fの割合で非同祖となるようなメイティング過程だったと仮定したときに得られる、ディプロタイプ分布がFを用いて表現できる、という点でFはIBDとつながってくる。

さらに進んで、近親交配を含む家系図があったときに、そのinbreeding coefficientの計算は家系図が作るループについて、F=¥sum_{L=¥{l_i¥}: loops} (¥frac{1}{2})^{distance(l_i)}¥times f( l_i)、ただし、f(l_i)はループl_iにおける最も祖先側の個体のinbreeding coefficint。