RNA-Seqデータを用いた系統解析 (2): Reference FASTA fileの作成（Trinityを使用）

以下のエントリーの続きです（じつに96日ぶり！）。

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

非モデル生物で、de novoに配列決定したRNA-Seqデータを系統解析に使用するには、いくつかのアプローチが考えられます。たとえば

(1) すべての種のデータをde novoでアセンブルして、共通のcontigを使って系統樹を書く
(2) 特定の1種でRNA-Seqデータをde novo assembleして、その種で作成したcontigsをリファレンス配列にして、他の種のデータはそのreference contigsにマッピングする。全種でmappngされたcontigで系統樹を書く
(3) 特定の1種で全ゲノム配列を決定し、それをreferenceにする

このエントリでは(2)のアプローチを取ります。(2)の利点は、de novo assembleが1回でいいので、計算時間が短縮されることと、orthologの同定が比較的容易であること（あとで説明します --> 別エントリ）、系統解析に使用する種がある程度近縁な場合は、効率よくデータを得られること、などがあります。あと、全ゲノム決定（3)が必要ないので比較的手軽にできる（コンピュータの計算能力が少なくて済む、経済的？であること）のも、現時点（2015年7月）では利点です。

一方、欠点としては、系統解析に使用する種があまり近縁でない場合、mappingがうまく行かない、または効率的にデータが取れない、ということがあります（このような場合、(1)を使う必要があります）。

1. .sraデータをgetする

このサイトから、まずは

X. birchmanni (SRR751378)
X. gordoni (SRR767711)
X. maculatus (SRR767714)
X. mayae (SRR767715)
X. xiphidium (SRR767782)

のデータを入手して下さい。.sra fileを入手したら、SRA toolkitの"fastq-dump"コマンドを使って.fastq fileを作って（展開して？）下さい。すべてsingle endになっています。

2. ReferenceとなるFASTA fileの作成

注意！：ここでの解析には、Linuxで動くコンピュータが必要です。Trinityを使ってde novo assembleを行うので、メモリは（できれば）128GB以上欲しいところです（可能なら256GB）。

参考URL: Trinity - 井上 潤
Trinityの前処理と動かし方については、OIST井上さんのこのページをほとんど参考にしています。

ここでは、X. birchmanni (SRR751378) のRNA-Seqデータを使用します。大体FASTQ fileで6.96GBあります。

(1) DynamicTrim.plによるデータ精度の評価、low quality readsの除去

まず最初に、FASTQ fileに含まれるqualityの低い配列を除きます。そのためのソフトとして、まずSolexaQAをダウンロード・インストールして下さい。

Shellで以下のコマンドを実行します。

DynamicTrim.pl SRR751378.fastq -h 20 &

(2) LengthSort.plによる短すぎるreads (<50bp)の除去

Shellで以下のコマンドを実行します

LengthSort.pl SRR751378.fastq.trimmed -length 50

(3) Trinityによるde nove assemble

Trinityは、イルミナ社のNGSによるRNA-Seqデータのde novo assembleに特化したassemblerで、比較的短いリード（100bpなど）を得意としているようです。結構頻繁に更新されるので、時々URLをチェックする必要があります。現在は、Web上でも実行できるようです（Linuxがない人は、こちらを使うのも一つの方法です）。

Trinityのダウンロード・インストールはURL参照して下さい（普通は難しくないです）。

とりあえず実行します。以下のコマンドを実行します。

Trinity --seqType fq --single SRR751378.fastq.trimmed.single --JM 200G --CPU 12

オプションの説明：

--seqType : sequence fileのタイプ（FASTQの場合はfq）
--single : single end readsを使用する場合　そのあとに.fastq file名　ペアエンドの場合は--left (seq1) --right (seq2)
--JM : 使用する最大のメモリ容量 // 100G-200Gに設定（使用するコンピュータのスペックによる）
--CPU : 使用するCPUの数（使用するコンピュータのスペックによる）

計算時間は使用するコンピュータの性能にもよりますが、このデータだと、私が研究室で使用しているLinux machine (2.1GHz x 12 core, 256GB RAM)で、大体24時間くらいでした。

計算がうまくいくと、outputとして複数のファイルを含むフォルダ (trinity_out_dir)ができます。

次のエントリに続きます（次は、Trinity結果の整理、reference用のFASTA fileの作成）。

MiSeqデータのMHC領域へのマッピング（17）：RによるVCF fileの操作 (3): SNPのゲノム上での位置と、対立遺伝子頻度との関係をグラフ化する

2015-06-22

3. 各SNPのゲノム上の位置と、allele frequencyとの関係をグラフ化する（散布図の作成）

各SNPごとのallele frequencyを求めたので、ゲノム上での各SNPのallele frequencyの分布を視覚化してみます。ここでは、Rの"plot"関数を使用します。

まず、以下のようにコマンドを入力します。

plot(MHC_test4_freq[,2],MHC_test4_freq[,11],pch=20,cex=0.5,col=4)

plot関数は、引数の1番目（MHC_test4_freqの2列目）と2番目（MHC_test4_freqの11列目）の散布図を作ります。

この状態では、6番染色体のかなり広い範囲でのSNPsのallele frequencyがplotされているので、範囲を絞ります。たとえばHLA-Aの位置はchr6の29942470-29945884なので、この範囲に絞ります。この場合、xlim=c(範囲開始, 範囲終わり)というオプションを指定します。

#位置29942001から29946000までを表示
plot(MHC_test4_freq[,2],MHC_test4_freq[,11],pch=20,cex=0.5,col=4,xlim=c(29942001,29946000))

その他のMHC (HLA)遺伝子の位置はこちら。

参考URL:
plot関数の基本的な使い方：
R言語で統計解析入門： plot()関数を使った散布図の作図体験 梶山 喜一郎
plot関数のオプションなど：
http://cse.naro.affrc.go.jp/takezawa/r-tips/r/53.html

MiSeqデータのMHC領域へのマッピング（16）：RによるVCF fileの操作 (2):VCF fileから必要な情報を抽出する

2015-06-22

2. VCF fileから必要な情報を含む列のみを抽出し、新しいtableを作る

Variant Effect PredictorでアノテーションをしたVCF fileには様々な情報が含まれていますが、そのままでは情報が多すぎて、ちょっと見づらいです。Rを使うことで、必要な情報のみを含む新しい表を作ることができます。

まず、実際にVCF fileのテーブルを開いてみます。前回の、対立遺伝子頻度および遺伝子型頻度データを付与した"MHC_test4_annotated"を使います。

たとえば、必要な情報の列のみからなる新しいtableを作るときは、以下のように入力します。

## 必要な部分だけを抽出したtableを作る
# 位置、SNP alleles, gene name, ENSEMBL ID, status, genotype, allele freq
MHC_test4_freq <- MHC_test4_annotated[,c(1,2,4,5,11,12,15,53,54,55,56,57,58,59)]

MHC_test4_annotatedの1,2,4,5,11,12,15,53,54,55,56,57,58,59列だけからなる新しいtable "MHC_test4_freq"が作成されます。

参考URL:
R：データフレームの列を入れ替える・削除する。 - Qiita

データのtable（データフレーム）の必要な列を抽出、削除、入れ替えなどを行う方法が紹介されています。

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"にすべて書いています。

.bed fileをUCSC Genome Browserで表示する

2015-06-22

Agilent SureDesignによる、captureによるエクソーム解析について、"SureSelect Human All Exon V5"はHLA領域の遺伝子をどのくらいカバーしているかを見てみます。

AgilentのSure Designへの登録と、.bed fileの入手方法についてはHSSの方に頂いた資料（Word file）を参照して下さい。

www.chem-agilent.com

このPDF fileを参考にしました。
Bed fileをgenome browserで確認する方法（User nameとパスワードはHSSの方に頂いたWord fileにあります）。

1. UCSC Genome Browserを開く

genome.ucsc.edu

2. ".bed" fileを読み込む

ここでは、SureDesignの登録サイトからダウンロードした"S04380110_Padded.bed"を読み込みます。画面左上の"Genome Browser"をクリックします（下図）。

左側、家のマークの下です。

移動先の画面の設定部分（下図）で、"add custom track"をクリックします。

移動先の画面で、「ファイルを選択」をクリックし、"S04380110_Padded.bed"を選択してから"Submit"をクリックします（下図）。

うまくいくとデータが読み込まれます（ちょっと時間がかかります）。

3. Genome Browserでプローブの位置を確認する

うまくファイルが読み込まれると、"Manage Custom Tracks"という画面に移動します。画面右側の"Go to Genome Browser"ボタンをクリックすると、実際のGenome Browserの画面に移ります。Browserの画面の上から3行目、"Padded"と書かれている行が、プローブの位置を示しています（下図）。Genome browser内で、"Padded"行の右側に、濃い青のバーで示されているところが、プローブのゲノム上での位置に対応します。

Genome browserに移動した最初の時点では、browserで見えている範囲が狭い（316bp）ので、見える範囲をもっと広げてみます。設定部分の"enter position, gene symbol or search items"と書かれている入力窓に、見たいゲノム上の位置を入れます（下図）。たとえば、1番染色体の1-1,000,000までを見る場合は

chr1:1-1,000,000

と入力して"go"ボタンをクリックします。

Genome上の見たい範囲を移動する場合は、">"や"<"のボタンをクリックします。">"の数が多いボタンは、より長い距離を移動します。

4. ヒトHLA領域に設計されたプローブを確認する

次に、HLA領域をUCSC Genome browser上で見てみます。HLA-Aのヒトゲノム上の位置は

chr6:29942470-29945884

なので、これを入力して"go"ボタンをクリックします。

その他のHLA遺伝子の位置は、このエントリに書いています。