マクロジェンにRNA-seqを委託した話とアダプター配列除去方法

今回はマクロジェンに委託したRNA-seqについての話。弊ラボではそこにRNA-seqを委託しているので、ラボメンに向けての備忘録も兼ねて。

シーケンスの状態

弊研究室ではRNAシーケンスをマクロジェンに委託している。シーケンサはNovaSeq6000で100bpもしくは150bpで読んでいる。相乗りで委託することで返ってくるまでの時間はかかるが、比較的安価に委託できる。質については他と比較したことがないのでわからないが、よく読めていると思うし、クォリティチェックもしっかりやってくれるので安心である。

ただ、100bpで委託したときはQC(クォリティチェック)のときに特別問題はなかったが、150bpで委託したときにはアダプター配列が含まれていると出た。仕様書を見るとどちらも「アダプター配列は除去していない」と書かれているので、本来は含まれるものなのであろう。100bpと150bpで含まれる・含まれない理由をマクロジェン側に尋ねると、

今回の条件では、アダプターの両側から、シーケンスプライマー直下の
150ntのシーケンスを行っております。
そのため、インサートサイズが150bpに満たない場合は、
反対側のアダプター配列がシーケンスされることとなります。
インサートサイズが150bpを越える場合は、反対側のアダプター配列に
到達しないため、アダプター配列がリード情報に含まれなくなります。

ということらしい。わかりやすく図説するとこのようになる。と思う。

この読みたい配列(DNA Insert)の長さはまちまちなので、長くシーケンスするとアダプター配列まで読まれてしまうものが出てくるということらしい。

なるほど。

アダプター配列を取り除く

で、「このアダプター配列の何が問題なの？」ってことだが、「マッピング率がすこぶる悪くなる」ってことだ。今回100bpやっていて、追加で150bpの長さで読んでもらったので気づけたが、明らかにマッピング率が悪くなる。

100bp長のときは90%あったマッピング率が150bp長でシーケンスしたときのをマッピングすると65%にまで落ち込んだ。はて。と思って今回のこれにつながったわけだが、何はともあれアダプター配列を取り除いてやらねばならない。

Quality Control のためのソフトウエアである“prinseq++”を用いて無理矢理その部分を一括60bp削ってやるとかいう荒業でとりあえずやってみたところ、マッピング率は90%弱まで回復した。

まあでもそんなことしてしまったら大事な配列の部分まで失ってしまうのでアダプター配列だけを取り除きたい。そこでTrimmomaticというソフトウエアを用いてアダプター配列だけを除去してやる。

Trimmomatic の使い方

公式ページからダウンロード。適当なディレクトリに入れておく。

他のサイトが解説していることから、今回はペアエンドの使い方というか、使用例のみ触れておく。

$ java -jar ~/Bioinformatics/tools/Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 24 -phred33 -trimlog ./log.txt ./INPUT_1.fastq.gz ./INPUT_2.fastq.gz paired_output_1.fq unpaired_output_1.fq paired_output_2.fq unpaired_output_2.fq ILLUMINACLIP:../adapter.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:36

で、この肝心の’adapter.fa’だが、使うライブラリに合わせて記述する。今回はIllumina社のTruSeq Stranded mRNA ということだったので、このマニュアルに記載されている配列をMulti fasta で記述した。

>adapter/1
AGATCGGAAGAGCACACGTCTGAACTCCAGTCA

>adapter/2
AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT

で、実行するとこんな感じの標準出力。

TrimmomaticPE: Started with arguments:
 -threads 24 -phred33 -trimlog ./log.txt ./INPUT_1.fastq.gz ./INPUT_2.fastq.gz paired_output_1.fq unpaired_output_1.fq paired_output_2.fq unpaired_output_2.fq ILLUMINACLIP:../adapter.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:15 MINLEN:36
Using Long Clipping Sequence: 'AGATCGGAAGAGCACACGTCTGAACTCCAGTCA'
Using Long Clipping Sequence: 'AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGT'
ILLUMINACLIP: Using 0 prefix pairs, 0 forward/reverse sequences, 1 forward only sequences, 1 reverse only sequences
Input Read Pairs: 24065880 Both Surviving: 23792505 (98.86%) Forward Only Surviving: 173500 (0.72%) Reverse Only Surviving: 82669 (0.34%) Dropped: 17206 (0.07%)
TrimmomaticPE: Completed successfully

これ、ログを出させているけど4GB近くになったので正直いるかどうか微妙。

その結果をFastQCで見てみると