pythonってすごいね

RNAseqを用いた遺伝子発現量解析、機械学習を用いた回帰、分類などの解析を中心に記事を書いていきたいです!

Trimomaticを用いたアダプター配列の除去 (RNAseq)

基本文

シングルエンド (SE)

java -jar <trimmomatic.jarのパス> SE -threads <スレッド数> -phred33 or -phred64 -trimlog <logの保存先> <inputのパス> <outputのパス> <option> 

現在は、だいたい-phred33だと思われる

e.g. 
/
java -jar trimmomatic-0.39.jar SE -threads 4 -phred33 -trimlog log.txt input.fa output.fa
/

ペアエンド (PE)

java -jar <trimmomatic.jarのパス> PE -threads <スレッド数> -phred33 or -phred64 -trimlog <logの保存先> <input 1のパス> <input 2のパス> <ペアのリード output 1のパス> <非ペアのリード output 1のパス> <ペアのリード output 2のパス> <非ペアの output 2のパス> <option>

e.g.
/
java -jar trimmomatic-0.39.jar PE -threads 4 -phred33 -trimlog log.txt  input_1.fa input_2.fa paired_output_1.fa unpaired_output_1.fa paired_output_2.fa unpaired_output_2.fa <option>
/


詳細な設定

アダプター配列の除去

ILLUMINACLIP:<fast形式のファイル.fa>:<seed mismatch>:<palindrome clip threshold>:<simple clip threshold>

e.g.
ILLUMINACLIP:adapter_sequence.fa:2:30:10

ファスタ形式のファイルの例

>adapter_sequence
AAAAAAAAAAAAGGGGGGGGGGGGG

Forward、Reverseかを指定したい場合、

/1を末尾につけることで、forwardのアダプター配列から取り除かれる
/2を末尾につけることで、reverceのアダプター配列から取り除かれる

>forward/1
CCCCCCCCAAAAATT

>forward/2
TTTTTTTGGGGGAAAA

>reverse/1
AAAGAATTGGGAAAA

>reverse/2
AGGGAATTGGGAAAA

ペアエンドで配列を指定する際は、Prefixを共通にする。(/1、/2以外の名前を同一にしておく)

seed mismatch

ミスマッチをいくつまで許容するか?

palindrome clip threshold

何塩基(?)ペアエンドのリード間でアダプター配列が一致しているか

simple clip threshold

何塩基アダプター配列が一致しているか?

リードのクオリティに応じたフィルタリング

SLIDINGWINDOW:<検索ウィンドウのサイズ>:<クオリティ>

指定したウィンドウサイズで検索を行い、指定したクオリティ以下の領域が存在した場合、それ以降を解析から除外する。

e.g.
SLIDINGWINDOW:4:15

LEADING:<リードのクオリティ>

5'末端から開始し、指定したクオリティを下回った場合、対象塩基を除去、次の塩基のクオリティを算出指定く。

TRAILING:<リードのクオリティ>

3'末端から開始し、指定したクオリティを下回った場合、対象塩基を除去、次の塩基のクオリティを算出指定く。

e.g.
LEADING:3
TRAILING:3

塩基長の指定

CROP:<塩基長>

リードの頭から、指定した塩基長を取得する。

HEADCROP:<塩基長>

リードの頭から指定の塩基を切り落とす

MINLEN:<塩基長>

最小長の設定

e.g.
CROP:50
HEADCROP:10
MINLEN:36

参考

Trimmomatic Manual: V0.32

USADELLAB.org