ラベル Python の投稿を表示しています。 すべての投稿を表示
ラベル Python の投稿を表示しています。 すべての投稿を表示

2013年7月27日土曜日

RNA-seqデータの解析パイプラインを作ろう(データの視覚化)

RNA-seqのデータ(fastqファイル)をBowtieやTopHatなどでマッピングした後、bamファイルやsamファイル(場合によってはbedファイル)として出力されます。

そのデータを用いて、Cufflinksなどのリードの集計ソフトで定量値(FRKM値やRPKM値)を算出し、各転写産物の発現量を見積もると同時に、実際にゲノム上にマッピングされたリードを確認したいケースもあると思います。

そこで、今回はbam / sam / bedファイルから視覚化データを作る方法例を示したいと思います。個人でマッピングされたリードを確認する際にはIGVが便利だと思いますが、他人とデータを共有する場合では、UCSC genome browser上で可視化させたほうが都合が良いと思います。

〈スクリプトのダウンロード〉
今回は、github上で公開されているPythonスクリプトを利用します。このスクリプトはメモリの消費が少なく、複数のファイル形式(sam / bed / bowtie)に対応しているので、使い勝手が良く個人的におすすめです。

(1)rdbio-scriptsのURLに行き、右下にある「Download ZIP」をクリックし、bed2wig.pyを含む全スクリプトをダウンロードする。もしくはbed2wig.pyのURLへ移動し、スクリプトをコピペして保存する。

bed2wig.pyのスクリプト
https://github.com/daler/rdbio-scripts/blob/master/sequenceFiles/bed2wig.py
rdbio-scripts
https://github.com/daler/rdbio-scripts

〈スクリプト使用例〉
基本的な使い方は、
bed2wig.py -h
でヘルプを見ることができます。
ここでは、bamファイルをwigファイルに変換する方法について紹介したいと思います。

(1)bamファイルをソートし、samファイルに変換する。(TopHatから出力されるbamファイルは予めソートされていたと思うのでこの作業は必要ないかも。)
samtools sort accepted_hits.bam accepted_hits_sorted.bam
samtools view -F 0x0004 accepted_hits_sorted.bam > accepted_hits_sorted.sam
⇒「-F」オプションで0x0004を指定することでマッピングされなかったリードを排除しています。

(2)samファイルからwigファイルに変換する。
bed2wig.py -i accepted_hits_sorted.sam --type sam -o  accepted_hits_sorted.wig --verbose
⇒「-i」オプションでインプットファイルを指定します。「-type」オプションでインプットファイルのファイル形式を指定します。「-o」オプションでアウトプットファイルを指定します。「--verbose」オプションで進捗状況を出力させます。

(3)wigファイルの1行目を書き換える。
lessコマンドやviでファイルの中身をみてみると、1行目は「track type=wiggle_0 alwaysZero=on」となっています。このままだと使いにくいので、下記のように、適当に書き換えておきます。(もちろん、bed2wig.pyを直接書き換えても良いと思います。たぶん、269行目を書き換えればOK.。)

track name='CTRL_TopHat' description='CTRL_TopHat' type='wiggle_0' color=255,128,128 yLineMark=0.0 yLineOnOff=on visibility=2 maxHeightPixels=40:40:20 priority=1

「track name」と「description」でデータの名称を各ファイルごとに書き換えておきましょう。UCSC genome browserにアップロードする場合は、ここで指定した名称が以前アップロードしたデータと重複していると、アップロードした時に前のデータが消えてしまうので注意が必要です。

(4)ファイルを圧縮する。
UCSC genome browserではgzip (.gz), compress (.Z), or bzip2 (.bz2)の3種類の圧縮形式に対応しており、ファイルが大きくてアップロードに時間がかかりそうな場合、ファイルを圧縮しておきましょう。
bzip2 -c accepted_hits_sorted.wig > DATA.wig.bz2
⇒「-c」オプションを指定し、標準出力でファイルを出力。

〈UCSC genome browserへのアップロード〉
(1)Homeへ移動し、左側にある「Genome Browser」をクリック。
http://genome.ucsc.edu/index.html

(2)「manage custom track」ボタンをクリック。(以前にアップロードしたファイルが残っていると、「Manage Custom Tracks」ページに移動するので、その場合は「add Custom Tracks」ボタンをクリック。)

(3)「add Custom Tracks」ページに移動するので、「Paste URLs or data:」の「Or upload:」の「ファイルを選択」をクリックし、アップロードしたいwigファイルを指定する。

(4)「Submit」ボタンをクリック。ファイルがアップロードされるので、アップロードが完了するまで放置する。(アップロードしている間、ページ移動をしないこと。)

(5)アップロード完了後、「Manage Custom Tracks」ページに移動します。アップロードしたファイルは「go to Genome browser」ボタンをクリックし、Genome browser上で確認することができます。

〈Custom tracksの維持・構成の保存〉
Custom trackをアップロードしても、一週間ぐらいすると設定がリセットされて見れなくなってしまいます。そこで、Custom tracksおよびその他のTrackの構成を保存するために、UCSC genome browserが提供している「Sesson」機能を利用しましょう。

「Sesson」では、アカウントを取得することで、Custom trackの維持やTrackの構成の保存を行うことができます。また、アップロードしたCustom trackを他人が見れるように設定することもできます。

(1)Homeへ移動し、上の「My Data」→「Sessions」をクリック。
http://genome.ucsc.edu/index.html

(2)「Creat an account」をクリックし、必要なデータを入力しアカウントを作成する。

(3)Sessionにログインした後、「Save Settings」で「Save current settings as named session:」にSessionの名称を入力。

自分のアカウントにログインしていない他人にデータが見られるようにしたい場合は「allow this session to be loaded by others」のボックスにチェックをつける。

そして、「Submit」ボタンをクリックし、Sessionの保存を行う。(もちろん、Custom trackの追加やgenome browser上の構成を自分なりにアレンジして保存したい組み合わせを事前にgenome browser上で組んでおく。)

⇒Sessionを上書きしたい場合は、上書きしたいSession名を入力した後、「Submit」ボタンをクリック。上書きしてよいか確認を求められるので、「Yes」ボタンをクリック。

2013年3月23日土曜日

Pythonのインストール方法(Linux OS)

Pythonのインストール方法のメモ。

(1)Pythonのホームページに移動。

(2)左にある「DOWNROAD」をクリック。

(3)Python 2.7.3 bzipped source tarball (for Linux, Unix or Mac OS X, more compressed)を選択し、ファイルをダウンロードする。

(4)ダウンロードしたファイルの解凍
$ tar xjvf Python-2.7.3.tar.bz2

(5)解凍したディレクトリの場所に移動。
$ cd Python-2.7.3/

(6)以下のコマンドに順番に実行。
$ ./configure

$ make

$ sudo make install

(7)Pythonがインストールされているか確認。
$ python -h
usage: python [option] ... [-c cmd | -m mod | file | -] [arg] ...
Options and arguments (and corresponding environment variables):
-B     : don't write .py[co] files on import; also PYTHONDONTWRITEBYTECODE=x
-c cmd : program passed in as string (terminates option list)
-d     : debug output from parser; also PYTHONDEBUG=x
-E     : ignore PYTHON* environment variables (such as PYTHONPATH)
-h     : print this help message and exit (also --help)
-i     : inspect interactively after running script; forces a prompt even
         if stdin does not appear to be a terminal; also PYTHONINSPECT=x
-m mod : run library module as a script (terminates option list)
-O     : optimize generated bytecode slightly; also PYTHONOPTIMIZE=x
-OO    : remove doc-strings in addition to the -O optimizations
-R     : use a pseudo-random salt to make hash() values of various types be
         unpredictable between separate invocations of the interpreter, as
         a defense against denial-of-service attacks
-Q arg : division options: -Qold (default), -Qwarn, -Qwarnall, -Qnew
-s     : don't add user site directory to sys.path; also PYTHONNOUSERSITE
-S     : don't imply 'import site' on initialization
-t     : issue warnings about inconsistent tab usage (-tt: issue errors)
-u     : unbuffered binary stdout and stderr; also PYTHONUNBUFFERED=x
         see man page for details on internal buffering relating to '-u'
-v     : verbose (trace import statements); also PYTHONVERBOSE=x
         can be supplied multiple times to increase verbosity
-V     : print the Python version number and exit (also --version)
-W arg : warning control; arg is action:message:category:module:lineno
         also PYTHONWARNINGS=arg
-x     : skip first line of source, allowing use of non-Unix forms of #!cmd
-3     : warn about Python 3.x incompatibilities that 2to3 cannot trivially fix
file   : program read from script file
-      : program read from stdin (default; interactive mode if a tty)
arg ...: arguments passed to program in sys.argv[1:]

Other environment variables:
PYTHONSTARTUP: file executed on interactive startup (no default)
PYTHONPATH   : ':'-separated list of directories prefixed to the
               default module search path.  The result is sys.path.
PYTHONHOME   : alternate <prefix> directory (or <prefix>:<exec_prefix>).
               The default module search path uses <prefix>/pythonX.X.
PYTHONCASEOK : ignore case in 'import' statements (Windows).
PYTHONIOENCODING: Encoding[:errors] used for stdin/stdout/stderr.
PYTHONHASHSEED: if this variable is set to 'random', the effect is the same
   as specifying the -R option: a random value is used to seed the hashes of
   str, bytes and datetime objects.  It can also be set to an integer
   in the range [0,4294967295] to get hash values with a predictable seed.

$ python --version
Python 2.7.3

以上。

2013年3月2日土曜日

Pythonモジュールのインストール方法(Linux OS)

Pythonモジュールのデータをダウンロードすると、「seyup.py」というPythonのスクリプトが入っています。
そこで、このスクリプトを利用して以下のコマンドを打ち込むことで、モジュールをインストールすることができます。

$ python setup.py install

基本的にこのコマンドを打てばインストールできるはずです。