2015-08-05
ヒトゲノムのリファレンス配列は、GRCh37またはhg19と呼ばれるデータセットがよく使われています（2013年12月以降は、GRCh38 (hg38)が最新のようです）。
GRCh37（またはGRCh38）は、UCSC Genome Browserのページからダウンロードできます。
UCSC Genome Browser -> Genome -> Downloads とメニューをたどっていくと、このページに行きます。このページの下の方に"Human genome"というカテゴリがあるので、そこからFull data setをクリックすると、ダウンロードのページに移動できます。
ここでは、"chromFaMasked.tar.gz"をダウンロードします（GRCh37の場合）。展開すると各染色体ごとのfasta fileの入ったフォルダができます。
（2015-08-06追記）あとでVariant Effect Predictorを使う場合は、対応するリファレンスゲノム配列がGRCh38.p3なので、こちらを使用する必要があります。
2015-08-06
まず最初に、マッピングを行うNGSのサンプルデータを入手します。
このデータを使用します。
HapMapで使用されている、東京に住む日本人のデータになります。Illumina HiSeq 2000で読まれたデータ（ペアエンド）で、.sraフォーマットで8.9GBあります。詳しくはこのURLも参照して下さい。
データはこのFTPサイトに置かれています。Firefoxなどのブラウザなら、直接FTPサイトからダウンロードできます。Safariを使用している場合は、直接FTPサイトに接続できないので、ターミナルからftpコマンドを使って接続します。
以下のようにします。
(1) まず、データをダウンロードしたい（自分のMacの）ディレクトリに移動します。
たとえば、今日の日付フォルダ（"150806"など）を適当な場所に作成して、cdコマンドでそのフォルダ内に移動します。
(2) ターミナルのシェルに、以下のように入力します。
ftp ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByRun/sra/ERR/ERR251/ERR251633/
うまくつながれば、一連の文字が表示されたあとに
… 250 CWD command successful ftp>
と、"ftp>"コマンド入力を受け付けるモードになります。
このとき、"ERR251633/"の最後の"/"を忘れると
221 goodbye
と表示されて、FTPサイトから抜けてしまうので要注意です。
(3) FTPサイトのディレクトリの位置を確認します。"dir"コマンドを使います。
dir
と入力して
229 Entering Extended Passive Mode (|||50373|) 150 Opening ASCII mode data connection for file list -r--r--r-- 1 ftp anonymous 8862493705 Apr 10 2013 ERR251633.sra 226 Transfer complete
のように表示されればOKです。
データのダウンロードは、"get"コマンドを使って以下のように行います。
get ERR251633.sra
8.9GBあるので、データのダウンロードにはそれなりに時間がかかります（一晩くらいかかる）。
ちなみに、データのダウンロードが終わってFTPサイトを抜けるときは
bye
と入力します。
追記（LFTPを使う）："LFTP"が使えれば、こちらのほうが早いと思います。LFTPはこのサイトから入手できます。Mac OSの場合、packageを展開するだけでインストールできました。使い方は、コマンドをftpからlftpにするだけで、あとは同じです。
2015-08-06
今回は、NGS (Illumina HiSeq) によって配列決定されたヒト全ゲノムのデータを、ヒトの参照ゲノム配列にマッピングします。ここでの目標は、以下の3つです。
1. BWAを用いて、NGS readsをヒトゲノム参照配列にマッピングする
2. SNP情報を.vcf fileとして抽出する。Variant Effect Predictorを使ってアノテーションを付ける
3. 遺伝子の塩基配列を得る。各遺伝子にアノテーションを付ける
4. NGSで決定したヒト個体のゲノムの塩基配列をFASTA形式で得る
以下のエントリで、それぞれの作業について考えていきます（次回へ続く・・）。
以下のエントリーの続きです（じつに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の作成）。
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
