NGSデータ解析まとめ

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

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 >"をクリックして解析を実行します。

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

RNA-Seqデータを用いた系統解析 (1): 解析の方針

RNA-Seqで得られる多数の異なる遺伝子座のデータをもとに系統推定をする方法について、以下では考えていきます。以下には、解析の方針を箇条書きにしてみます(変更の可能性あり)。

(1) RNA-Seqによるデータの入手(Illumina MiSeqなどを使用)

(2) Reference配列の作成(非モデル生物の場合:Trinityによるde novo assembly)

(3) Reference配列に対する近縁種データのマッピングSTAMPYを使用)

(4) Mapping結果からcontigを抽出(samtoolsを使用)、contigのquality checkなど

(5) Contigのアラインメント、遺伝子の読み枠(open reading frame, ORF)の同定、BLAST検索によるアノテーション(いくつかのPerl scriptを使用)

(6) すべてのcontigを統合したデータから系統樹を作成する(RAxMLを使用)

(7) Bayesian concordance analysisによるconcordance treeの作成、種系統樹と異なる分岐を示す遺伝子の同定(BUCKyなどを使用)

以下のような内容でまとめていきます。

解析を行う上でのサンプルデータは、Evolutionの以下の論文(Cui et al. 2013)で発表されたXiphophorus fishの複数種データを使用する予定です。解析の流れについても、大体この論文を参考にしています。データはここからダウンロードできます(.sra file)。Brain由来のRNA-Seq readのデータ(Illumina HiSeq 2000による)になります。

FASTA形式ファイルの(一括)変換

FASTA形式の塩基配列(アミノ酸配列)ファイルをNEXUS形式やPHYLIP形式に変換したい、ということはよくあります。単一のファイルならClustalwMacCladeを使用して変換できますが、多数のFASTAファイルを扱う場合(数百〜数千遺伝子を用いた系統解析など)、一括でファイル形式を変換するスクリプトがあれば便利です。

ここでは、FASTA形式 --> NEXUS形式、またはFASTA形式 --> PHYLIP形式へ変換するPerlスクリプトを公開します。アラインメント自体は行わないので、すでにアラインメントされたFASTA形式ファイルのみが対象になります。また、配列名は9文字まで、という制約があります。

(1) fasta_to_nex.pl

#!/usr/bin/perl
# fasta_to_nex.pl
# 2015-04-07 Y Hashiguchi

# fasta format をnexus formatに変更する
# interleaveしない
# sequenceの開始位置をそろえる
# ntax, ncharをカウント

# このversionはMrBayes blockを書き加えない // 150407
# comment outはあり

## usage ##
# ./fasta_to_nex.pl [infile.fa|.aln|.fasta|.fas] > [outfile.nex]

#

use strict;
use warnings;


### Define variables ##

my $fname = $ARGV[0];

my @lines = ();		# infileの各行を入れる配列
my @seq_with_name = ();	#名前付きsequenceを入れる配列

my $ntax = 0;		# taxa数
my $nchar = 0;		# character数(gapやNも含むsequence長さ)
my $seq_data = 0;	# sequence with name
my $datatype = 'dna';	#data type (dna or protein)

my %sequence;	#sequence nameおよびsequenceを入れるhash

##

### MAIN ###
# input file (FASTA) の読み込み @linesに入れる
open (FASTA, "$fname") || die "Cannot open file:$!";
while(<FASTA>) {
	chomp;
	push(@lines, $_);
}
close(FASTA);

# taxa数をcount (= ">"の数を数える)
foreach my $line (@lines) {
	$ntax++ if($line =~ /\>/);
}

$sequence{"Name"} = shift(@lines);	#最初の1配列の名前(1行目)は$sequence{"Name"}に入れる
$sequence{"Name"} =~ s/\>//;	# ">"を除く

while ($#lines > -1) {
	my $data = shift(@lines);
	if($data =~ /\>/) {
		$nchar = length($sequence{"Sequence"});	# test
		while(length($sequence{"Name"}) < 10) {		#長さが10文字まで // ここでは、nameの長さが10文字より少ないことを仮定している
			$sequence{"Name"} .= " ";	# 後ろに" "を加える
		}
#		print $ntax, "\t", $nchar, "\n";		# test
		$seq_data = $sequence{"Name"} . $sequence{"Sequence"};	#nameとsequence統合
		push(@seq_with_name, $seq_data);	#配列 @seq_with_nameに結果入れる
#		print $sequence{"Name"}, $sequence{"Sequence"}, "\n";		# test
		$sequence{"Sequence"} = "";		# sequence中身を削除
		$sequence{"Name"} = $data;	#sequence_name代入
		$sequence{"Name"} =~ s/\>//;	# ">"を除く
	} else {
		$sequence{"Sequence"} .= $data;
	}
}

# 最後の1配列について
while(length($sequence{"Name"}) < 10) {		#長さが10文字まで 
	$sequence{"Name"} .= " ";	# 後ろに" "を加える
}
$seq_data = $sequence{"Name"} . $sequence{"Sequence"};	#nameとsequence統合
push(@seq_with_name, $seq_data);	#配列 @seq_with_nameに結果入れる

## print out as ".nex" format
print '#NEXUS', "\n\n";
print '[', "\n", 'Data from:', "\n", $ARGV[0], "\n", ']', "\n\n\n";
print 'begin data;', "\n    ", 'dimensions ntax=', $ntax,  ' nchar=', $nchar, ';', "\n";
print '    format datatype=', $datatype, ' interleave=no gap=-;', "\n";
print '    matrix', "\n";
foreach(@seq_with_name) {		# test
	print $_, "\n";
}
print '    ;', "\n";
print 'end;', "\n";

(2) fasta_to_phy.pl

#!/usr/bin/perl
# fasta_to_phy.pl
# 2015-04-07 Y Hashiguchi

# fasta format をPHYLIP formatに変更する
# fasta_to_nex.plの副産物としてできた(おまけ)
# interleaveしない
# sequenceの開始位置をそろえる
# ntax, ncharをカウント
# alignment済みのfasta file (.aln)のみ対応、これ自体はalignmentは行わない

# comment outはあり

## usage ##
# ./fasta_to_nex.pl [infile.fa|.aln|.fasta|.fas] > [outfile.phy]

#

use strict;
use warnings;


### Define variables ##

my $fname = $ARGV[0];

my @lines = ();		# infileの各行を入れる配列
my @seq_with_name = ();	#名前付きsequenceを入れる配列

my $ntax = 0;		# taxa数
my $nchar = 0;		# character数(gapやNも含むsequence長さ)
my $seq_data = 0;	# sequence with name

my %sequence;	#sequence nameおよびsequenceを入れるhash

##

### MAIN ###
# input file (FASTA) の読み込み @linesに入れる
open (FASTA, "$fname") || die "Cannot open file:$!";
while(<FASTA>) {
	chomp;
	push(@lines, $_);
}
close(FASTA);

# taxa数をcount (= ">"の数を数える)
foreach my $line (@lines) {
	$ntax++ if($line =~ /\>/);
}

$sequence{"Name"} = shift(@lines);	#最初の1配列の名前(1行目)は$sequence{"Name"}に入れる
$sequence{"Name"} =~ s/\>//;	# ">"を除く

while ($#lines > -1) {
	my $data = shift(@lines);
	if($data =~ /\>/) {
		$nchar = length($sequence{"Sequence"});	# test
		while(length($sequence{"Name"}) < 10) {		#長さが10文字まで 
			$sequence{"Name"} .= " ";	# 後ろに" "を加える
		}
#		print $ntax, "\t", $nchar, "\n";		# test
		$seq_data = $sequence{"Name"} . $sequence{"Sequence"};	#nameとsequence統合
		push(@seq_with_name, $seq_data);	#配列 @seq_with_nameに結果入れる
#		print $sequence{"Name"}, $sequence{"Sequence"}, "\n";		# test
		$sequence{"Sequence"} = "";		# sequence中身を削除
		$sequence{"Name"} = $data;	#sequence_name代入
		$sequence{"Name"} =~ s/\>//;	# ">"を除く
	} else {
		$sequence{"Sequence"} .= $data;
	}
}

# 最後の1配列について
while(length($sequence{"Name"}) < 10) {		#長さが10文字まで 
	$sequence{"Name"} .= " ";	# 後ろに" "を加える
}
$seq_data = $sequence{"Name"} . $sequence{"Sequence"};	#nameとsequence統合
push(@seq_with_name, $seq_data);	#配列 @seq_with_nameに結果入れる

#print $sequence{"Name"}, $sequence{"Sequence"}, "\n";	#最後のsequenceをprint

print $ntax, " ", $nchar, "\n";
foreach(@seq_with_name) {
	print $_, "\n";
}

複数ファイルを一括で変換する場合は、以下のshell scriptを使用して下さい。カレントディレクトリにあるすべての".aln"拡張子を持つFASTA形式のファイルを変換します(注意:clustal formatには対応していません! 拡張子.alnでも)

(1) FtoN.sh

#!/usr/bin/sh
# FtoN.sh

# FASTA -> NEXUSへのformat変換
# Current directoryに存在するすべてのfileで実行
# alignment済みのファイル (.aln)のみ、alignmentそれ自体は行わない
# mrbayes blockは追加しない
# 結果はphy_files directoryに入る

mkdir nex_files

for fname in *.aln; do

fname=`echo "$fname" | sed -E 's/\.aln//g'`		# 拡張子を除く

# .fasta(.aln) --> .nexに変換 // fasta_to_nex.pl使用
./fasta_to_nex.pl $fname.aln > $fname.nex
mv $fname.nex nex_files

done

(2) FtoP.sh

#!/usr/bin/sh
# FtoP.sh

# FASTA -> PHYLIPへのformat変換
# Current directoryに存在するすべてのfileで実行
# alignment済みのファイル (.aln)のみ、alignmentそれ自体は行わない
# mrbayes blockは追加しない
# 結果はnex_files directoryに入る

mkdir phy_files

for fname in *.aln; do

fname=`echo "$fname" | sed -E 's/\.aln//g'`		# 拡張子を除く

# .fasta(.aln) --> .phyに変換 // fasta_to_phy.pl使用
./fasta_to_phy.pl $fname.aln > $fname.phy
mv $fname.phy phy_files

done

ファイル形式の変換、一般的に、FASTAから別の形式への変換を行うことは多いですが、逆のパターンは少ない気がします。