Research Tips (Only Japanese)

DDBJのスパコン上でDe novoアセンブリ解析をするための方法についての備忘録（随時更新）

目次

事前準備

1.  DDBJスパコンのウェブサイトにて、ユーザーアカウント発行基準及び利用規定を確認

2.  SSH公開鍵・秘密鍵の生成、及びスパコンのアカウント新規登録を行う（通常は一般解析区画）（方法はこちらを参照）

3. 申請者（責任者）宛に誓約書のPDFファイルが添付されたメールが届くので、誓約書に電子署名した上で返信する（方法はこちら）

4. アカウント登録証が申請時に登録した住所に郵送で届く

5. アカウント登録証に書いてあるアカウント名と初期パスワードを確認する（このパスワードはアカウント申請システムのマイページにログインするときに使う）

6. アカウント名と（SSH公開鍵・秘密鍵を作成したときに設定した）パスワードを用いて、スパコンにログインする

DDBJスパコンへのログイン方法（Macの場合）

*DDBJスパコンのウェブサイトを参照

1. ターミナルを開き、「ssh (アカウント名)@(ゲートウェイノード名)」と入力（一般解析区画には2種類のゲートウェイノードがあるが、通常は@gw.ddbj.nig.ac.jpでいい）。

例： ssh usrname@gw.ddbj.nig.ac.jp

2. "Enter passphrase for key ..."と表示されるので、SSH公開鍵・秘密鍵を作成した際に設定したパスフレーズを入力

3. 以下のような文章が表示されるはずなので、「qlogin」と入力し、自分のアカウントのホームディレクトリにログインする

---------------------------------------------------------------------

Thank you for using NIG supercomputer system.

This is the gateway node, do not run program here.

Please use 'qlogin' to login to a login node.

---------------------------------------------------------------------

（メモ）ログインノード上でのメモリ最大利用量を指定してqloginする場合は以下のように入力

   qlogin -l s_vmem=50G -l mem_req=50G　　（50Gは50GBのこと）

4. 「exit」または「logout」コマンドを2回入力すると、DDBJのスパコンから完全にログアウトする（1回目のexitでqloginからログアウト、2回目でDDBJスパコンからログアウトしている）

De novoアセンブリと発現量推定

　次世代シーケンサー（NGS: next generation sequencing）から出力された、ゲノム未解読の生物の全mRNAの塩基配列の大量の断片を、Trinityと呼ばれるプログラムを用いてde novoでアセンブルし、その結果を元に、RSEMというソフトウェアツールを用いて各mRNA配列の発現量を推定する。

（注）ペアエンドリード（Pair-end read）を読んだサンプルで、かつトリミング済みのFASTQファイルを受け取ってから始める場合を以下に記す

①Trinityのインストール

 　Trinity自体はDDBJスパコンのSingularityコンテナの中に入ってはいるため、TrinityだけするのであればわざわざDDBJのスパコンの自分のホームディレクトリに入れる必要はないのだが、RSEMもSingularityコンテナ内のTrinityの機能を使って実行しようとすると、なぜかエラーが出て実行できない。そのため、自分のホームディレクトリにもTrinityをインストールしておく。

1. qloginしてホームディレクトリにアクセス

2. 次のコマンドを入力し、Trinityの最新版 v2.15.2（2024年11月現在）をダウンロード

wget https://github.com/trinityrnaseq/trinityrnaseq/releases/download/Trinity-v2.15.2/trinityrnaseq-v2.15.2.FULL.tar.gz

3. 次のコードでダウンロードファイルを解凍

tar xvzf trinityrnaseq-v2.15.2.FULL.tar.gz

4. ホームディレクトリ内に"src"ディレクトリ（名前はなんでもいい）を作成し、ダウンロードファイルを解凍したものをsrcディレクトリに配置

mkdir src

mv trinityrnaseq-v2.15.2 ~/src/

5. "trinityrnaseq-v2.15.2"ディレクトリに移動

cd src/trinityrnaseq-v2.15.2

6.  lsと入力し、"trinityrnaseq-v2.15.2"ディレクトリ内に"Makefile"があることを確認

7. 「make」と入力し、コンパイルする

8. 「make plugins」と入力し、プラグインもコンパイルする

（注）trinityrnaseq-v2.15.2"ディレクトリ内にREADMEファイルがあるはずなので、それを読んで「make」でどこまでコンパイルすべきなのかは確認する

②発現量推定に必要なソフトウェア群をインストール

　発現量推定に必要なsamtools、bowtie2、bowtie2-build、rsem-calculate-expressionもホームディレクトリにインストール。

1. samtoolsをインストール

cd local #ホームディレクトリに"local"ディレクトリに移動（なければ新規作成してから移動）

mkdir bin #"local"ディレクトリ内に"bin"ディレクトリを作成

wget https://github.com/samtools/samtools/releases/download/1.21/samtools-1.21.tar.bz2 #プログラムをダウンロード

tar -xvf samtools-1.21.tar.bz2 #.tar.bz2形式であるダウンロードファイルを解凍

cd samtools-1.21 #解凍してできた"samtools-1.21"ディレクトリに移動

./configure --prefix=/home/usrname/local #samtoolsのコンパイルの準備

make -j 8 #コンパイルする

make install #インストール（これで/home/usrname/local/bin内にsamtoolsのソフトウェアツール群が入る）

cd #ホームディレクトリへ戻る

2. bowtie2とbowtie2-buildをインストール

cd local/bin #"local"ディレクトリ内の"bin"ディレクトリへ移動

wget https://sourceforge.net/projects/bowtie-bio/files/latest/download #最新版のプログラムをダウンロード

unzip download #ダウンロードした"download"という名のファイル（.zip形式）を解凍

cd bowtie2-2.5.4 #解凍してできたディレクトリへ移動（2024年11月30日時点でver. 2.5.4が最新）

make #"Makefile"という名のファイルがあるのを確認し、コンパイルする

cd #ホームディレクトリへ戻る

（注：他のウェブサイトを見るとbowtie2のソフトウェアツール群はmakeしなくてもダウンロードしただけで使えるようなことを書いてある場合があるが、2024年11月時点ではmakeしないと、解凍した"bowtie2-2.5.4"ディレクトリ内にbowtie2-buildやbowtie2-alignなどが含まれておらず、RSEM実行中にプログラムが途中でエラー吐いて止まってしまう）

3. rsem-calculate-expressionをインストール

cd local/bin #"local"ディレクトリ内の"bin"ディレクトリに移動

wget https://github.com/deweylab/RSEM/archive/v1.3.3.tar.gz #2024年11月時点での最新版のv1.3.3のRSEMをダウンロード

tar -xvf v1.3.3.tar.gz #ダウンロードファイルを解凍（これはmakeする必要なさそう）。"RSEM-1.3.3"ディレクトリができる

cd #ホームディレクトリへ戻る

③PATHを設定

　DDBJスパコンのsingularityコンテナ内にも、Trinity（とその関連ツールのalign_and_estimate_abundance.pl）、samtools、bowtie2、bowtie2-build、rsem-calculate-expressionの実行ファイルが格納されているため、デフォルト設定では「which」コマンドでこれらのソフトウェアツール群を探索すると、singularityコンテナへのパス（/usr/bin）に繋がる。しかし、最近singularityに格納されている方のソフトウェアツール群を使ってRSEMによる発現量推定を行おうとすると何故かエラーを吐いてプログラムが中断されてしまうため、自分のホームディレクトリにわざわざインストールした上記のソフトウェアツール群を用いる。

　自分のコンピューターのコンソール上で「Export PATH=...」と入力して新たなコマンド探索パスを追加しても、DDBJスパコンからログアウトしたら追加したパスはリセットされてしまう。ログアウトしても追加パスが有効になるようにするため、ホームディレクトリ内の隠しファイルである".bash_profile"にパスを直接追加する（.bashrcへパスを追加しても良い）。

1. cd #ホームディレクトリへ移動

2. nano .bash_profile #".bash_profile"の編集画面を開く

3. 以下のスクリプトをコピペ（.bash_profileに自分で設定した優先するパスに繋げる、かつ.bashrcを優先的に読み込むよう以下のスクリプトを書き込むという内容。.bashrcにパスを追加する場合は"# User specific environment and startup programs"以降の行のみコピペする）

# .bash_profile

# Get the aliases and functions

if [ -f ~/.bashrc ]; then

. ~/.bashrc

fi

# User specific environment and startup programs

PATH=/home/ymitaka/local/bin:\

/home/ymitaka/src:\

/home/ymitaka/src/util:\

/home/ymitaka/local/bin/bowtie2-2.5.4:\

/home/ymitaka/local/bin/RSEM-1.3.3:\

$PATH:$HOME/.local/bin:$HOME/bin

export PATH

4. 「Ctrl + o」で保存し、「Ctrl + x」でnanoを閉じる

5. ホームディレクトリに戻り、「source .bash_profile」と入力した後、一旦ログアウトする

6. 再度qloginして、「which」コマンドを使ってTrinity、align_and_estimate_abundance.pl、samtools、bowtie2、bowtie2-build、rsem-calculate-expressionへのパスが、自分のホームディレクトリから始まるパス（/home/usrname/ ...）が表示されるか確認する

＜追加したパスについての説明＞

/home/ymitaka/src # "Trinity"があるディレクトリへのパス

/home/ymitaka/src/util # "align_and_estimate_abundance.pl"があるディレクトリへのパス

/home/ymitaka/local/bin # "samtools"があるディレクトリへのパス

/home/ymitaka/local/bin/bowtie2-2.5.4 # "bowtie2"と"bowtie2-build"があるディレクトリへのパス

/home/ymitaka/local/bin/RSEM-1.3.3 # "rsem-calculate-expression"があるディレクトリへのパス

④Trinityの実行

　あらかじめホームディレクトリに転送した各サンプルのペアエンドリードファイル（.fq.gz形式圧縮ファイル）を使い、Forwardのリードファイル同士まとめた圧縮ファイルと、Reverseのリードファイル同士まとめた圧縮ファイルの2つを作成。

1. Forwardどうしまとめたファイル（merged_f.fq.gz）作成

cat File01_val_1.fq.gz File02_val_1.fq.gz File03_val_1.fq.gz ... File10_val_1.fq.gz > merged_f.fq.gz

2. Reverseどうしまとめたファイル（merged_r.fastq）作成

cat File01_val_2.fq.gz File02_val_2.fq.gz File03_val_2.fq.gz ... File10_val_2.fq.gz > merged_r.fq.gz

3. Trinityを実行するためのシェルスクリプト（"Trinity.sh"ファイル）を新規作成

nano Trinity.sh

4. 編集画面が出てきたら、以下のコマンドを入力

#!/bin/sh

#$ -S /bin/bash

#$ -N Trinity

#$ -l medium

#$ -cwd

#$ -l d_rt=1440:00:00

#$ -l s_rt=1440:00:00

#$ -l s_vmem=50G

#$ -l mem_req=50G

#$ -pe def_slot 5

#echo original path: $PATH

ulimit -s unlimited

echo running on `hostname`

echo starting at

date

singularity exec /usr/local/biotools/t/trinity\:2.15.2--pl5321hdcf5f25_0 Trinity --seqType fq --max_memory 50G --left merged_f.fq.gz --right merged_r.fq.gz --CPU 5 --verbose  --no_version_check

echo ending at

date

（メモ：この行以降はTrinity.shへは入力しない）

• 上記はsingularityコンテナ内に元々あるTrinity（trinity\:2.15.2--pl5321hdcf5f25_0）を実行する場合のもの（Trinityをやるだけならこれでも構わない）。自分のホームディレクトリにあるTrinityを使う場合は以下のように書く。

/home/usrname/src/Trinity --seqType fq --max_memory 50G --left merged_f.fq.gz --right merged_r.fq.gz --CPU 5 --verbose  --no_version_check

• Trinity.shでJobを投げる際、以下のパラメーターをいじればCPU数は増えてJobは早くなるが、スパコンが混んでるとJobが待ち状態になりやすい。

#$ -pe def_slot 5

--CPU 5

• メモリが足りないとコアダンプを起こしてjobが止まり、core.***というファイルができる。この場合は、以下のメモリの値を増やせば大丈夫だが、これもjobが入りにくくなる。

   #$ -l s_vmem=50G

   #$ -l mem_req=50G

• "#$ -l medium" でスパコンのどの計算ノードにJobを投稿するか指定できる。"medium" でプログラムがうまく動かない場合は"medium-ubuntu" に変更するのもアリ。

5. 保存してnanoを閉じる

6. 以下を入力してTrinityによるde novo assemblyを実行

qsub Trinity.sh

7. 以下のように表示されればJobのSubmitは完了したことになる（"XXXXXXXX"は本来は数字でJob IDを表す）

Your job XXXXXXXX (“Trinity”) has been submitted

（メモ）

Job実行中に利用可能なコマンドは以下の通り（こちらも参照）

qstat #現在のJobの状態が確認可能（stateが "r" だとJobが走っていることを示す）

less Trinity.eXXXXXXXX 出力されたエラーの内容が確認できる

less Trinity.oXXXXXXXX 出力された途中経過を確認できる

qdel XXXXXXXX Jobの一時中断

8. Trinityが完了すると、カレントディレクトリ内に"trinity_out_dir.Trinity.fasta"という名の、各contigの塩基配列がFASTA形式で書かれているファイルや"trinity_out_dir"ディレクトリなどが生成されていることを確認する。

9. Trinityによるアセンブル結果のQuality checkを出力するため、以下のコマンドを実行

singularity exec /usr/local/biotools/s/seqkit\:2.8.2--h9ee0642_0 seqkit stat -a trinity_out_dir.Trinity.fasta

（メモ：N50やGC contentなどの結果はターミナルのウィンドウ内に表示されるので、テキストファイルなどにコピペして保存）

⑤各contigの塩基配列についてORFを探索＆アミノ酸配列に翻訳して出力

　singularityコンテナ内に入っているtransdecoderを用いて、"trinity_out_dir.Trinity.fasta" 書かれている各contigの塩基配列をアミノ酸配列に翻訳する

1. 以下のコマンドを入力

singularity exec -e /usr/local/biotools/t/transdecoder:5.7.0--pl5321hdfd78af_0 TransDecoder.LongOrfs -m 100 -t trinity_out_dir.Trinity.fasta

2. ディレクトリ内に "trinity_out_dir.Trinity.fasta.transdecoder_dir" という新たなディレクトリが作られる

3. その中の "longest_orfs.pep" というファイルにアミノ酸配列に変換された全アイソフォームの配列が書かれている

⑥RSEMの実行

　自分のホームディレクトリ内で発現量解析を実行する。事前準備として、ホームディレクトリ内に元のペアエンドリードファイル（各サンプルに対応する.fq.gz形式のファイル）と、Trinityによって出力された "trinity_out_dir.Trinity.fasta" ファイルや "trinity_out_dir.Trinity.fasta.transdecoder_dir" ディレクトリなどがあることを確認する。

1. RSEMを実行するためのシェルスクリプトを新規作成

nano RSEM.sh

2. 編集可能な画面が出てきたら以下を入力（サンプルごとの発現量推定が同時並行で行われるようにする）

#!/bin/sh

#$ -S /bin/bash

#$ -N RSEM

#$ -l medium-ubuntu

#$ -cwd

#$ -l d_rt=1440:00:00

#$ -l s_rt=1440:00:00

#$ -l s_vmem=50G

#$ -l mem_req=50G

#$ -pe def_slot 5

#$ -t 1-10:1  

SAMPLES=(File01 File02 File03 ... File10)

SAMPLE="${SAMPLES[${SGE_TASK_ID}-1]}"

#echo original path: $PATH

ulimit -s unlimited

echo running on 'hostname'

echo starting at

date

/home/ymitaka/src/util/align_and_estimate_abundance.pl \

--transcripts trinity_out_dir.Trinity.fasta \

--seqType fq \

--left ${SAMPLE}_val_1.fq.gz \

--right ${SAMPLE}_val_2.fq.gz \

--est_method RSEM \

--aln_method bowtie2 \

--trinity_mode \

--prep_reference \

--coordsort_bam \

--thread_count 5 \

--output_dir ${SAMPLE}_RSEM_outdir

echo ending at

date

（メモ：この行以降はRSEM.sh内に書き込まない）

• "#$ -t 1-10:1"の1-10は繰り返し数。RNA抽出液のサンプル数の数に合わせる

• SAMPLES=("" "" ...) この中にサンプル名を入れていく。そうすると次のSAMPLE="${SAMPLES[${SGE_TASK_ID}-1]}"の中にサンプル名が格納され、"--left ${SAMPLE}_val_1.fq.gz --right ${SAMPLE}_val_1.fq.gz"の部分でサンプル名が自動で入る。${SAMPLE}以外はファイル名を共通させる必要がある。

• --output_dir ${SAMPLE}_RSEM_outdirによって、サンプル名ごとに結果がアウトプットされる。

3. 保存して、nanoを閉じる

4. 以下のように入力して、Indexファイルを作成

bowtie2-build -f trinity_out_dir.Trinity.fasta trinity_out_dir.Trinity.fasta.bowtie2

5. 以下のように入力して、各サンプルについてJobを同時並行で実行（サンプル数が10なら、qstatでJobの進捗を確認した際に10個のJobが進行しているのが見えるはず）

qsub RSEM.sh

6. RSEMが完了すると、「(サンプル名)_RSEM_outdir」という名前のディレクトリがサンプルごとに生成される。その中に、contigごとに発現量を計算した結果のファイルである "RSEM.genes.results" と、各contigのisoformごとに発現量を計算した結果のファイルである "RSEM.isoforms.results" が入っているのを確認する。

7. 以下のコマンドを入力し、全サンプルのRSEM.genes.resultsをマージして一つの表（CSVファイル）を作成

/home/usrname/local/bin/RSEM-1.3.3/rsem-generate-data-matrix \

/home/ymitaka/File01_val_1_RSEM_outdir/RSEM.genes.results \

/home/ymitaka/File02_val_1_RSEM_outdir/RSEM.genes.results \

... /home/ymitaka/File10_val_1_RSEM_outdir/RSEM.genes.results > GeneExpression_Count.csv

8. 以下のコマンドを入力し、全サンプルのRSEM.isoforms.resultsをマージして一つの表（CSVファイル）を作成

/home/usrname/local/bin/RSEM-1.3.3/rsem-generate-data-matrix \

/home/ymitaka/File01_val_1_RSEM_outdir/RSEM.isoforms.results \

/home/ymitaka/File02_val_1_RSEM_outdir/RSEM.isoforms.results \

... /home/ymitaka/File10_val_1_RSEM_outdir/RSEM.isoforms.results > IsoformExpression_Count.csv

9. 作成したCSVファイルやcontigの塩基配列のFASTAファイル、自分のコンピューターにダウンロードし、アノテーションなど諸々の解析に使う。

（メモ）

CSVファイルの中を見てみると、行が各contig（またはisoform）、列名がサンプルになっている。このままでは列名が上記のパスになって長すぎるため、余計な部分を消す。

また、上記の2つのCSVファイルはカンマ区切りになっていないため、エクセルで表示すると各行の1列目のセルに長いタブ区切りの文字列が入ってしまっている。そのため、エクセルの"データ"タブの"データファイル指定（Power Query）"を用いて、カンマ区切りになるよう編集した上で再度CSV形式で保存し直す必要あり。