RNA-seqデータ解析(2群間比較, Unstranded, Single-end)
概要
RNA-seqのデータを用いて、siRNAによるノックダウンや特定の刺激により発現量が変動したRNA群を同定する方法(2群間のデータの比較)を紹介します。
ここでは、データのダウンロード方法からLinuxを用いたデータ解析まで、RNA-seqのデータ解析の詳細をStep by Stepで説明します。
ただし、基本的なLinuxのコマンド操作は理解している前提で説明を行います。 パソコン(もしくはスパコン)のメモリは16GB以上を確保することが望ましいです。
使用するデータ
RNA-seq (Single-end, 101bp; Control siRNA, UPF1 siRNA)
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE86148
使用するソフトウェア
- FastQC (v0.11.5)
http://www.bioinformatics.babraham.ac.uk/projects/fastqc/ - FASTX-toolkit (v0.0.14)
http://hannonlab.cshl.edu/fastx_toolkit/ - Bowtie1 (v1.1.2)
http://bowtie-bio.sourceforge.net/index.shtml - TopHat (v2.1.1)
https://ccb.jhu.edu/software/tophat/index.shtml - Samtools (v1.3.1)
http://samtools.sourceforge.net/ - RSeQC (v2.6.4)
http://rseqc.sourceforge.net/ - bedtools (v2.26.0)
http://bedtools.readthedocs.io/en/latest/ - Cufflinks (v2.2.1)
http://cole-trapnell-lab.github.io/cufflinks/ - featureCounts (subread v1.5.0.post3)
http://subread.sourceforge.net/ - edgeR (v3.14.0)
https://bioconductor.org/packages/release/bioc/html/edgeR.html
解析ワークフロー
RNA-seqデータを2群間比較する方法として、Cuffdiffを用いた方法とfeatureCount-edgeRを用いた方法を紹介したいと思います。
TopHat-Cuffdiff workflow
TopHat-featureCounts-edgeR workflow
サンプルコード
Github上に必要なファイルはすべて置いてあります。
https://github.com/Imamachi-n/NGS-Tutorial/tree/master/RNA-seq_Tutorial
RNA-seq_unstr_single-end_1_mapping.sh
TopHatによるマッピング・featureCountsによるリードカウントまでを行うスクリプト。FASTQファイルを引数として指定。
$ ./RNA-seq_unstr_single-end_1_mapping.sh hoge.fastq
10行目: gtfFile
にGTFファイル(アノテーション情報)を指定。
11行目: indexContamFile
にrRNAなどのコンタミ配列のBowtie用のIndexファイルを指定。
12行目: indexGenomeFile
にGenome配列のBowtie用のIndexファイルを指定。
31行目: -r
オプションでHouse-keeping genesの遺伝子領域のBEDファイルを指定。
RNA-seq_unstr_single-end_2_cuffdiff_calc.sh
Cuffdiffによる2群間比較を行うスクリプト。
$ ./RNA-seq_unstr_single-end_2_cuffdiff_calc.sh
8行目: filename
に保存先のディレクトリ名を指定。
9行目: gtfFile
にGTFファイル(アノテーション情報)を指定。
13-14行目: カンマ(,
)区切りでBAMファイル(比較したい2種類)を指定。
17行目: gene_list
にgencode.v19.annotation_filtered_symbol_type_list.txt
(後述)を指定。
RNA-seq_unstr_single-end_2_edgeR_calc.sh
edgeRによる2群間比較を行うスクリプト。
$ ./RNA-seq_unstr_single-end_2_edgeR_calc.sh
8行目: annoList
にgencode.v19.annotation_filtered_symbol_type_list.txt
(後述)を指定。
11行目: saveDir
に保存先のディレクトリ名を指定。
~/custom_command/<hoge.py>
となっている箇所はカスタムのPythonスクリプトが置いてあるパス・ファイル名を指定する。
シェルスクリプトを使用するときの注意点
権限の問題
実行権限をシェルスクリプトに与えた後、実行すること。
$ chmod 755 hoge.sh
改行コードの問題
WindowsとLinuxで改行コードが異なるため、Windows-Linux間でファイルのやり取りをしていると、改行コードの違いに起因する問題が生じることがある。対処法としては、テキストエディタの設定を変えることでLinuxの改行コードを使用する。
ここではAtom
と呼ばれるテキストエディタでの設定方法を説明する。
- [Ctrl] + [ , ](カンマ)キーを押して、Settingsの画面を呼び出す。
- 左のメニューから、
Packages
をクリック。 - Installed Packagesで
line-ending-selector
と検索。 - Core Packagesの項目から、line-ending-selectorの
Settings
をクリック。 - 設定画面で、Default line endingを
LF
に変更。
これで、新規に作成したファイルの改行コードがデフォルトでLF(Linuxの改行コード)になります。
各ステップの説明
0. 必要なソフトウェアのインストール
まず、前述したソフトウェアをbiocondaを用いて自身の解析環境にインストールします。 biocondaのインストールに関しては、下記の記事を参照のこと。
$ conda install fastqc $ conda install fastx_toolkit $ conda install bowtie=1.12 $ conda install tophat $ conda install samtools $ conda install rseqc $ conda install bedtools $ conda install cufflinks $ conda install subread $ conda install bioconductor-edger
注意!!
tophatによるマッピングで最新のbowtieを使うとエラーが出るので、以前のバージョンのbowtie v1.12を代わりにインストールする。
1. SRAファイルのダウンロード
NGSデータはシークエンスされたリード(塩基配列)の情報とクオリティスコアのデータをまとめたFASTQ
形式のファイルとして保存されています。
FASTQファイルの詳細についてはこちらを参照。
https://en.wikipedia.org/wiki/FASTQ_format
FASTQデータは非常に大きいファイルなので、NCBI GEOでは、FASTQファイルを圧縮したSRA
ファイルと呼ばれるファイル形式でデータを公開しています。
まず、下記のURLにアクセスすると、
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE86148
GSEと呼ばれるIDで管理された研究プロジェクトごとのページで、研究の概要や引用されている論文情報を見ることができます。
また、そのプロジェクトで解析したNGSデータがSamples
という項目にリストアップしてあります。
今回は、35サンプルのうち上6つのサンプルを使用します。
GSMから始まるIDをクリックすると、各サンプルの情報を見ることができます。 サンプル名、使用したCell lineや抗体、簡単な実験プロトコルなどの情報が見ることができます。
そこから、さらに下を見ていくと、SRAファイルの置き場所があります。 Downloadの(ftp)をクリックして、SRAファイルをダウンロードします。
SRR4081222をクリック。 ブラウザによって見え方が異なります。ここではGoogle Chromeを使っています。
SRR4081222.sra
のURLをコピー。
wgetで該当するSRAファイルをダウンロードします。
$ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081227 $ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081226 $ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081225 $ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081224 $ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081223 $ wget https://sra-download.ncbi.nlm.nih.gov/srapub/SRR4081222
SRRのIDだけだとわかりにくいので、ファイル名をリネームします。
$ mv SRR4081222 SRR4081222_Control_1.sra $ mv SRR4081223 SRR4081223_Control_2.sra $ mv SRR4081224 SRR4081224_Control_3.sra $ mv SRR4081225 SRR4081225_UPF1_knockdown_1.sra $ mv SRR4081226 SRR4081226_UPF1_knockdown_2.sra $ mv SRR4081227 SRR4081227_UPF1_knockdown_3.sra
SRAファイルからFASTQファイルに変換
先ほどインストールしたSRA Toolkitを使って、NCBI GEOからダウンロードしたSRAファイルをFASTQファイルに展開します。
$ fastq-dump SRR4081222_Control_1.sra $ fastq-dump SRR4081223_Control_2.sra $ fastq-dump SRR4081224_Control_3.sra $ fastq-dump SRR4081225_UPF1_knockdown_1.sra $ fastq-dump SRR4081226_UPF1_knockdown_2.sra $ fastq-dump SRR4081227_UPF1_knockdown_3.sra
2. Quality check
FastQCを使って、FASTQファイルに記載されている各リードのクオリティをチェックします。
FASTQファイルはPhred Quality Scoreと呼ばれる指標によって、各塩基のクオリティを記述しています。
詳しくはこちらのサイトを参照のこと。
https://en.wikipedia.org/wiki/Phred_quality_score
実行例
$ mkdir fastqc_SRR4081222_Control_1 $ fastqc -t 8 -o ./fastqc_SRR4081222_Control_1 \ ./SRR4081222_Control_1.fastq -f fastq
mkdir
コマンドでfastqc_SRR4081222_Control_1
という名前のディレクトリをカレントディレクトリに作ります。
次に、fastqc
を使ってSRR4081222_Control_1.fastq
のデータのクオリティをチェックします。
-t
オプションで使用するCPUのコア数を指定します。コア数を増やせば単純に計算速度が向上しますが、FastQCではメモリを1コアあたり250MB消費するので、要求できるメモリ量を超えるコア数を指定するとエラーが起こります。
-o
オプションで出力先のディレクトリを指定し、-f
オプションでInputファイルのファイル形式を指定しています。
FastQCから出力されるデータが、fastqc_SRR4081222_Control_1
ディレクトリの中のSRR4081222_Control_1_fastqc.html
になります。
ファイルをブラウザで開くと以下のようなグラフを確認できます。
Basic Statistics
基本的な統計量が示されています。
- 全リード数
- リード長
- GC contents
などなど
Per base sequence quality
リードの各塩基のクオリティスコアを示しています。 Phred quality scoreがだいたいグリーンの領域(Scoreが28以上)に収まっているかどうか確認します。 結果として、クオリティが低いリードは含まれていないことが確認できます。
Per tile sequence quality
フローセルの各タイルごとのクオリティスコアを示しています。
Illumina社製の次世代シーケンサーでは、「フローセル」と呼ばれるガラス基板上でDNA合成反応を行います。このフローセルは「タイル」と呼ばれる単位に区切られており、各タイルごとに塩基配列を決定しています。
シーケンスをかけたときに、例えば特定のタイルに気泡やゴミが入っているとクオリティの低下が見られることがあります。
特定のタイルで著しいクオリティの低下が見られる場合は、シークエンス時に上記のような問題があったと考えられます。
詳しくはこちらの資料を参照のこと。
http://nagasakilab.csml.org/ja/wp-content/uploads/2012/04/Sato_Sugar_JPNreview.pdf
Per sequence quality scores
各リードの全長のクオリティスコアの分布を示しています。
Per base sequence content
各塩基のA, T, G, Cの含有量を示しています。 RNA-seqの場合、それぞれの含有量はほぼ25%ずつになりますが、PAR-CLIPのようにRNA結合タンパク質と結合しているRNA配列を抽出してきている場合、それぞれの含有率に偏りが見られます。
Per sequence GC content
リードのGC contentsの分布を示しています。
Per base N content
各塩基中に含まれるN
の含有率(塩基を読めなかった箇所)を示しています。
Sequence Length Distribution
リード長の分布を示しています。
Sequence Duplication Levels
Duplidate readsの含まれている数を示しています。
Overrepresented sequences
頻出する特徴配列が示されています。リード中にアダプター配列などが混入している場合、その配列が示されます。
Adapter Content
各塩基ごとに見たときのリード中に含まれているアダプターの割合を示しています。 あくまで、FastQCに登録されているアダプター配列しか確認していないので、登録されていないアダプター配列を使っていた場合、そのアダプター配列がリード中に混入していても確認できないことがあります。
Kmer Content
特定の塩基配列のパターンがリード中に頻出していないかどうかチェックしています。
3. Quality filtering
FASTQファイル中にクオリティスコアの低いリードが含まれていることが確認された場合、FASTX-toolkitを用いてクオリティの低い塩基及びリードを除去します。
今回のデータでは、クオリティの低いリードが含まれていないのでこの工程をパスします。
Illumina製のシークエンサーでいうと、HiSeq以降のデータはデータのクオリティが非常に高いので、いちいちクオリティフィルターにかける必要がなくなりました。
一応、実行例を示すと以下の通りになります。
実行例
$ fastq_quality_trimmer -Q33 -t 20 -l 18 -i ./SRR4081222_Control_1.fastq | fastq_quality_filter -Q33 -q 20 -p 80 -o ./SRR4081222_Control_1_1_filtered.fastq
fastq_quality_trimmerのパラメータ
-Q
: Phred quality scoreのバージョン指定(Solexa, Illumina 1.3+, Illumina 1.5+といったかなり古いバージョンのFASTQファイルでは'64'を指定。最近のFASTQファイルは'33'と指定しておけばOK)。-t
: 指定したPhred quality scoreを下回る塩基をトリミングする。-l
: トリミング後の最低リード長。
fastq_quality_filterのパラメータ
-q 20 -q 80
: リードの80%
でPhred quality scoreが20
を下回る場合、そのリードを排除する。
4. Quality check
フィルタリング後のFASTQファイルをFastQCを用いて再確認します(今回のデータはフィルタリングを行わないのでこの工程もパス)。
実行例
$ mkdir fastqc_SRR4081222_Control_1_filtered $ fastqc -o ./fastqc_SRR4081222_Control_1_filtered ./SRR4081222_Control_1_1_filtered.fastq -f fastq
5. Mapping
TopHatを用いてGenomeとTranscriptomeにリードをマッピングを行います。
注意!!
RNA-seqのデータ解析については近年、マッピングソフトとしてSTARがデファクトスタンダードになっています(あくまで個人的な意見です。ENDODEプロジェクトで使われて一気にスタンダードになった気がします)。
確かに、STARのほうがTopHatと比較してマッピングの精度が向上しているので、現時点ではSTARがファーストチョイスかと思います。
しかし、STARはメモリを大量に消費し、ヒトゲノムだと32GB以上は必要となります。メモリ不足が懸念される解析環境では、メモリ消費が少ないTopHatを選択するのが無難だと思います(昔はこれがデファクトだったので)。
Indexファイルの作成
TopHatの内部で使用されているBowtieでゲノムへマッピングするために、ヒトゲノムに対するIndexファイルを用意する必要があります。
Bowtieのサイトから、Indexファイルはダウンロードすることもできます。右端のカラムに各ゲノムのIndexファイルが羅列されているので、任意のIndexファイルをダウンロードしてきます。 http://bowtie-bio.sourceforge.net/index.shtml
もしくは、自身でIndexファイルを作成することも可能です。以下では、その方法を説明します。
まず、Human Genome(hg19)のFASTAファイルのダウンロードしてきます。 次に、catコマンドで各染色体ごとのFASTAファイルをマージして、hg19.faという名前で出力します。
$ for file in 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 X Y M do wget http://hgdownload.soe.ucsc.edu/goldenPath/hg19/chromosomes/chr${file}.fa.gz gunzip chr${file}.fa.gz done $ cat chr1.fa chr2.fa chr3.fa chr4.fa chr5.fa chr6.fa chr7.fa chr8.fa chr9.fa chr10.fa chr11.fa chr12.fa chr13.fa chr14.fa chr15.fa chr16.fa chr17.fa chr18.fa chr19.fa chr20.fa chr21.fa chr22.fa chrX.fa chrY.fa chrM.fa > hg19.fa
最後に、hg19.fa
ファイルからBowtie用のIndexファイルを作成します。
$ bowtie-build ./hg19.fa hg19
- 1つ目の引数: リファレンスゲノムのFASTAファイルを指定。
- 2つ目の引数: Indexファイルの名前を指定します。ここでは、Indexの名前を
hg19
とします。
アノテーション情報(GTFファイル)の準備
TopHatでは、GTF形式のファイルのTranscriptomeのアノテーション情報を使って、Splice junctionサイトに対してもマッピングを行うことができます。
今回はTranscriptomeの情報として、Gencodeのアノテーション情報を利用します。
https://www.gencodegenes.org/releases/19.html
上記のサイトにアクセスすると、さまざまなデータが置いてあるが基本的にgencode.v19.annotation.gtfというメインのアノテーション情報を利用すればOKです。
ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz
wgetコマンドで該当のファイルをダウンロードし解凍します。
$ wget ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz $ gunzip ftp://ftp.sanger.ac.uk/pub/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz
カレントディレクトリにgencode.v19.annotation.gtfが保存されます。
GenomeとTranscriptomeへのマッピング
ここまでで、ゲノムのIndexファイルと、アノテーション情報の入ったGTFファイルが用意できたので、いよいよTopHatを用いて各シークエンスリードをゲノムとトランスクリプトームにマッピングします。
以下のコマンドを各サンプル(FASTQファイル)に対して行います。
実行例
$ tophat --bowtie1 -p 8 -o ./tophat_out_SRR4081222_Control_1 \ -G ./gencode.v19.annotation.gtf \ ./hg19 ./SRR4081222_Control_1_1_filtered.fastq
--bowtie1
: マッピングソフトにBowtie1
を使用。-p
: 使用するコア数を指定。-o
: 出力先のディレクトリ名を指定(ディレクトリは作成していなくても自動的に作成してくれる)。-G
: アノテーション情報として、GTFファイルを指定。- 1つ目の引数: 先ほど作成したヒトゲノム(hg19)のBowtie1のIndexファイルを指定。
- 2つ目の引数: マッピングしたいFASTQファイルを指定。
TopHatによるマッピングが終わると、
- align_summary.txt
- accepted_hits.bam
というファイルが出力されます。
align_summary.txt
がマッピング結果をまとめたテキストファイルで、中身を見てみると、
Reads: Input : 13282516 Mapped : 12196063 (91.8% of input) of these: 1887950 (15.5%) have multiple alignments (985 have >20) 91.8% overall read mapping rate.
上記のように、マッピングされたリードの数・割合や、Multi-hitしたリードの数・割合などが記載されています。
accepted_hits.bam
はマッピングされたリードの情報、具体的には、各リードがマッピングされた染色体座標軸、アライメントスコア、リファレンスゲノムと比較したときのIndelやMutationの有無などの情報が記載されています。
ただし、BAMファイルはバイナリファイルなので、ファイルを開いて直接中身を見ることができません。
ファイルの中身を確認したい場合はsamtools
を用いることで、SAM形式に変換してファイルの中身を確認することができます。
$ samtools view ./tophat_out_SRR4081222_Control_1/accepted_hits.bam | less
SAMファイルの中身。
HWI-ST1075L:319:C42GNACXX:1:1315:16479:46448 256 chr1 12137 0 36M * 0 0 CCTGCATGTAACTTAATACCACAACCAGGCATAGGG @@@FFFFFHHHHHJJJJJJIJJJJJGIGEHIGGIEG XA:i:1 MD:Z:0T35 NM:i:1 XS:A:+ NH:i:6 CC:Z:chr15 CP:i:102518998 HI:i:0 HWI-ST1075L:319:C42GNACXX:1:1214:15035:37072 272 chr1 12333 0 36M * 0 0 GGCTGTGACTGCTCAGACCAGCCGGCTGGAGGGAGG JIJJIIIIIJJIGIIIIJJJJJIHHHGHFFFFFCCC XA:i:0 MD:Z:36 NM:i:0 XS:A:+ NH:i:6 CC:Z:chr15 CP:i:102518802 HI:i:0 HWI-ST1075L:319:C42GNACXX:1:1301:7343:9547 272 chr1 12693 1 36M * 0 0 CTCTTGTGAGTGTCCCCAGTGTTGCAGAGGTGAGAG JJIHIJJJJJJIGJJJJJJJJJJHHHHHFFFFFCCC XA:i:0 MD:Z:36 NM:i:0 XS:A:+ NH:i:4 CC:Z:chr15 CP:i:102518442 HI:i:0 HWI-ST1075L:319:C42GNACXX:1:2204:8186:34501 16 chr1 12693 1 36M * 0 0 CTCTTGTGAGTGTCCCCAGTGTTGCAGAGGTGAGAG JJJJJJJJIJJJIIJJJJJJIJJHHHHHFFFFFCCC XA:i:0 MD:Z:36 NM:i:0 XS:A:+ NH:i:4 CC:Z:chr15 CP:i:102518442 HI:i:0 HWI-ST1075L:319:C42GNACXX:1:1305:5797:18762 272 chr1 12957 0 36M * 0 0 CATCCCCTTCACTCCCAGCTCAGAGCCCAGGCCAGG DEHIGGIFF?2F;GGFAIIHIGGHDFHDDFDFD@@@ XA:i:0 MD:Z:36 NM:i:0 XS:A:+ NH:i:6 CC:Z:chr12 CP:i:92621 HI:i:0 HWI-ST1075L:319:C42GNACXX:1:2213:9952:39080 272 chr1 13024 0 36M * 0 0 GTGCATGAAGGCTGTCAACCAGTCCATAGGCAAGCC JJJJJJJJJJIIJIIHIGJJJJJHHGHHFFFFFCCC XA:i:0 MD:Z:36 NM:i:0 XS:A:+ NH:i:7 CC:Z:chr12 CP:i:92554 HI:i:0
BAM/SAMファイルについての詳細を知りたい場合は、下記のURLを参照のこと。
https://samtools.github.io/hts-specs/
https://samtools.github.io/hts-specs/SAMv1.pdf
https://samtools.github.io/hts-specs/SAMtags.pdf
6. データの品質チェック
RNA-seq用のRNAサンプルはバイオアナライザーなどの機械であらかじめRNAの品質をチェックしてからcDNAライブラリを作成しています。
ただ、古いデータの中にはRNAの品質に問題があるものもあり、遺伝子領域を見たときに、例えば、5'側から3'側でのCoverageがガタガタするケースがあります。
その他にも、PolyAセレクションを行ったサンプルの場合、RNAの品質が悪いと3'側に極端にリードが偏る傾向にあることが知られています。
そこで、マッピングされたリードから遺伝子領域上のシークエンスリードのCoverageを調べることによって、大まかにRNAの品質も問題がないか確認を行います。
今回は、RSeQC
のスクリプトの1つであるgeneBody_coverage.py
を用いて、遺伝子領域上でのCoverageを調べてみようと思います。
このスクリプトの実行には、BAMファイルとそのIndexファイル、BED形式で記述された遺伝子領域の情報(アノテーション情報)が必要になります。
まず、アノテーション情報に関しては、基本的にどの細胞でも発現しているHouse-keeping genesのみに絞って解析したほうがよいです。
そこで、RSeQCの開発元が配布しているHouse-keeping genesのBEDファイルのリストを入手してきます。
$ wget https://sourceforge.net/projects/rseqc/files/BED/Human_Homo_sapiens/hg19.HouseKeepingGenes.bed
次に、結果の保存先のディレクトリの用意と、BAMファイルのIndexファイルを作成しておきます。
$ mkdir geneBody_coverage_SRR4081222_Control_1 $ samtools index ./tophat_out_SRR4081222_Control_1/accepted_hits.bam
mkdir
コマンドで新しいディレクトリを作成。samtools index
でBAMファイルのIndexファイルを作成します。引数にはBAMファイルを指定するだけです。
それでは、House-keeping genesの遺伝子領域のCoverageを調べてみたいと思います。
実行例
$ geneBody_coverage.py -r ./hg19.HouseKeepingGenes.bed \ -i ./tophat_out_SRR4081222_Control_1/accepted_hits.bam \ -o ./geneBody_coverage_SRR4081222_Control_1/SRR4081222_Control_1_RSeQC_output
-r
: 先ほどダウンロードしたHouse-keeping genesの情報が格納されているBEDファイルを指定。-i
: Coverageを調べたいBAMファイルを指定(加えて、Indexファイルを同じディレクトリにおいておく必要がある)。-o
: 出力先のディレクトリとファイル名のPrefix。
結果として、指定した名前.geneBodyCoverage.curves.pdf
というファイルが出力されます。ファイルの中身を見るとこんな感じ。
RNAの品質に問題がある(例えば、RNAが壊れている)と、3'側にリードが偏っていたり、cDNAライブラリの構築に使ったRNA量が少なく、PCR時に特定の配列が異常に増幅され(PCR bias、Duplicate readsが多い)、ガタガタしたCoverageの分布を取ることがあります。
こういったデータはRNAの定量性に問題があるため、使用しないように気をつけましょう。
このステップで問題がありそうだと判断されるデータに関しては、下記のUCSC genome browserでの可視化により、個々の遺伝子上にマッピングされたリードの分布に注意して見ていく必要があります。
7. Visualization (For UCSC genome browser)
BedGraphファイルの準備
Bedtools
を利用して、TopHatから得られたBAMファイルを、UCSC genome browser上で可視化できるファイル(Bedgraphファイル)に変換します。
まず、mkdir
コマンドでデータを保存するディレクトリを作っておきます。
$ mkdir UCSC_visual_SRR4081222_Control_1
bedtools genomecovで特定のゲノム領域でのCoverageを計算します(BAMファイルからBedgraph
ファイルを作成)。
Bedgraphファイルの詳細については、下記のURLを参照のこと。
https://genome.ucsc.edu/goldenpath/help/bedgraph.html
$ bedtools genomecov -ibam ./tophat_out_SRR4081222_Control_1/accepted_hits.bam -bg -split \ > ./UCSC_visual_SRR4081222_Control_1/SRR4081222_Control_1_4_result.bg
-ibam
: BAMファイルを入力ファイルとして指定する。-bg
: 出力ファイルの形式をBedGraphファイルに指定。-split
: Splicing juctionにマッピングされたリード(リードが分割されてマッピングされているリード)を考慮する。
標準出力でBedGraph
ファイルが出力されます。
UCSC genome browserへデータをアップロードする際に、ファイルの種類やTrackの名前などを提示する必要があります。
https://genome.ucsc.edu/goldenpath/help/customTrack.html
それらの情報をヘッダー行(ファイルの一行目)に記載します。
$ echo "track type=bedGraph name=RNA-seq_Control_1 description=RNA-seq_Control_1 visibility=2 maxHeightPixels=40:40:20" \ > ./UCSC_visual_SRR4081222_Control_1/tmp.txt
echo
コマンドを使って、ヘッダー行の情報をtxtファイルとして出力しています。
そのあとに、cat
コマンドを用いて、先ほど作ったBedGraphファイルと結合させます。
$ cat ./UCSC_visual_SRR4081222_Control_1/tmp.txt ./UCSC_visual_SRR4081222_Control_1/SRR4081222_Control_1_4_result.bg \ > ./UCSC_visual_SRR4081222_Control_1/SRR4081222_Control_1_4_result_for_UCSC.bg
出力されたファイルをそのままアップロードしてもOKですが、容量が大きいためアップロードに時間がかかります。そのため、アップロードするファイルのサイズを小さくするために、ここではbzip2ファイルに圧縮します。
$ bzip2 -c ./UCSC_visual_SRR4081222_Control_1/SRR4081222_Control_1_4_result_for_UCSC.bg \ > ./UCSC_visual_SRR4081222_Control_1/SRR4081222_Control_1_4_result_for_UCSC.bg.bz2
データのアップロード
UCSC genome browserにデータをアップロードします。
ホーム画面から、[MyData] -> [Custom Tracks]をクリックします。
Add Custom Tracksというページに移動するので、ファイルを選択をクリックし、アップロードするファイルを選択し、submitボタンをクリックします。
アップロードが完了すると、Manage Custom Tracksというページに移動します。右端にあるgoボタンをクリックすると、データをみることができます。
8. RNA-seqデータの2群間比較(Cuffdiffのケース)
Cuffdiffを用いて、siCTRLとsiUPF1の2群間のRNA-seqデータを比較します。
cuffdiff -p 8 --multi-read-correct -o ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode \ /home/akimitsu/database/gencode.v19.annotation_filtered.gtf \ ./tophat_out_SRR4081222_Control_1/accepted_hits.bam,./tophat_out_SRR4081223_Control_2/accepted_hits.bam,./tophat_out_SRR4081224_Control_3/accepted_hits.bam \ ./tophat_out_SRR4081225_UPF1_knockdown_1/accepted_hits.bam,./tophat_out_SRR4081226_UPF1_knockdown_2/accepted_hits.bam,./tophat_out_SRR4081227_UPF1_knockdown_3/accepted_hits.bam
-p
: 使用するコア数を指定。--multi-read-correct
: 定量の精度を上げる目的で、複数箇所にマッピングされたリードに重み付けを行う。-o
: 出力先のディレクトリを指定。- 1つ目の引数: アノテーション情報のGTFファイルを指定。
- 2つ目・3つ目の引数: 2種類のBAMファイルをそれぞれ指定。今回のように、siCTRLとsiUPF1の条件でn=3でデータがある場合、それぞれのBAMファイルのパスをカンマ区切りで並べる(例えば、
CTRL_1.bam,CTRL_2.bam,CTRL_3.bam KD_1.bam,KD_2.bam,KD_3.bam
)。
パラメータの詳しい解説は以下のサイトを参照のこと。
http://cole-trapnell-lab.github.io/cufflinks/cuffdiff/index.html
次に、Gene symbolなどの情報をCuffdiffから出力されるファイルに加える。
遺伝子リストの準備
まず、extract_gene_symbol_type_from_gtf.py
スクリプトを用いて、GTFファイルから遺伝子リストを作成します(Cuffdiffで集計したデータからmRNAとlncRNAに分類分けするために用います)。
以下のURLに必要なスクリプトが置いてある。
https://github.com/Imamachi-n/NGS-Tutorial/blob/master/RNA-seq_Tutorial
$ python2 extract_gene_symbol_type_from_gtf.py gencode.v19.annotation_filtered.gtf gencode.v19.annotation_filtered_symbol_type_list.txt
- 1つ目の引数: GTFファイルを指定する。
- 2つ目の引数: 遺伝子リストが出力される。
次に、先ほど作成した遺伝子リストをもとに、cuffdiffの出力ファイルからmRNAとlncRNAの情報を抽出し、それぞれ別のファイルとして保存する。
# mRNA $ python ~/custom_command/cuffdiff_result.py \ ./gencode.v19.annotation_filtered_symbol_type_list.txt \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/gene_exp.diff \ mRNA \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/cuffdiff_result_mRNA.txt # lncRNA $ python ~/custom_command/cuffdiff_result.py \ ./gencode.v19.annotation_filtered_symbol_type_list.txt \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/gene_exp.diff \ lncRNA \ ./cuffdiff_out_RNA-seq_UPF1_Knockdown_Gencode/cuffdiff_result_lncRNA.txt
- 1つ目の引数: 先ほど用意した遺伝子リストを指定。
- 2つ目の引数: cuffdiffの出力ファイルである
gene_exp.diff
ファイルを指定。 - 3つ目の引数:
mRNA
もしくはlncRNA
を指定(指定したRNA群のデータを取得)。 - 4つ目の引数: mRNAもしくはlncRNAのリストを出力。
8. RNA-seqデータの2群間比較(featureCounts-edgeRのケース)
各遺伝子に対してのリード数のカウント
まず、featureCountsで各遺伝子領域にマッピングされたリード数をカウントします。
$ mkdir featureCounts_result_SRR4081222_Control_1 $ featureCounts -T 8 -t exon -g gene_id -s 2 -a ./gencode.v19.annotation_filtered.gtf \ -o featureCounts_result_SRR4081222_Control_1/featureCounts_result_SRR4081222_Control_1.txt \ ./tophat_out_SRR4081222_Control_1/accepted_hits.bam
mkdir
コマンドを使って、保存先のディレクトリを作り、そのあとfeatureCountsを実行します。
-T
: 使用するコア数の指定。-t exon
: リードのカウントに使用するExonの情報を指定。ここでは、GTFファイル内で3列目のfeature
がexon
である行を使用するように指定している。-g gene_id
: 各行のメタ情報(遺伝子ID、Transcript ID、遺伝子名など)のうち、名前として使用する情報を指定する。ここでは、gene_id
を指定することで、Gencodeで決められた遺伝子IDを使っている。-a
: アノテーション情報(GTFファイル)を指定。-o
: 出力先のパスとファイル名を指定。- 1つ目の引数: BAMファイルを指定。
featureCountsから出力されたリードの集計データのリストのうち、遺伝子IDとリード数の列だけを抽出します。
$ sed -e "1,2d" featureCounts_result_SRR4081222_Control_1/featureCounts_result_SRR4081222_Control_1.txt \ | cut -f1,7 - > featureCounts_result_SRR4081222_Control_1/featureCounts_result_SRR4081222_Control_1_for_R.txt
1, 2行目はHeader行なので、sed
コマンドで除きます。次に、cut
コマンドで遺伝子IDとリード数の列だけを抽出して、データを標準出力で得ます。
edgeRを用いて2群間比較を行う
edgeR_test.R
スクリプトを使って、edgeRを用いた2群間比較を行います。
以下のURLに必要なスクリプトが置いてある。
https://github.com/Imamachi-n/NGS-Tutorial/blob/master/RNA-seq_Tutorial
まず、mkdir
コマンドで保存先のディレクトリを用意する。
mkdir edgeR_result_UPF1_Knockdown
edgeR_test.R
スクリプトを用いて、2群間比較を行う。
Rscript ./edgeR_test.R \ featureCounts_result_SRR4081222_Control_1_ss/featureCounts_result_SRR4081222_Control_1_for_R.txt,\ featureCounts_result_SRR4081223_Control_2_ss/featureCounts_result_SRR4081223_Control_2_for_R.txt,\ featureCounts_result_SRR4081224_Control_3_ss/featureCounts_result_SRR4081224_Control_3_for_R.txt,\ featureCounts_result_SRR4081225_UPF1_knockdown_1_ss/featureCounts_result_SRR4081225_UPF1_knockdown_1_for_R.txt,\ featureCounts_result_SRR4081226_UPF1_knockdown_2_ss/featureCounts_result_SRR4081226_UPF1_knockdown_2_for_R.txt,\ featureCounts_result_SRR4081227_UPF1_knockdown_3_ss/featureCounts_result_SRR4081227_UPF1_knockdown_3_for_R.txt \ Control,Control,Control,Knockdown,Knockdown,Knockdown \ edgeR_result_UPF1_Knockdown \ Yes
- 1つ目の引数: 先ほど用意したリード数を集計した各データをカンマ区切りで指定する。
- 2つ目の引数: 1つ目の引数で指定したサンプルに対して、ラベル付け(例えば、
Control
とKnockdown
)を行う。データはカンマ区切りで指定する。 - 3つ目の引数: 保存先のディレクトリを指定する。
- 4つ目の引数: 複数のサンプル(n=2以上)を取っているかどうか(Duplicateのチェック)。
Yes
もしくはNo
と記入。
アノテーション情報を追加する
edgeRから出力されたデータでは、遺伝子IDしか情報がないので遺伝子名などの情報を付け加えます。
$ python2 ./annotate_gene_symbol_type.py \ ./gencode.v19.annotation_symbol_type_list.txt \ ./edgeR_result_UPF1_Knockdown/edgeR_test_result.txt \ ./edgeR_result_UPF1_Knockdown/edgeR_test_result_anno_plus.txt \ $ python2 ./split_into_each_gene_type.py \ ./edgeR_result_UPF1_Knockdown/edgeR_test_result_anno_plus.txt
annotate_gene_symbol_type.py
スクリプトでアノテーション情報を付加しています。
- 1つ目の引数: GTFファイルから作成した遺伝子リストを指定します。
- 2つ目の引数: 先ほどのスクリプトを使って得られた出力ファイルを指定します。
- 3つ目の引数: アノテーション情報を付加したデータを出力します。
split_into_each_gene_type.py
スクリプトでmRNAやlncRNAなどRNAの種類に応じてデータを分割します。
- 1つ目の引数: 上記のスクリプトから得られた出力ファイルを指定します。
(おまけ)CuffdiffとedgeRの結果を比較する
遺伝子名(ヘッダー行を除く)でソートする。
(head -n +1 cuffdiff_result_mRNA.txt && tail -n +2 cuffdiff_result_mRNA.txt | sort -k 3 ) > cuffdiff_result_mRNA_sorted.txt
(head -n +1 edgeR_test_result_anno_plus_mRNA.txt && tail -n +2 edgeR_test_result_anno_plus_mRNA.txt | sort -k 3 ) > edgeR_test_result_anno_plus_mRNA_sorted.txt
Excel上でマージして比較する。
CuffdiffとedgeRから得られたDEG比較
CuffdiffとedgeRのFold-changeの比較
発現量が低い遺伝子については、cuffdiffとedgeRの結果にばらつきが見られる。
CuffdiffとedgeRの使い分けについて
結局どちらを使ったらよいかという話ですが、好みの問題なのでどちらでも構わないというのが私の答えです(どちらかと言えば、edgeRを勧めますが)。
edgeRを使用する場合、Rのスクリプトを少し書かなければいけないので、プログラミングがわからない人にとっては少しハードルが高いのかもしれません。
その点では、Cuffdiffはコマンドを打つだけでOKなので楽です(プログラミングゼロで結果を出せる!)。
ただ、Cuffdiffに関しては少しクセのあるソフトなので、使用するときには注意が必要だと考えています。
あと、Cuffdiffは他の手法と比較してFalse-positiveが多いという解析結果も出ています(どのデータでも同様の傾向があるかどうかはわかりませんが)。
https://genomebiology.biomedcentral.com/articles/10.1186/gb-2013-14-9-r95
また計算速度もCuffdiffが数時間かかるのに対して、edgeRだと数分で終了します(記事を書いてて、なんかedgeRがいい気がしてきた)。
今までの経験に基づく個人的な意見を以下にまとめました。
mRNAのみを解析対象としているならCuffdiff・edgeRを選択
mRNAのみを解析対象にする場合、Cuffdiffの結果はかなり信頼性が高い(qPCRなどで再現がほぼ取れる)と実感しています。
まぁ、edgeRを使用してもかまわないと思いますが…。
long-noncoding RNAを解析対象としているならedgeRを選択
long-noncoding RNA(lncRNA、PROMPTやeRNAを含む)を解析対象にしている場合、Cuffdiffだとおかしな結果が返ってくることが多々あります。
qPCRで再現が全く取れなかったり、そもそも集計されたリード数とUCSC genome browserで可視化したときのCovarageの具合がどうも一致していない気がする(カウントしてくれたリード数がおかしいんですが、Cuffdiffさん…)などなど…。
GC-contentsが高かったり、Repetitive elementsを含んでいたり、Transcript長が短かったり(1kb未満)するとどうもCuffdiffがうまく計算してくれない気がします(理由の1つとして、Cuffdiff内で行わるBiasの補正が問題なのではと思っていますが、詳しくはよくわかっていません…)。
なので、lncRNAを解析しているなら、CuffdiffではなくedgeRを使うようにしましょう。
Pythonパッケージ: 便利なコマンドラインパーサClickを使ってみる
概要
コマンドラインパーサとして、argparse
と呼ばれるPythonの標準パッケージが存在するが、使い勝手がいまいちだと感じています(あくまで個人的な感想ですが…)。
今回は、たまたま見つけたclick
と呼ばれるサードパーティ製のPythonパッケージを用いてコマンドラインの引数を解釈するコードを書いてみようと思います。
argparse
と比較してclick
の書き方だと、コードが書きやすく、あとから読み返しやすくなると感じます。
Clickのドキュメントはこちら。
http://click.pocoo.org/6/
内容
すでに、こちらのブログに詳しく解説されているのでここでは忘備録的な記載に留めようと思います。
インストール
AnacondaもしくはBiocondaがインストールされていることが前提。
conda install click
基本的な使い方
@click.command()
デコレータを用いて関数を修飾するだけでOK。
print
関数の代わりにclick.echo
を用いて文字列を標準出力することもできる。
@click.command() def cmd(): click.echo("Command was Done.")
続けて、@click.option()
デコレータに任意のオプションを追加していくことができる。argparseと同様に、-i
や--ifastq
などのような-
から始まるコマンドライン引数を指定することができる。
@click.command() @click.option('-i', '--ifastq') def cmd(): click.echo("Command was Done.")
また、作成したオプションに代入された値を取得するために、ifastq
のように-
を付けない形で文字列を指定することで、関数の中で値が代入された変数として扱えます。
以下のように、cmd
関数の引数にifastq
を指定しておく必要がある。
@click.command() @click.option('-i', '--ifastq', 'ifastq') def cmd(ifastq): click.echo(ifastq)
help
でHelpが呼び出されたときの説明文を指定できる。
@click.command() @click.option('-i', '--ifastq', 'ifastq', help='Input file in FASTQ file.') def cmd(ifastq): click.echo(ifastq)
実行例
#/usr/bin/env python import click @click.command() @click.option('-i', '--ifastq', 'ifastq', help='Input file in FASTQ file.') @click.option('-o', '--output-prefix', 'output', help='Prefix of output file.') def cmd(ifastq, output): click.echo(ifastq) def main(): cmd() if __name__ == '__main__': main()
実行すると、
$ python test.py -i test.fastq test.fastq $ python test.py --help Usage: test.py [OPTIONS] Options: -i, --ifastq TEXT Input file in FASTQ file. -o, --output-prefix TEXT Prefix of output file. --help Show this message and exit.
コマンドライン引数のタイプを指定したい
INT型やFLOAT型など、引数のタイプがあらかじめわかっている場合、引数のタイプを指定しておきたいケースが出てくる。
そういった場合、以下のようにargparseのときのように@click.option()
デコレータ内で指定することが可能。
type=<type>
で指定。
@click.option('-m', '--minimum-length', 'minimum_length', type=int, default=18, help='Discard trimmed reads that are shorter than LENGTH. Default: 18')
パラメータタイプは、int、float、bool、strが指定できるようです。
http://click.pocoo.org/6/parameters/
デフォルトの値を決めておきたい
特定のオプションでデフォルト値を指定しておきたいケースがある。 その場合、デフォルト値の設定は以下のように行う。
default='hogehoge'
で指定。
@click.option('-nl', '--nucleotide-trimming-list', 'nucleotide_trimming_list', default='A,N', help='Trimmed nucleotide list. Default: A,N')
指定する引数を限定したい
オプションの引数として、特定の文字列や数値しか受け付けたくない場合、以下のように指定できる引数を限定することができる。
type=click.Choice(['hoge', 'huga'])
で指定。
@click.option('-nd', '--nucleotide-trimming-direction', 'nucleotide_trimming_direction', type=click.Choice(['three', 'five']), help='Trimmed nucleotide direction. Default: three')
Helpパラメータをカスタマイズしたい
Clickのデフォルトでは、Help文は--help
を指定することで出力される。
一般的に、Help文は--help
だけでなく-h
でも出力されるケースが多い。そこで、-h
でもHelp文が出力されるように設定を変更する。
@click.command()
デコレータ内で指定。
CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) @click.command(context_settings=CONTEXT_SETTINGS)
サブコマンド(後述)を利用している場合は、@click.group
デコレータ内で指定する。
CONTEXT_SETTINGS = dict(help_option_names=['-h', '--help']) @click.group(context_settings=CONTEXT_SETTINGS)
コールバック関数で値を検知・Helpを表示
オプションの引数を評価したい場合、コールバック関数を利用すると便利です。
@click.option()
デコレータ内でcallback=<function名>
を指定する。
また、コールバック関数内で、ctx.get_help()
でHelp文を表示させることができる。ctx.exit()
でPythonスクリプトの実行を止めることができる。
def validate_ifastq(ctx, param, value): if not value: click.echo('Error: Do not choose your FASTQ file.\n') print(ctx.get_help()) ctx.exit() return(value) @cmd.command() @click.option('-i', '--ifastq', 'ifastq', callback=validate_ifastq, help='Input file in FASTQ file. [required]')
サブコマンドを用意したい場合
@click.group()
デコレータを用いて、サブコマンドを利用することができます。
#/usr/bin/env python import click @click.group() def cmd(): pass @cmd.command() @click.option('-i', '--ifastq', 'ifastq', help='Input file in FASTQ file.') @click.option('-o', '--output-prefix', 'output', help='Prefix of output file.') def trim_polya(ifastq, output): click.echo('Start trimming polyA sequence from reads...') @cmd.command() @click.option('-i', '--ifastq', 'ifastq', help='Input file in FASTQ file.') @click.option('-o', '--output-prefix', 'output', help='Prefix of output file.') def get_pass_read(ifastq, output): click.echo('Start getting polyA site-supporting reads...') def main(): cmd() if __name__ == '__main__': main()
Help文は以下のとおり。
$ python test.py --help Usage: test.py [OPTIONS] COMMAND [ARGS]... Options: --help Show this message and exit. Commands: get_pass_read trim_polya
サブコマンドごとのオプションをHelp文で出力することもできる。
$ python test.py get_pass_read --help Usage: test.py get_pass_read [OPTIONS] Options: -i, --ifastq TEXT Input file in FASTQ file. -o, --output-prefix TEXT Prefix of output file. --help Show this message and exit.
ENCODEプロジェクトのNGSデータを活用する
概要
ENCODEプロジェクトとは、ポストゲノム研究戦略としてヒトゲノムの機能性エレメントの同定とその全体像を解明するために組織された国際プロジェクトの1つです。
RNA-seq、ChIP-seq、DNase-seq、Hi-C、ChIA-PET、eCLIP-seqなどの解析手法を用いて、毎月膨大なデータを産出しています。
ENCODEプロジェクトでは、共通のプロトコルに基づいた実験と解析によりデータが産出され、即座に一般公開される仕組みを取っています。また、プロトコルも公開されているので自身で内容を確認・再解析することもできます。
ENCODEプロジェクトのデータは、プロジェクトに参加していないラボ・研究機関であっても、ウェブサイト上で公開されているNGSデータに関しては、自由にダウンロードして解析し論文のデータとして使用することもできます。
Data Use, Software, and Analysis Release Policies – ENCODE
実際に、ENCODEのデータを利用して投稿された論文はこれまで1700本にのぼります。
Publications using ENCODE data – ENCODE
ENCODEプロジェクトのデータベースには、有用な情報やデータがいっぱい公開されているのでこれらを使わない手はありません。
そこで今回は、ENCODEプロジェクトに登録されているデータを活用して手軽にデータ解析する方法をまとめました。今回は主に、RNA結合タンパク質(RBP)が結合するRNA領域を網羅的に同定する方法であるeCLIP-seqデータについて見ていきたいと思います。
また、RBPをIPした時に使用した抗体の情報など即座に使える情報も載っているのでそれらについても紹介しようと思います。
欲しいデータを探す
キーワードで検索
まず、ENCODEプロジェクトのウェブサイトに行きます。
ENCODE: Encyclopedia of DNA Elements – ENCODE
右上の検索欄に興味のあるRBPの名前を入力します。
試しにSFPQ
というRBPで検索にかけると、eCLIP-seq、Bind-n-seq、RNA-seq等のデータや、IPに使用した抗体の情報が出てきます。
データの種類から検索
他にも、[Encyclopedia]->[about]から各種データにアクセスできます。
この他にも、いろいろな検索方法が用意されていますが割愛します。
抗体の情報
右端にAntibody
と記載されているのが、抗体の情報です。
その中でも、緑色のマークがついているものは、ENCODEでStandardな抗体として認められた抗体です(抗体の特異性など幾つかの基準をもとに選定されている)。
クリックすると、その抗体でIPしたときのWestern blottingのデータや、shRNAでノックダウンしたときに目的のバンドが減弱しているかどうかWestern blottingで確認したデータなどが公開されています。
以上の情報から、各メーカーの抗体の良し悪しを知ることができます。 ちなみに、以下がENCODEプロジェクトで確認した抗体に関する論文になります。
http://yeolab.github.io/papers/2016/sundararaman_molcell_2016.pdf
Bind-n-seqデータ
Bind-n-seqについては、下記の論文を参照のこと。方法の詳細については、ここでは割愛します。
https://www.ncbi.nlm.nih.gov/pmc/articles/PMC4142047/
リコンビナントタンパク質のQCデータと、Enrichされたモチーフ配列の情報が載っています。
さらに、Processed data
の項目のところに、Enrichされた4, 5, 6, 7ntのモチーフ配列のデータがそれぞれ置いてあります(output typeがenrichment
となっているやつ)。
ダウンロードボタン(Accession
にある小さいアイコン)をクリックすると、該当のファイルをダウンロードできます。
同じページにRawデータ(FASTQファイル)とプロトコルが置いてあるので、自身の手で再解析することも可能です。
eCLIP-seqデータ
1. 必要なファイルをダウンロード
eCLIPのデータを選択し、該当するデータのページに飛ぶと下のほうにProcessed data
の項目があります。
ゲノム上にマッピングされたリードを可視化するためのデータとしてbigWig
ファイルを、各Peakの染色体座標軸に関するデータとしてbed
ファイルをそれぞれダウンロードします。
また、右上のGRCh38
(hg38)をhg19
に変更することで、異なるバージョンのヒトのリファレンスゲノムをもとにして作成したファイルにアクセスできます。
今回は、hg38 (GRCh38)
のデータをダウンロードします。
$ wget https://www.encodeproject.org/files/ENCFF823MPG/@@download/ENCFF823MPG.bigWig $ wget https://www.encodeproject.org/files/ENCFF960OTE/@@download/ENCFF960OTE.bed.gz $ gunzip ENCFF960OTE.bed.gz
ファイル名がわかりにくいのでリネームします。
$ mv ENCFF823MPG.bigWig SFPQ_eCLIP_rep1_ENCFF823MPG.bigWig $ mv ENCFF960OTE.bed SFPQ_eCLIP_rep1_ENCFF960OTE.bed
2. ファイルフォーマットの変換
bedgraphファイルの用意
bigwig
ファイルではUCSC genome browserにアップロードできないので、bedgraph
ファイルと呼ばれる形式にフォーマット変換します。
変換作業を行うために、今回はbigWigToBedGraph
というプログラムを使います。
下記のURLから、目的のプログラムのバイナリファイルを入手します。
http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/
$ wget http://hgdownload.cse.ucsc.edu/admin/exe/linux.x86_64/bigWigToBedGraph $ chmod 755 bigWigToBedGraph
ダウンロードした後、実行権限を与え、パスの通っているディレクトリに入れます。
続いて、bigWig
ファイルをbedgraph
ファイルにフォーマット変換します。
# bigWigToBedGraph <bigWigファイル> <bedGraphファイル> $ bigWigToBedGraph SFPQ_eCLIP_rep1_ENCFF823MPG.bigWig SFPQ_eCLIP_rep1_ENCFF823MPG.bg
- 1つ目の引数: フォーマット変換したい
bigwig
ファイルを指定。 - 2つ目の引数: 出力される
bedGraph
ファイルの名前を指定。
次に、UCSC genome brower用にデータを整えます。ENCODEのデータの中にはハプロタイプのシークエンスにもマッピングしたデータが含まれているので、UCSC genome browserでデータを可視化するために、Chr1-22, chrX, chrY, chrM
からなるCanonical genomeにマッピングされたデータのみを抽出する必要があります。
下記のような適当なスクリプトを書いて、データを整えます。
$ python ENCODE_eCLIP_mod.py SFPQ_eCLIP_rep1_ENCFF823MPG.bg SFPQ_eCLIP_rep1_ENCFF823MPG_for_UCSC.bg
最後に、UCSC genome browserへデータをアップロードする際に、ファイルの種類やTrackの名前などを提示する必要があります。それらの情報をヘッダー行(ファイルの一行目)に記載します。
$ echo "track type=bedGraph name=SFPQ_eCLIP_rep1_ENCFF823MPG description=SFPQ_eCLIP_rep1_ENCFF823MPG visibility=2 maxHeightPixels=40:40:20 color=255,0,0" > ./header_tmp.txt $ cat header_tmp.txt SFPQ_eCLIP_rep1_ENCFF823MPG_for_UCSC.bg > SFPQ_eCLIP_rep1_ENCFF823MPG_for_UCSC_header_plus.bg $ bzip2 -c SFPQ_eCLIP_rep1_ENCFF823MPG_for_UCSC_header_plus.bg > SFPQ_eCLIP_rep1_ENCFF823MPG_for_UCSC_header_plus.bg.bz2
echo
コマンドのところでヘッダー行の情報をtxtファイルとして出力しています。
そのあとに、cat
コマンドを用いて、先ほど作ったBedGraphファイルと結合させます。
出力されたファイルをそのままアップロードしてもOKですが、容量が大きいためアップロードに時間がかかります。そのため、アップロードするファイルのサイズを小さくするために、ここではbzip2
ファイルに圧縮しています。
bedファイルの用意
ENCODEのサイトから得られるbedファイルを、bigWigファイルと同様にUCSC genome browserにアップロードできるファイルに加工します。
下記のような適当なスクリプトを書いて、データを整えます。
$ python ENCODE_eCLIP_mod_bed.py SFPQ_eCLIP_rep1_ENCFF960OTE.bed SFPQ_eCLIP_rep1_ENCFF960OTE_for_UCSC.bed
最後に、UCSC genome browserへデータをアップロードするために、ファイルの種類やTrackの名前の情報をヘッダー行(ファイルの一行目)に記載します。ここも先程と同様です。
$ echo "track type=bed name=SFPQ_eCLIP_rep1_ENCFF960OTE_peak description=SFPQ_eCLIP_rep1_ENCFF960OTE_peak visibility=2 maxHeightPixels=40:40:20 color=255,0,0" > ./header_tmp.txt $ cat header_tmp.txt SFPQ_eCLIP_rep1_ENCFF960OTE_for_UCSC.bed > SFPQ_eCLIP_rep1_ENCFF960OTE_for_UCSC_header_plus.bed $ bzip2 -c SFPQ_eCLIP_rep1_ENCFF960OTE_for_UCSC_header_plus.bed > SFPQ_eCLIP_rep1_ENCFF960OTE_for_UCSC_header_plus.bed.bz2
データのアップロード
UCSC genome browserにデータをアップロードします。
ホーム画面から、[MyData] -> [Custom Tracks]をクリックします。
Add Custom Tracksというページに移動するので、ファイルを選択
をクリックし、アップロードするファイルを選択し、submit
ボタンをクリックします。
アップロードが完了すると、Manage Custom Tracksというページに移動します。右端にあるgo
ボタンをクリックすると、データをみることができます。
可視化したデータをみてみると、以下のようになります。
ざっくりとですが、以上になります。
HoloLensを使ってPDBに登録されているタンパク質の立体構造を現実空間に投影する
概要
HoloLens向けのアプリをビルドして実機にデプロイするところをまとめてみました。
Microsoftが用意しているHello World的なチュートリアルをなぞってみます。
Holograms 100
ただのCubeを表示させるだけだと味気ないので、
PDBからタンパク質の立体構造のデータを取ってきて、PyMolとBlenderで適当に加工した3Dデータを表示させます。
動画をキャプチャすると見た目がしょぼく見えてしまいますがこんな感じになります。
このように、HoloLensを使えば、空間を認識してマーカーなしでその場にオブジェクトを置くことができます。
必要なもの
HoloLensのアプリ開発環境
- Visual Studio 2015 Update 3
http://dev.windows.com/downloads - HoloLens Emulator
http://go.microsoft.com/fwlink/?LinkID=823018 - Unity 5.5
https://store.unity.com/download
PDBデータの加工
- Anaconda3 (Python3.5, Scipy)
- PyMOL
- Blender 2.78a
https://www.blender.org/download/
内容
1. HoloLensのアプリ開発に必要なツールのインストール
詳しくは、MicrosoftのInstall the toolsを参照のこと。
https://developer.microsoft.com/ja-jp/windows/holographic/install_the_tools
- Visual Studio 2015 Update 3
http://dev.windows.com/downloads - HoloLens Emulator
http://go.microsoft.com/fwlink/?LinkID=823018 - Unity 5.5
https://store.unity.com/download
をそれぞれインストールする。
インストール時の留意点
Visual Studio 2015 Update 3
Tools (1.4)
とWindows 10 SDK (10.0.10586)
を併せてインストールする必要がある。
HoloLens Emulator
Hyper-V
の機能をあらかじめONにしておく必要がある。
方法は以下の通り。
1. CortanaでHyper-V
と検索。
2. Windowsの機能の無効化と有効化(コントロールパネル)をクリック。
3. Hyper-V
の項目にチェックを付け、パソコンを再起動
Unity
Windows Store .NET Scripting BackEnd
を併せてインストールする。
2. PDBデータの加工に必要なツールのインストール
Anacondaのインストール
PyMolがPythonで動作するソフトであるため、Python3をあらかじめパソコンにインストールしておく。
Download Anaconda Now! | Continuum
Python2.7
でも構わないが、今回はPython3.5
のほうをダウンロードしてインストールした。
PyMOLのインストール
下記の記事に紹介されているので詳しくはそちらを参照のこと。
PyMOL インストール
1. 非公式のPythonパッケージのインストーラーを配布しているサイトへ行き、
Python Extension Packages for Windows - Christoph Gohlke
- Pmw-2.0.1-py3-none-any.whl
- pymol-1.8.5.0-cp35-cp35m-win_amd64.whl
- pymol_launcher-1.0-cp35-cp35m-win_amd64.whl
をそれぞれダウンロードする。
2. コマンドプロンプトを開き、先ほどダウンロードしたファイルからPyMOLをインストールする。
$ pip install Pmw-2.0.1-py3-none-any.whl $ pip install pymol-1.8.5.0-cp35-cp35m-win_amd64.whl $ pip install pymol_launcher-1.0-cp35-cp35m-win_amd64.whl
Anacondaをインストールしていれば、Scipy
とWheel
がインストールされているはずですが、万が一上記のライブラリがインストールされていない場合は下記のように自身でインストールする。
Anaconda Navigatorを起動して、Environmentsをクリックし、Allを選択して、それぞれのライブラリを検索する。 インストールしたいライブラリにチェックを付け、右下のApplyボタンをクリックして該当するライブラリをインストールする。
3. PyMOLの実行ファイルのショートカットを作っておく。
PyMOLはC:\Users\自分のユーザー名\Anaconda3\PyMOL.exe
にあるので、ショートカットを作っておいてデスクトップなどに置いておくと便利です。
Blenderのインストール
下記のURLからインストーラーをダウンロードし、自身の環境にインストールする。
https://www.blender.org/download/
3. Unityに取り込ませるタンパク質の立体構造のデータを用意する
PDBのサイトから適当なタンパク質の立体構造をダウンロード
下記のURLからPDBのサイトへアクセスします。
RCSB Protein Data Bank - RCSB PDB
PDBのサイトから、好きなタンパク質の名前で検索をかけ、適当なタンパク質のデータを選択する。
各立体構造のページに飛んだら、右上のDownload Files
ボタンをクリックして、PDB Format
のファイルをダウンロードする。
ちなみ、今回使用したデータはこちら。
https://files.rcsb.org/download/3Q0L.pdb
PyMOL上でリボン構造とSurface構造を保存する
PyMOLをWindows10上で起動し、[File]->[Open]から先ほどダウンロードしてきたPDBファイルを開く。
開いたデータ(ここでは「3qOl」)の欄の横にあるH
ボタンをクリックし、メニューの先頭のeverything
を選択する。
すると、右図のように立体構造の表示がクリアされる。
次に、S
ボタンをクリックし、メニューのcartoon
を選択する。
すると、リボン構造としてタンパク質の立体構造が表現される。
これをBlenderで取り込めるVRML2
ファイルとして保存する。
Surface構造も同様の手順で取得する。
まず、先ほどと同様にH
ボタンをクリックし、メニューの先頭のeverything
をクリックすることでリボン構造を消す。
次に、S
ボタンをクリックし、メニューの先頭のsurface
を選択する。
すると今度は、Surface構造としてタンパク質の立体構造が表示される。
これをBlenderで取り込めるVRML2
ファイルとして保存する。
Blenderで3Dデータを加工する
Blenderを起動する。起動した初期画面にはCubeがあらかじめ設置されているので、まずはこれを消す。 右クリックでオブジェクトを選択して、キーボードの「Delete」ボタンを押すと、ウインドウ上に「Delete」するか求められるので「Delete」をクリックしてCubeを消す。
BlenderがVRML2
ファイルを取り込めるかどうか確認する。
[File]->[User Preferences]をクリックして、Blender User Preferences
を開く。
Add-onsのタブをクリックし、VAML2を検索する。Web3D X3D/VAML2 format
の項目にチェックが付いていない場合は、チェックを付けて左下のSave User Settings
をクリックする。
次に、Blenderにタンパク質の立体構造のデータ(リボン構造とSurface構造)を取り込ませる。
[File]->[Import]->[X3D Extensible 3D(.x3d/.wrl)]をクリックする。そして、先ほどPyMOLで用意したリボン構造のほうのVAML2ファイルを選択しBlenderにインポートする。
- 取り込んだリボン構造に対して右クリックし、オブジェクトを選択する。
- キーボードの
Tab
キーを押して、Object Mode
からEdit Mode
に切り替える。 - 右端のカラムから
Material
ボタンを押す。 +
ボタンをクリックして、新しいMaterialMaterial.001
を作成する。- 新しく作ったMaterialを選択して、
Assign
ボタンをクリックする。 - Diffuseの項目で適当な色を選択する。
Tab
キーを押して、Object Mode
に戻る。- もとからあった
Shape
を選択し、-
ボタンをクリックして消す。
結果として、このようにリボン構造に色がつくはずです。
同じ要領でSurface構造をBlenderにインポートする。
- 先程と同様に、右クリックでSurface構造のみを選択した状態で、
Tab
キーを押しEdit Mode
に切り替える。 - 右端のカラムから
Material
ボタンを押す。 - Diffuseの項目で適当な色を選択する。
Transparency
にチェックを付ける(リボン構造と異なり、Surface構造は透明になるように設定します)。Alpha
値を0.3
に設定する(透明度を設定。0に近いほど透明度が増す)。Tab
キーを押し、Object Mode
に切り替える。- キーボードの
A
キーを押し、オブジェクトを全選択。 S
キーを押し、マウスを動かしてオブジェクトのサイズが小さくなるように調節する。 もとのサイズだと保存したファイルが大きくなるので、できるだけ容量を抑えるためにサイズを小さくします。
最後に、完成したオブジェクトを、UnityにインポートできるFBX
ファイルとしてエクスポートする。
[File]->[Export]->[FBX(.fbx)]をクリックし、任意の場所・名前で保存。
4. Unity上でHoloLensアプリをつくる
新しいプロジェクトの作成
Unityを起動して、右上のNEW
をクリック。プロジェクト名をPDBViewer
として、Create project
ボタンをクリックして、新しいプロジェクトを作成する。
タンパク質の立体構造のデータをインポート
上記で作成したタンパク質の立体構造のFBX
ファイルをUnityにインポートする。
方法は簡単で、Project
タブのAssetsの直下に、該当のFBXファイルをドラッグ&ドロップでインポートする。
[Ctrl]+[S]で新しいSceneを適当な名前で保存する。今回はPDBViewer
という名前で保存した。
HoloLens向けにMain cameraの設定
次に、HoloLens向けにMain cameraの設定(Inspector)を変更する。
- TransformのPositionを
(X: 0, Y: 1, Z: -10)
から(X: 0, Y: 0, Z: 0)
に変更。 - CamemraのClear Flagsを
Skybox
からSolid color
に変更。 - CameraのBackgroundのRGBA値を
(0, 0, 0, 0)
に変更。 - CameraのNear Clip Planeを
0.3
から0.85
に変更(85 cmより近い地点ではオブジェクトの表示を消す)。
タンパク質の立体構造のデータをSceneに追加
Assetsから、Hierarchy
にタンパク質の立体構造のデータをドラッグ&ドロップで配置する。
TranformでPositionを(X: 0, Y: 0, Z: 0.5)
、Scaleを(X: 0.1, Y: 0.1, Z: 0.1)
に設定する。
HoloLens向けにビルドするための設定変更
1. Unity Performance settings
[Edit]->[Project Settings]->[Quality]をクリックして、Qualityの設定画面を出す。
Default
の横にある下三角ボタンをクリックし、Fastest
を選択する。
QualityをFantastic
からFastest
にしておかないとFPSの低下を招き、結果としてHoloLens上での動作がもたつく。
2. Unity Build settings
HoloLens向けのビルドの設定を行う。
まず、[File]->[Build Setting…]をクリックして、ビルドの設定画面を出す。
- SDKを
Universal 10
に変更。 - UWP Build Typeを
D3D
に変更。 Unity C# Projects
にチェックを付ける。
最後に、Windows Storeを選択して、左下のSwitch Platform
をクリックする。
すると、UnityのマークがWindows Storeに移動し、Build対象が切り替わったことが確認できます。
Player Settings...
をクリックし、Settings for Windows StoreのOther settingsのVirtual Reality Supported
にチェックを付ける。
Virtual Reality SDKsの項目で、Windows Holographic
を選択する。
先ほど行ったビルドの設定を先にしておかないと、Virtual Reality Supported
の設定でWindows Holographic
が選択肢に出てこないので注意。
3. ビルドを行う
[File]->[Build Settings…]をクリックし、ビルドの設定画面を開き、Build
ボタンをクリック。
どこにアプリをビルドするか聞かれるので、App
フォルダを新規作成してそこにアプリをビルドする。
5. Visual StudioからHoloLensへアプリをデプロイする
ビルドが終わるとプロジェクトのウインドウが開くので、AppフォルダからPDBViewer.sln
をクリックする。
Visual Studioが起動する。Package.appxmanifest
を右クリックして、View Code
をクリックする。
TargetDeviceFamily
の名前のところをWindows.Universal
からWindows.Holographic
に変更し、ファイルを保存。
HoloLensをUSB経由でパソコンに接続し、Visual Studio上で以下のように設定する。
Debug
をRelease
に変更。ARM
をx86
に変更。Local Machine
をDevice
に変更。
[Debug]->[Start Without Debugging]をクリックし、HoloLensにアプリをデプロイする。
以上になります。よくよく見ると今回プログラミングを一切していません…。
参考
詳しくはこちらを参照のこと。
Windows10でNGS解析をやってみる - Bash on Ubuntu on Windows × Bioconda
概要
次世代シーケンス(NGS)のデータ解析というとLinux OS上で作業するのが一般的で、Windowsで行うためにはVMwareやVirtualBoxでLinuxの仮想環境を構築する必要がありました。
ただ、この方法だと環境構築が面倒で、場合によってはうまく動かないケースがありました。また、LinuxとWindowsの両方のシステムを動かすことになるので、パソコンへの負荷が大きいのも欠点でした。
しかし、2016年8月3日に配信開始された大型アップデート「Windows 10 Anniversary Update」で状況が変わりました。
このアップデートでLinux(Ubuntu)がWindows上で使用可能になる「Windows Subsystem for Linux (WSL)」が追加され、「Bash on Ubuntu on Windows」と呼ばれるBashがWindows上でも使えるようになりました。
今回はこのBash on Ubuntu on Windows
とBioconda
を組み合わせることで、NGSのデータ解析を行ってみたいと思います。
内容
1. Bash on Ubuntu on Windowsのインストール
下記のブログ記事ですでに紹介されているのでそちらを参照のこと。
インストール後、以下のようにメニュー画面からBash on Ubuntu on Windows
を呼び出すことができるようになります。
起動するとこんな感じです。
2. Biocondaのインストール
以前のブログで記事ですでに紹介したのでそちらを参照のこと。
Bash on Ubuntu on Windows
上でもminiconda
やBioconda
をインストールすることができます。
Bash on Ubuntu on WindowsからWindows内のフォルダにアクセスする
Bash on Ubuntu on Windows上からだと、WindowsのCドライブなどが/mnt
上にマウントされているように見えます。
$ cd /mnt/c $ cd /mnt/d
などとそれぞれのドライブにアクセスできます。
もちろん、Windowsのシステムファイルに直接アクセスすることはできませんが、ファイル操作などは行うことができます。 そのため、LinuxのサブシステムとWindows10の間でのファイル共有も容易です。
$ touch test.txt $ mv test.txt /mnt/d
RNA-seqのデータをいじってみる
Biocondaを使って、sra-tools
とfastqc
をインストールします。
$ conda install sra-tools $ conda install fastqc
今回、下記のRNA-seqのSRAファイルを使用する。
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSM1606099
Bashのホームディレクトリではなく、Dドライブ上にデータを配置するため、Dドライブに移動する。
cd /mnt/d
Bash on Ubuntu on Windows上でSRAファイルをダウンロードする。
$ wget ftp://ftp-trace.ncbi.nlm.nih.gov/sra/sra-instant/reads/ByExp/sra/SRX/SRX870/SRX870796/SRR1795409/SRR1795409.sra
SRAファイルを解凍してFASTQファイルを用意する。
fastq-dump SRR1795409.sra
FastQCを使って、FASTQファイルの中身をチェックする。
fastqc SRR1795409.fastq
Dドライブ上にFastQCのデータが保存されているので、Windows上からアクセスしてファイルの中身を見ることができます。
こんな感じ。
ガッツリ解析しようと思ったらスパコンを利用するほうが良いと思いますが、Windows上にテスト環境を作って練習するには手軽かもしれません。
NCBI GEOからまとめてSRAファイルを取得する
概要
NCBI GEOからまとめてSRAファイルを取得したい。
準備
Entrez Direct
のコンパイル済みのバイナリファイルをダウンロードして、パスを通すだけ。
* Entrez Direct
ftp://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/
https://www.ncbi.nlm.nih.gov/news/02-06-2014-entrez-direct-released/
インストール例
$ wget ftp://ftp.ncbi.nlm.nih.gov/entrez/entrezdirect/edirect.tar.gz $ tar zxvf edirect.tar.gz $ echo 'export PATH=~/softwareedirect:${PATH}' >> ~/.bashrc
もしくは、biocondaでインストールする。
$ conda install entrez-direct
実行例
下記のコマンドを実行。
esearch
の-query
オプションに適当なBioProjectのIDを指定する。
$ esearch -db sra -query PRJNA328218 | efetch --format runinfo | cut -d ',' -f 10 > srafile.txt $ wget -i srafile.txt
SRAファイルのダウンロード先を取得して、wgetでまとめてダウンロードする。
esearch
の-db
オプションからさまざまなデータベース(SRA、PubMedなど)を指定でき、-query
オプションでIDなどをもとに検索を行うことができる。
efetch
で検索した項目の情報を取得することができる。--formatの指定は、runinfo
とnative
の2種類ありruninfo
ではカンマ区切りで各サンプルの情報が出力される(すべての情報が得られるわけではなく概要のみ?)。
一方、native
を指定すると、XMLファイル形式ですべての情報が出力される(すべての情報が揃っているが、中身を調べるにはXMLをパースして整理する必要がある)。
SRAのファイル名にサンプル名を加えたい
SRAファイルをダウンロードした後、SRA ID
とサンプル名
を併記したファイル名にリネームする。
シェルスクリプト
以下のようにスクリプトを実行する。
./sra_download.sh <任意のBioProjectのID>
実行に必要なPythonスクリプト
parse_sra_xml.py (1.4 kB)
rename_sra_files.py (352 B)
使用する場合、
filepath="/path/to/python_scripts"
の部分に実行に必要なPythonスクリプトを置いた場所を指定する。
やっていることは単純で、
1. 特定のProject IDに含まれる複数のサンプル(SRAファイル)をダウンロード&XMLファイルを取得する。
2. XMLファイルをパースして、SRAファイル名とサンプル名の対応表を作る。
3. SRA IDとサンプル名が併記されたファイル名にリネームする。
という作業をPythonのスクリプトも使いながら(スマートじゃないですが)行っています。
実行すると、以下のようにSRAファイルが一括ダウンロードできます。
参考
- Tutorial: How to download raw sequence data from GEO/SRA
https://www.biostars.org/p/111040/ - NCBI SRA から FASTQ をダウンロードする方法
http://bi.biopapyrus.net/transcriptome/rnaseq-data/download-from-sra.html - prefetch すらっと落とす SRA
http://blog.amelieff.jp/?eid=231191
biocondaを利用してNGS関連のソフトウェアを一括でインストールする
概要
biocondaと呼ばれるパッケージマネージャーを用いて、バイオインフォマティクス関連のソフトウェアのインストールから、バージョン管理までを行うことができます。
バイオインフォマティクスの解析環境を用意しようと思っているなら、自力でソフトウェアをインストールするのではなく、biocondaを利用したほうが圧倒的に楽です(MEMEやCircosなどインストールに手間取るソフトウェアが一行のコマンドでインストールできちゃう)。
biocondaには現在、15,000以上のBiology関連のソフトウェア・パッケージが登録されており、おおむね、どのソフトウェアに関しても最新バージョンがインストールできる環境が整備されているみたいです。
ただし、パッケージが充実しているのはLinux版で、MacOS版はbiocondaの(インストールの)レシピが充実していない印象です。 なので、biocondaを導入する対象はLinux OSであることが望ましいと思います。
内容
基本的に、こちらのbioconda公式サイトにインストール方法が記述されているので、詳しくはこちらを参照のこと。
https://bioconda.github.io/
インストール方法
1. minicondaをインストールする(Anacondaというのもありますが、minicondaのほうはその最小構成版になっています。余計なパッケージをはじめから入れさせないためにminicondaを選択したほうがよいかもしれません)。
Python3がおすすめと書いてあるが、NGS関連のソフトウェアはpython2でコーディングされているものも多いため、python2の方を選択する。
http://conda.pydata.org/miniconda.html
2. 上記のURLからシェルスクリプトをダウンロードし、実行する。
ここでは、64bit版のPython2の環境を選択してインストールしています。
$ wget https://repo.continuum.io/miniconda/Miniconda2-latest-Linux-x86_64.sh $ bash Miniconda2-latest-Linux-x86_64.sh
いろいろ聞かれるが、基本的にEnter
を押すか、Yes
を入力すればいい。
最終的に、カレントディレクトリにminiconda2
というディレクトリができており、.bashrc
に以下の文が追加されているはずです。
> .bashrc # added by Miniconda2 4.2.12 installer export PATH="/home/user/miniconda2/bin:$PATH"
最後の質問で、No
を選択すると、上記の文が追加されません。
その場合、自身で上記の文を記述して、miniconda2/bin
にパスを通す必要があります。
3. インストール終了後、.bashrc
を再読み込み。
$ source ~/.bashrc
4. condaが正常にインストールされたかどうかチェック。 エラーが返ってこなければOK。
$ conda -h
5. biocondaのチャンネルを追加。biocondaチャンネルが最上位に来るように設定。
$ conda config --add channels conda-forge $ conda config --add channels defaults $ conda config --add channels r $ conda config --add channels bioconda
6.任意のパッケージ・ソフトウェアをインストール。
$ conda install cutadapt $ conda install meme $ conda install fastqc $ conda install fastx_toolkit $ conda install bowtie $ conda install bowtie2 $ conda install prinseq $ conda install circos $ conda install tophat $ conda install cufflinks $ conda install samtools $ conda install bedtools $ conda install sra-tools $ conda install trimmomatic $ conda install paralyzer $ conda install star $ conda install rsem $ conda install bwa $ conda install picard $ conda install subread $ conda install entrez-direct $ conda install art $ conda install macs2 $ conda install htseq $ conda install rseqc
biocondaに登録されているパッケージ群のリストは下記のサイトからチェックできる。
https://bioconda.github.io/recipes.html
インストール時の注意点
biocondaをインストールする際に、PYTHONPATH
が.bashrc
に記載されていると、うまくインストールされない。
インストールチェックを行うと、エラーが発生する。
$ conda -h Traceback (most recent call last): File "/home/akimitsu/miniconda3/bin/conda", line 4, in <module> import conda.cli ImportError: No module named 'conda'
対処法としては、.bashrc
中でPYTHONPATH
をExportしている箇所をコメントアウトする。ターミナルを再起動・再ログインしてから、再インストールする。
# コメントアウト # export PYTHONHOME=/usr/local/package/python2.7/current # export PYTHONPATH=/home/user/software/python_path27/lib/python2.7/site-packages # export PATH=/home/user/software/python_path27/bin:${PYTHONHOME}/bin:${PATH} # export LD_LIBRARY_PATH=/home/user/software/python_path27/lib:${PYTHONHOME}/lib:${LD_LIBRARY_PATH}
biocondaに登録されているパッケージのインストールに関する注意点
python2に依存しているパッケージは、python3環境下ではインストールできない。 そのため、以下で説明するPython2用の仮想環境を構築するか、python2向けのminicondaを再インストールする必要がある。
$ conda install macs2 Fetching package metadata ............. Solving package specifications: .... UnsatisfiableError: The following specifications were found to be in conflict: - macs2 - python 3.5* Use "conda info <package>" to see the dependencies for each package.
python2に依存しているパッケージ例
conda install macs2 conda install htseq conda install rseqc
biocondaの環境をチェック
biocondaでインストールしたパッケージは、conda list
で確認できる。
$ conda list # packages in environment at /home/akimitsu/miniconda2: # art 3.19.15 1 bioconda bcftools 1.3.1 1 bioconda bedtools 2.26.0 0 bioconda bowtie 1.2.0 py27_0 bioconda bowtie2 2.3.0 py27_0 bioconda bwa 0.7.15 0 bioconda bx-python 0.7.3 py27_0 bioconda bzip2 1.0.6 3 cairo 1.14.6 0 cffi 1.9.1 py27_0 circos 0.69.2 0 bioconda ...
python2向けのbiocondaでpython3を使う
python3を動かすための仮想環境を構築することができる。
以下では、py35
という名前のpython3.5がインストールされた仮想環境を作っている。
$ conda create --name py35 python=3.5
--name
で仮想環境の名前を指定、
python=
でPythonのバージョンを指定する。
condaでは現在、Python 2.7, 3.4, 3.5に対応している。
Managing Python — Conda documentation
登録した環境設定を確認する。
$ conda info --envs # conda environments: # py35 /home/your_name/miniconda2/envs/py35 root * /home/your_name/miniconda2
環境の切り替え。
# py35に入る source activate py35 # py35から抜ける source deactivate py35
これを利用して、プロジェクトごとに環境を用意して特定のパッケージ・バージョンで解析を行うこともできます。
CRANに登録されているRパッケージをインストールする
biocondaに登録されていないPythonやRのパッケージは、pipやinstall.packages()で正常にインストールできないことがある。 とりあえず、Rのパッケージは以下の通りに入力するとインストールされるみたいです。
$ R > install.packages('任意のパッケージ', repos='http://cran.us.r-project.org')
参照
- bioconda
https://bioconda.github.io/