SSブログ

日本人と中国人の違いを「マンハッタンプロット」で調べる [ゲノム解析で古代史]

古代人や縄文人のゲノム解析ですが、なぜか「マンハッタンプロット」を使った分析は見たことがありません。

スクリーンショット 2022-05-05 101534.png

上の図にある説明では、「疾患との関連が高い」とありますが、要するに「異なるグループの遺伝子(SNP)の違い」が偶然かどうかの確率を示しているわけです。

誰でもすぐに思いつきそうなことですが、本当に誰もやってないのかな?

ここでは、テストケースとして、前回調べた1000 Genome Projectのデータを使った日本人と中国人の違いをマンハッタンプロットを使って見てみることにします。
第9染色体.png
前回の分析結果

下の図は、ABO血液型がある第9染色体後半部分のマンハッタンプロットを描いてみたものです。
赤い↓がABO血液型の場所なのですが、残念なことに意外と差が出てないですね。

manhattanplot3.png

点が上に行くほど差があることを示しているのですが、あまりパッとしません。

なお、水平の赤い線以上が統計的に有意なことを示しています。

もう一つ残念なのは、ノイズが多いのため、このままでは使いものにならないことです。

どうしようかな…。

--☆----☆----☆--

なお、具体的なやり方は次のとおりです。

★PLINKについて★

【1】PEDファイルにPhenotypeを追加する

前回までにできた、マンハッタンプロットをしたいPEDファイル(テキストファイルです)にPhenotypeを追加します。

PLINKの使い方 1 〜フォーマット整形〜
http://bio-and-info.blogspot.com/2013/05/plink1.html

2. デフォルトのフォーマット
(1) PED FORMAT:サンプルとジェノタイプ情報
ここの、列#6 Phenotypeを、EmEditorで追加(変更)します。
(デフォルトは「-9」のようです…)
・例 日本人=1、外国人=2

【2】関連分析

plink --file inputfile --assoc --allow-no-sex --out outputfile
→outputfile.assocに結果が出力されます。
(--allow-no-sexは性別を無視するオプションで、元データに性別が入ってないため、指定しないとエラーになります)

《参考》
・PLINKでGWAS(basicなcase/control study)をするコマンド
 https://jojoshin.hatenablog.com/entry/2016/09/22/094319
 plink --bfile データ名 --assoc --out アウトプットするデータ名
 →上のファイルはバイナリではないので、「--bfile」を「--file」に変更する

★Rについて★

カラフルなマンハッタンプロットを描く
https://staffblog.amelieff.jp/entry/2018/05/15/173432

《Rのスクリプトの例》

rm(list=ls())
setwd("C:\\データがあるフォルダ名")
#必要なパッケージをロードする。
library("qqman") #ライブラリ読み込み
library("data.table")
dat.raw <- fread("データファイル名.assoc", stringsAsFactors=F,header=T)
dat.ylim <- c(0, max(10, max(-log10(dat.raw$P) * 1.04, na.rm=TRUE)))
manhattan(dat.raw, ylim=dat.ylim, genomewideline = -log10(1e-07))

コメント(0) 

コメント 0

コメントを書く

お名前:[必須]
URL:
コメント:
画像認証:
下の画像に表示されている文字を入力してください。

※ブログオーナーが承認したコメントのみ表示されます。

Facebook コメント