23.miRNA-seq解析: 既知miRNAの発現量比較

小川

『はじめに』

miRNA(microRNA)は遺伝子の発現を調節する機能を有し、癌や神経疾患などのバイオマーカーとして多くの研究がされています。miRNA-seq(microRNA sequencing)解析は、次世代シーケンサーを利用したRNA解析の一つで、発現量の比較やターゲット遺伝子の予測などが行われています。
今回は弊社が提供するVirtual Machine (VM)上で行ったmiRNAの発現量比較についてご紹介します。miRNAの検出には、DDBJ上で公開されている健常者の正常細胞と癌患者の癌細胞から得られたFASTQデータを利用しました。 VMについては【Windows上で動作するNGSデータ解析環境】(https://www.dynacom.co.jp/product_service/development/win_ngs.html)を参照ください。
*ファイルパスやユーザー名などはVM環境でのものになりますので、環境が異なる場合は適宜修正ください。

1.データのセットアップ

DDBJに公開されているmiRNA-seqのFASTQファイルをダウンロードします。
(https://trace.ddbj.nig.ac.jp/DRASearch/submission?acc=DRA001043)
$ wget ftp://ftp.ddbj.nig.ac.jp/ddbj_database/dra/fastq/DRA001/DRA001043/DRX011925/\
DRR013035.fastq.bz2
*DRR013035.fastq.bz2のデータを例として扱っています。本解析を行うのにはDRR013035.fastq.bz2を含めた10サンプル分のFASTQファイルのダウンロードおよびクオリティーコントロール、マッピング作業が必要になります。

健常者の正常細胞のmiRNA-seq (DRR013035-DRR013039)データと癌患者の癌細胞のmiRNA-seq (DRR013040-DRR013044)データの計10サンプル利用します。

次にmiRNAのリファレンス配列をmiRBaseからダウンロードします。
(http://www.mirbase.org/ftp.shtml)

$ wget ftp://mirbase.org/pub/mirbase/CURRENT/mature.fa.gz
$ zcat mature.fa.gz | grep -A 1 '^>hs' | grep -v '^--' > mature_hsa.fa
$ fasta_nucleotide_changer -i mature_hsa.fa -o mature_hsa_changeU.fa -d
$ bowtie-build -f mature_hsa_changeU.fa mature_hsa_changeU

今回は完成型(mature)miRNAからヒトの配列を抽出し、それをリファレンス配列にします。またマッピングできるようにリファレンス配列の「U」を「T」に変換した後に、Bowtie用のインデックスを作成します。塩基の変換はFASTX-Toolkitを利用しました。

2.クオリティーコントロール

マッピングを行う前にクオリティーコントロールを行います。
$ fastqc --nogroup -o ./ DRR013035.fastq.bz2
$ cutadapt \
-b TACAGTCCGACGATCTCGTATGCCGTCTTC \
-b CTACAGTCCGACGATCTCGTATGCCGTCTT \
-b TCGTATGCCGTCTTCTGCTTGAAAAAAAAA \
-m 16 DRR013035.fastq.bz2 > DRR013035_noadapt.fastq
$ java -Xmx4g -jar /opt/bio/local/Trimmomatic-0.33/trimmomatic-0.33.jar SE -phred33 \
DRR013035_noadapt.fastq DRR013035_noadapt_trim.fastq \
CROP:28 LEADING:20 TRAILING:20 AVGQUAL:20 SLIDINGWINDOW:4:15 MINLEN:16
$ fastqc --nogroup -o ./ DRR013035_noadapt_trim.fastq

各サンプルの品質をFastQCを用いてチェックした後に、cutadaptおよびTrimmomaticを用いてクオリティーコントロールを行います。今回は、Illumina社の次世代シーケンサーでよく使用されている3つのアダプター配列を除去し、塩基長およびquality scoreで配列のフィルタリングを行いました。

3.マッピング

Bowtieを使用してマッピングを行います。
$ bowtie -S -m 1 -v 2 -a --best --strata mature_hsa_changeU DRR013035_noadapt_trim.fastq | samtools view -Sb -o DRR013035_noadapt_trim_bowtie.bam -
$ samtools view -b -F 4 DRR013035_noadapt_trim_bowtie.bam > DRR013035_noadapt_trim_bowtie_mapped.bam
$ samtools sort DRR013035_noadapt_trim_bowtie_mapped.bam DRR013035_noadapt_trim_bowtie_mapped.sort
$ samtools index DRR013035_noadapt_trim_bowtie_mapped.sort.bam

マッピング結果はSAM形式で出力されます。出力される結果をパイプで渡し、Samtoolsを利用して、SAM形式のファイルをBAM形式に変換します。また、保存されたBAM形式のファイルからマッピングされたmiRNAのみを抽出します。さらにソートおよびインデックス処理を行います。ソートおよびインデックス処理を行うことにより、ファイルサイズを圧縮し、効率的にデータを利用することができます。

4.発現量比較

発現量の比較を行います。
$ samtools idxstats DRR013035_noadapt_trim_bowtie_mapped.sort.bam > DRR013035_bowtie_expression.txt

Samtoolsのidxstatsオプションを使用し、各miRNAのリード数をカウントします。これらの作業をサンプル分行います。今回は10サンプル分の*bowtie_expression.txtファイルを作成しました。
次に作成されたファイルから各サンプルごとのmiRNAの発現量を一つにまとめたタブ区切りの一覧を作成します。今回の解析で、2005個のmiRNAの発現量の一覧が作成されました。一覧の作成にはperlを利用しました。

次にサンプルを正常細胞の群と癌細胞の群に分けて、2群間の発現量の比較を行います。
解析には統計解析ソフトRedgeRパッケージ、また散布図の作成にはggplot2パッケージを利用しました。


図は正常細胞の発現量と比較した癌細胞の発現量になります。プロットされている点は、正常細胞のmiRNA発現量/癌細胞のmiRNA発現量の比率を表しています。図中の赤点は癌細胞で有意(logFC > 2.0 かつ FDR < 0.01)に発現量が上がっているものを示し、青色の点は有意(logFC < -2.0 かつ FDR < 0.01)に発現量が下がっているものを示しています。この結果、正常細胞と比較して発現量が有意に上昇している160個のmiRNAと有意に減少している84個のmiRNAが検出されました。

今回はmiRNAの発現量比較まで行いました。miRNA-seq解析は発現量比較の他に新規miRNA候補の予測やターゲット遺伝子の予測などがあります。 弊社では、NGS解析の出張レクチャーサービスやVMを利用したLinux/NGS解析の基本セミナーなども行っております。ご興味のある方はお問い合わせください。

最後に、当コラムを作成するにあたって使用したmiRNA-seqのサンプルデータを公開してくださいました鹿児島大学および千葉大学の研究者の皆様方に深謝いたします。

参考文献:
T. Itesako1, N. Seki, H. Yoshino, T. Chiyomaru, T. Yamasaki, H. Hidaka, T. Yonezawa, N. Nohata, T. Kinoshita, M. Nakagawa, H. Enokida: The MicroRNA Expression Signature of Bladder Cancer by Deep Sequencing: The Functional Significance of the miR-195/497 Cluster. PLOS ONE. February 2014 | Volume 9 | Issue 2 | e84311