NGSデータ解析まとめ

サカナ研究者の手探りNGS解析(おもに進化生物学)

MiSeqデータのMHC領域へのマッピング(15):RによるVCF fileの操作 (1): SNPごとの対立遺伝子頻度(allele frequency)および遺伝子型頻度(genotype frequency)の計算

2015-06-22

ここでは、Rを使ってVCF fileに含まれる情報を整理する方法を考えます。具体的には、

  • VCF fileに含まれる各sampleのgenotype情報から、SNPごとの対立遺伝子頻度(allele frequency)および遺伝子型頻度(genotype frequency)を求める
  • VCF fileから必要な情報だけを抽出し、SNPsについての新しい表を作成する
  • SNP positionとallele frequencyとの間の関係(散布図)を作成する

以上3つのことを行います。

VCF fileのRへの読み込み、および基本的な操作については、このエントリを参考して下さい。

(復習)VCF fileのRへの読み込み

読み込むVCF fileですが、前回Variant Effect PredictorアノテーションをしたVCF file(MHC_test4.var.annotated.vcf)を読み込みます。読み込む方法はこのエントリを参考して下さい。

もし"MHC_test4.var.annotated.vcf"がない場合はこちらへ。

1. SNPごとの対立遺伝子頻度(allele frequency)および遺伝子型頻度(genotype frequency)の計算

VCF fileでは、sampleごとのgenotype情報は、Genotype fieldの列に書かれています。そこでは、

1/1: ALT alleleのホモ接合体
0/1: REFとALT allelesのヘテロ接合体
0/0: REF alleleのホモ接合体

を示しています(このエントリも参照)。これはこれでOKなのですが、この表示のままでは、対立遺伝頻度や遺伝子型頻度の計算には不便なので、これを以下のように変換します。

1/1 -> 2
0/1 -> 1
0/0 -> 0

また、genotype fieldの各列における、各alleleの統計的な評価(参照ページ)の数字やコンマも邪魔なので除きます。

(1) Genotyping fieldの変換(正規表現による)

まず、Rで以下のコマンドを入力します。

# V50 genotype1
genotype1 <- sub("1/1:[0-9]+,[0-9]+,[0-9]+:[0-9]+","2", MHC_test4_annotated[,50])
genotype1 <- sub("0/1:[0-9]+,[0-9]+,[0-9]+:[0-9]+","1", genotype1)
genotype1 <- sub("0/0:[0-9]+,[0-9]+,[0-9]+:[0-9]+","0", genotype1)
genotype1 <- as.numeric(genotype1)

ここでは、Rのコマンド"sub"と正規表現を使って、genotype fieldの列(50列目)の値を変換したベクトル(複数の数値からなる変数)"genotype1"を作成しています。最後のコマンドの"as.numeric"は、変換した値のデータ型を「文字列」から「数値」に変換する命令です。コマンドは1行ずつ入力しても良いですし、Rの「ファイル」-> 「新規文書」でエディタを開き、そこにペーストしてからコマンド全体を選択し、「command」を押しながら「return」を押すと、すべての行を上から順に実行します。

参考URL:
R.4.51. データ型変換 | R Financial & Marketing Library

同じことを51列目、52列目(別のsampleのgenotype fields)にも行い、それぞれgenotype2, genotype3という名前のベクトルを作ります。

(2) 変換したgenotypes(数値)をtable("MHC_test4_annotated")に結合する

まず、genotype1からgenotype3を、Rの"cbind"コマンドを使って結合し、一つのtable(m x nの表)を作成します。以下のコマンドを実行します。

# genotype1からgenotype3を表として結合
genotypes <- cbind(genotype1, genotype2)
genotypes <- cbind(genotypes, genotype3)

次に、以下のコマンドを実行することで、このtableをVCF fileのtableである"MHC_test4_annotated"に結合します。

# table genotypesをもとのVCF fileに結合
MHC_test4_annotated <- cbind(MHC_test4_annotated,genotypes)

結合されていることを、headコマンドで確認します。

# 確認
head(MHC_test4_annotated,10)

参考URL:
R言語で統計解析入門: サンプルの結合(行追加)と変数の結合(列追加) 梶山 喜一郎


(3) ALT_alleleの「対立遺伝子頻度」を計算する

ここではALT alleleの対立遺伝子頻度を計算します。ALT alleleの対立遺伝子頻度は、変換したgenotypeの数値を使って


(ALT_allele_freq) = (genotype1 + genotype2 + … + genotypeN) / (2 x number_of_individuals)


で計算できます。ここではN = 3なので、以下のコマンドを入力します(ついでに結果を"MHC_test4_annotated"に結合しています)。

# allele頻度 // referenceと異なるSNP alleleの頻度を求める
alt_allele_freq <- (genotype1 + genotype2 + genotype3) / 6	#allele freqを計算
# VCF fileに結果を結合
MHC_test4_annotated <- cbind(MHC_test4_annotated,alt_allele_freq)

(4) 各Genotypeごとの「遺伝子型頻度」を計算する

遺伝子型頻度の計算は、ちょっとややこしいです。各SNPごとに、特定のgenotypeの個体数を数える必要があります。とりあえず以下のコマンドを実行して下さい(ALT_alleleのホモ接合体の頻度をSNPごとに計算)。

# Tableで行ごとに特定のgenotype個数を数える
# 結果はベクトルで

# number of samples
no_samples <- 3

# genotype 1/1 -> alt_homo
freq_alt_homo <- c()
x <- 0
for (i in 1:length(genotype1)) {
	x <- length(which(genotypes[i,]==2))/no_samples
	freq_alt_homo <- c(freq_alt_homo,x)
}

ヘテロ接合体、REF alleleのホモ接合体も同様にして計算できます。

結果を結合してtableを作成、VCF fileに結合します。

# 結果を結合 // tableにする
freq_genotypes <- cbind(freq_alt_homo, freq_hetero)
freq_genotypes <- cbind(freq_genotypes,freq_ref_homo)

# 結果をVCF fileに結合
MHC_test4_annotated <- cbind(MHC_test4_annotated,freq_genotypes)

長くなったので、続きは次のエントリで。ここまでのRコマンドはこちらの"genotype.R"にすべて書いています。

MiSeqデータのMHC領域へのマッピング(14):Variant Effect PredictorによるSNPのアノテーション(その3)

2015-05-18

ここでは、Variant Effect Predictorで得られたアノテーション情報を含むVCF formatファイルをダウンロードして、Rで読み込む方法を解説します。

4. アノテーション情報を含むVCFファイルをダウンロードする

解析結果の画面(前回の内容参照)で、設定部分の右側"Download"(下図)から、"VCF"をクリックします。

f:id:hashiyuki:20150518160119j:plain

ダウンロードされるファイルを別名で保存します。ここでは"MHC_test4.var.annotated.vcf"とします。このファイルをCotEditorなどで開くとわかりますが、このファイルは改行コードがCR/LFになっているので、Terminal等で開くことも考えて、LFに変換しておきます。

ちなみに"Download"の下の方に見えている"Excel file"のような絵(上図)をクリックすると、.csv形式でファイルが保存され、Excelで開くことができます(Excelを使う場合、これも悪くない)。

5. アノテーション付き".vcf"ファイルをRで読み込む

ダウンロードしたVCFファイル"MHC_test4.var.annotated.vcf"は、Rで読み込むことができます。ただし、Rに読み込む前に一手間かける必要があります。

(1) まず、"MHC_test4.var.annotated.vcf"をCotEditorで開きます。

(2) 「検索」メニューから「検索…」を選択します。

(3)「検索と置換」ウィンドウの「正規表現」をチェックし、検索文字列の窓(上)に"\t"と入力、置換文字列の窓(下)に"|"と入力します(下図)。「すべて置換」をクリックして、"MHC_test4.var.annotated.vcf"のすべてのtabを"|"に変換します。

f:id:hashiyuki:20150518165449j:plain

(4) ファイルを上書き保存します。改行コードがLFになっていることを確認。

(5) Rを起動し、「その他」- 「作業ディレクトリの変更…」で、"MHC_test4.var.annotated.vcf"のあるディレクトリを選択します。

(6) 以下のコマンドで、"MHC_test4.var.annotated.vcf"の内容をRに読み込みます。

MHC_test4_annotated <- read.table("MHC_test4.var.annotated.vcf",header=F, sep="|")

エラーが出なければ、以下のコマンドで読み込まれていることを確認します。

head(MHC_test4_annotated, 10)
              • -

データを読み込めれば、Rを使ってTableの必要な部分のみを抜き出したり、特定の遺伝子や特定のタイプの変異を選択的に得ることができます。

たとえば、"HLA_DPB1"のSNPのうち、アミノ酸を変えるもの (missense_variant)のみを抜き出したい! という場合には、

HLA_DPB1_missense <- subset(MHC_test4_annotated,V11=="HLA-DPB1" & V9=="missense_variant")

とします。

条件に一致するSNP数をカウントするには

length(HLA_DPB1_missense[,1])

とします。

MiSeqデータのMHC領域へのマッピング(13):Variant Effect PredictorによるSNPのアノテーション(その2)

2015-05-18

Variant Effect PredictorによるSNPのアノテーション結果について説明します。

3. アノテーション結果を見る

(1) Summary statistics
解析がうまく行くと、このような結果の画面を見ることができます。

f:id:hashiyuki:20150511184240j:plain

ここでの表と円グラフは、結果全体の要約になっています。"Consequences (all)"は、遺伝子以外の領域も含めたすべてのSNPの位置の存在比率を、"Coding consequences"は遺伝子コード領域のSNPにおける同義・非同義、ナンセンス変異などの存在比率を示しています。

(2) 結果の絞り込み
これら"Summary statistics"の下の部分には、各SNP(多型)ごとのアノテーション結果が、表として示されています。ただし、SNPの数が多い場合には、結果そのままの状態では必要な情報を得ることが難しいです。そこで、表の上にある設定を行う部分(下図)を使用して、必要なSNP情報の絞り込みや、必要なデータのダウンロードを行います。

f:id:hashiyuki:20150511190328j:plain

"Navigation"では、一度に表示するSNPsの数を選択します。Defaultでは5になっているので、もっと多く表示したい場合は、50などに変更します。ただし、SNPsの数が多い(1000以上あるような)場合に"all"を選択してしまうと、表示に時間がかかってしまう場合があります。

次に、結果の表の各列について説明します。

Consequence: 多型variantがどのようなタイプの変異を引き起こすのか、についての情報です。intron_variantはintron内の変異、synonymous_variantは同義変異(アミノ酸を変えないコード領域内の変異:コドンの3rd positionの多くや1st positionの一部の変異)、missense_variantはアミノ酸を変える変異、upstream_gene_variantは遺伝子上流(5' UTRなど)の変異、splice_region_variantはexon-intron境界のsplicing siteに生じた変異、stop_gainedはナンセンス変異になります。
Impact: SNPの表現型に対する影響? LOWは表現型にはあまり関係ない(同義変異など)、MODERATEはある程度関係する(アミノ酸置換など)、MODIFIERは現時点ではわかりません(<--あとで調べること)。
Gene: ENSEMBL gene IDのようです。リンクになっています。
Amino acids, Codons: コード領域のSNPについて、アミノ酸(ミスセンス変異の場合)およびコドンの変異を示しています。
Existing variations: そのSNPが既知の変異に対応する場合、そのSNP IDなどへのリンク。
SIFT: アミノ酸配列の変異が、タンパク質の機能に与える影響を予測するアルゴリズム。予測の方法は、類似した配列を持つ複数のタンパク質の配列の比較に基づいているらしい(詳細はWeb siteを参照)。有害な(deleterious)あるいは許容されない(intolerant)変異は赤で表示されます。
PolyPhen: SIFTと同様に、coding領域の変異がタンパク質の機能に与える影響を予測するツール。SIFTより情報量が多いように見えます。(Website)。
GMAF: 1000 Genomes projectのデータにおけるminor alleleの頻度
AFR-MAF, etc.(その後8列):1000 Genomes projectで調べられている人類集団ごとのSNP allele頻度(existing allele --> majorなほうのallele?)。

これら各列の情報をもとに、表の上部にある設定部分の"Filters"(以下)を使って、必要な情報の絞り込みを行います。

f:id:hashiyuki:20150512112753j:plain

"Consequence"が"missense_variant"である多型のみを表示したい場合、"Uploaded variation"のプルダウンメニューから"Consequence"を選択し、"is" はそのままで、入力部分に"missense_variant"と入力して、"add"をクリックします。

f:id:hashiyuki:20150512113754j:plain

データが絞り込まれ、上のようにFilter(絞り込み)の内容が表示されます。ここから、さらに絞り込むことができます。たとえば、missense_variantのうち、PolyPhemの数値が0.9より大きいものを見たい場合は、"Uploaded variation"のプルダウンメニューから"PolyPhem"を選択し、"is"から">"を選択し、入力欄に0.9を入力して、"add"をクリックします。

Filterの条件を変えるときは、条件の右にある鉛筆を、Filterを削除する場合は、×をクリックします。

アノテーション結果のダウンロード、Rでの編集については、次回に続きます。

MiSeqデータのMHC領域へのマッピング(12):Variant Effect PredictorによるSNPのアノテーション(その1)

2015-05-18

以下の記事も参考にしています。

ショートリードの憂鬱 - 次世代シーケンサー: Ensembl - Variant Effect Predictor でSNPアノテーション

ここでは、データをリファレンスにマッピングすることで得られた多型(SNP)サイト(7回目の内容も参照)の情報について、個々のSNPに付随する情報(ゲノム上の位置、遺伝子、exon or intron, 同義 or 非同義置換、ミスセンス/ナンセンス変異など)、また既知のSNP IDと一致するかどうかなどを得る方法を説明します。ここでは、VCF formatで示された多型サイト情報を使用します。

1. Variant Effect Predictor

Variant Effect PredictorEnsemblが提供するウェブ上の多型アノテーションツールで、ヒト多型サイトについて以下のような情報を得ることができます。

  • 多型によって影響される遺伝子 (gene)または転写産物 (transcript)
  • 多型の位置 (e.g.転写産物の上流, コード領域, non-coding RNA, 調節領域)
  • タンパク質コード領域に与える影響 (e.g. 終止コドンの生成, アミノ酸の置換, 終止コドンの喪失, フレームシフト)
  • 既知のSNPにおけるSNP ID, 1000 genome project から見つかったSNPのminor allele頻度
  • タンパク質配列の変化に対するSIFT及びPolyPhen scores(<-- あとで説明)

2. Variant Effect Predictorを実行する

(1) まず、"Launch"をクリックします。
(2) 入力画面に行きます。以前に他の解析を行っている場合は、以前の結果に対するリンクが出るので、その場合は"New JEP job"をクリックします。
(3) "Or upload file"で「ファイルを選択」します。ここでは7回目で作った"MHC_test4.var.vcf"を選択します。
(4) 結果のフィルタリングをします。"Filters"の"Restrict results"の部分で、"Show one selected consequence per variant"を選択します。これを選択しないと、1つのSNPについて複数アノテーションが出てしまい、見づらいです。
(5) "Run >"をクリックして解析を実行します。

解析結果の解釈は次のエントリで説明します。

MiSeqデータのMHC領域へのマッピング(11):VCF fileについて(1)

2014-04-13

VCF fileのフォーマットについては、以下のURLにあるPDF fileに説明されています。ここでは、VCF formatによるDNA多型(SNPおよびindel)の記述方式を読み解いていきます。

http://samtools.github.io/hts-specs/VCFv4.1.pdf

まず、VCF formatについて、上の方にあるheader行を除くと、以下のような形式になっています。(前回の内容も参照)

CHROM	POS	ID	REF	ALT	QUAL	FILTER	INFO	FORMAT	DRR003760.sort.bam	DRR003761.sort.bam	DRR003762.sort.bam
6	3427377	.	ATTACTTA	ATTA	71.2	.	INDEL;DP=2;VDB=0.0061;AF1=1;AC1=6;DP4=0,0,0,2;MQ=60;FQ=-35	GT:PL:GQ	1/1:110,6,0:7	0/1:0,0,0:3	0/1:0,0,0:3
6	29829653	.	G	T	35.8	.	DP=13;VDB=0.0599;AF1=0.531;AC1=3;DP4=1,3,3,6;MQ=21;FQ=-8.43;PV4=1,7e-08,1,1	GT:PL:GQ	0/1:0,0,0:3	1/1:76,27,0:24	0/0:0,12,29:10
6	29829659	.	T	C	77.4	.	DP=13;VDB=0.0459;AF1=1;AC1=6;DP4=0,0,4,8;MQ=22;FQ=-32.7	GT:PL:GQ	1/1:0,0,0:8	1/1:88,27,0:34	1/1:25,9,0:16

上の例について、Web上の表示では列があっていないので、コピペしてRで読み込んで開いて下さい。(前回の内容も参照)

まず、一番上のheader行(各列の内容を示す)を説明します。

(1) FIXED fields
以下8つのカテゴリは必ず存在する。抜けている値(missing value)は"."で示される。
CHROM: Chromosome番号(ex. 6 <-- chr.6を示す)
POS: 多型サイトの位置(ex. 3427377 <-- 3427377番目の塩基)
ID: 既知の多型ID、たとえばdbSNPの情報などがある場合に、既知の多型サイトであれば、そのID番号が入ります)
REF: Reference配列(ヒトゲノム配列など)における塩基
ALT: Mappingした配列にみられる、referenceと異なる塩基
QUAL: その塩基サイトのPhred scoreを示します。Phred scoreについては英文Wikipediaの説明が参考になります。数字が大きい方がqualityは高い(20 or 30が信頼性の閾値となることが多い)です。
FILTER: (わからない。調べること)
INFO: ここ重要です! ここは、多型サイトについてのいろんな情報を含んでいます。ここでは、重要な部分だけ説明します(詳しくは別エントリで挑戦します)。

INDELIndelの多型であることを示す(特に記述がない場合はSNP)
DPDepth(そのサイトが何回読まれているか)を示す。DP=13なら13回読まれている
VDBVariant distance bias. Referenceにreadがマッピングされたときに、多型がどれくらいランダムに分布するかを見る指標(らしい)。アラインメント(マッピング)の正確性の指標(と思われる)。以下も参照
AF1ALT alleleの頻度(allele frequency)。この値はread自体のデータから推定されるもので、callされた各個体のgenotypeから推定されるものではないことに注意
AC1ALT alleleの数
DP4"="のあとの4つの数字はそれぞれ、REFのread depth(forward側), REFのread depth(reverse側), ALTのread depth(forward), ALTのread depth(reverse). 例:DP4=0,0,4,8(ATL allele のforwardからのread, reverseからのreadがそれぞれ4, 8)
MQRMS mapping quality (Root-mean-square mapping quality of covering reads): リンク先も参照
FQ"Phred probability of all samples being the same"
PV4P-values for strand bias, baseQ bias, mapQ bias and tail distance bias(ちょっとわからない <--あとで調べる

(2) Genotype fields
複数のサンプル(個体)が含まれるとき、サンプルごとに記述される部分。
FORMAT: genotype fieldsのデータ形式とデータ表示の順序を定義します。これではわかりにくいので、ちょっと具体的に説明します。

FORMAT       DRR003760.sort.bam   DRR003761.sort.bam   DRR003762.sort.bam
GT:PL:GQ     1/1:110,6,0:7        0/1:0,0,0:3          0/1:0,0,0:3

FORMAT列のGTというのが、DRR003760.sort.bam列の1/1に、PLが110,6,0に、GQが7に、それぞれに対応しています。Colon ":"が区切りになっています。他のサンプル(個体)も同じです。GT, PL, GQそれぞれの意味は以下です。

GTGenotype: 0がREF alleleで、1がALT alleleを示す。0/0はREFのホモ接合体、0/1はREFとALTのヘテロ接合体、1/1はALTのホモ接合体を示す
PLPhred-scaled genotype likelihoods(を整数値に丸めたもの)。要するに、各genotypeのもっともらしさを数値化したもの。おそらく左から1/1, 0/1, 0/0に対応
GQ条件付き (conditional) genotype quality