高知大学　R5年度　DS/DXセミナー　

統計基礎、bulk RNA-Seq解析ハンズオントレーニング

日時

2024年3月5日(火) 10:30−18:00\ 2024年3月6日(水) 10:30−18:00

場所

高知大学　物部キャンパス　農場棟１階講義室

講義の概要と学習目標

公共データベースに格納されているbulk RNA-Seqデータを取得し、発現変動遺伝子（Differentially expressed genes; DEGs）の検出、DEGsを用いた各種エンリッチメント解析を行う。

用意するもの

インターネット接続可能なコンピュータ１台（OS: Windows 11, or macOS）

スケジュール

１日目

10:30-12:00　統計学基礎\ 13:00-18:00　bulk RNA-Seq解析ハンズオントレーニング事前説明

２日目

10:30-12:00　bulk RNA-Seq解析ハンズオントレーニング\ 13:00-18:00　bulk RNA-Seq解析ハンズオントレーニング

統計基礎講義資料

下記よりダウンロード下さい。\ https://github.com/yosui-nojima/bulk_RNA-Seq_Kochi_Univ._2024-0305-0306/blob/main/Kochi_DSDX_Seminar_20240305_compressed.pdf

bulk RNA-Seq解析ハンズオントレーニング

1. 環境構築

Windows 11

Windows 11でUnix環境を構築するため、Cygwinを使用します。

- Cygwinのインストール

https://www.cygwin.com/install.htmlにアクセス
「setup-x86_64.exe」をクリックし、setup-x86_64.exeをダウンロード
ダウンロードディレクトリに保存された「setup-x86_64.exe」をダブルクリック
「はい」をクリック

「次へ」をクリック

下記を選択し、「次へ」をクリック

下記のように選択し、「次へ」をクリック

下記のように選択し、「次へ」をクリック。ただし、画面上のxxxは任意のユーザー名に変更すること。

下記のように選択し、「次へ」をクリック

日本のミラーサーバーを選択し、「次へ」をクリック

Cygwinはデフォルトでは、必要最低限のパッケージしかインストールしてくれません。解析に必要なパッケージを個別に指定してインストールします。指定するパッケージは下記の８つです。該当するパッケージのDefaultをプルダウンでInstallに変更します。
- Admin
- Archive
- Base
- Devel
- Editors
- Graphics
- Libs
- Perl

下記のように選択し、「次へ」をクリック

「次へ」をクリック

ダウンロードが開始されます。

下記のようにダウンロードが完了しなかったパッケージが表示されることがありますので、ダウンロードが進んでいるか定期的にチェックして下さい。表示された場合は、「Retry」を選択し、再度インストールを開始します。

再びダウンロードが開始されます。合計で２、３時間かかります。

完了をクリック。

- 作業ディレクトリの作成

Cygwinを起動します。
ホームディレクトリ下に作業ディレクトリを作成するため、下記のコマンドを実行
```
cd
mkdir bulksem
```

- NCBI SRA Toolkitのインストール

下記のコマンドを実行

cd ~/bulksem
curl -OL https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/3.0.10/sratoolkit.3.0.10-win64.zip
unzip ./sratoolkit.3.0.10-win64.zip

今回使用するsratoolkitのプログラムは、prefetchとfastq-dumpの２つです。\ 解凍後、下記コマンドを実行します。
```
./sratoolkit.3.0.10-win64/bin/prefetch -V
```
```
./sratoolkit.3.0.10-win64/bin/fastq-dump -V
```
以下のメッセージがそれぞれ表示されれば、NCBI SRA Toolkitのインストールは完了です。

- FastQCのインストール

下記コマンドを実行

cd ~/bulksem
curl -OL https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.12.1.zip
unzip fastqc_v0.12.1.zip
cd ./FastQC
./fastqc -v

FastQC v0.12.1のメッセージが表示されれば、FastQCのインストールのインストール完了

#2024年3月6日更新 Windowsユーザーは下記サイトのFastQC v0.11.9 (Win/Linux zip file)をクリックして、ファイルをダウンロードして下さい。\ https://kb.10xgenomics.com/hc/en-us/articles/360048465032-When-to-run-FastQC\ ダウンロードしたファイルはbulksemディレクトリには移動させず、ダウンロードディレクト下に置いたまま解凍し、run_fastqc.batをダブルクリックして実行して下さい。アプリケーションが起動されます。

- javaのインストール

Windowsにjavaがインストールされていない場合は、https://www.java.com/ja/download/にアクセスし、jre-8u401-windows-x64.exeをダウンロードする（2024年2月28日時点）。
ダウンロードディレクトリに保存された「jre-8u401-windows-x64.exe」をダブルクリック
「はい」をクリック

「インストール」をクリック

インストール中

「閉じる」をクリック

- Trimmomaticのインストール

下記コマンドを実行

cd ~/bulksem
curl -OL http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip
java -jar ./Trimmomatic-0.39/trimmomatic-0.39.jar -version

以下のメッセージが表示されれば、Trimmomaticのインストール完了

- HISAT2のインストールとコンパイル

https://cloud.biohpc.swmed.edu/index.php/s/fE9QCsX3NH4QwBi/downloadにアクセスして、hisat2-2.2.1-source.zipをダウンロードします。
ダウンロードディレクトリに保存されたhisat2-2.2.1-source.zipを、作業ディレクトリであるC:\cygwin64\home\xxx\bulksemにドラッグ・アンド・ドロップでして下さい（xxxは各端末で設定されている任意のユーザー名です。）
下記コマンドを実行
```
cd　~/bulksem
unzip ./hisat2-2.2.1-source.zip
cd ./hisat2-2.2.1
make
```
makeでのコンパイルには30分ほどかかります。
コンパイル後、下記コマンドを実行します。
```
./hisat2 -h
```
以下のメッセージが表示されれば、hisat2のインストール・コンパイルは完了です。\ （最初の数行のみ表示させています。）

コンパイル後、下記コマンドを実行します。
```
./hisat2-build -h
```
以下のメッセージが表示されれば、hisat2-buildのインストール・コンパイルは完了です。\ （最初の数行のみ表示させています。）

- featureCountsのインストール

https://sourceforge.net/projects/subread/files/subread-2.0.6/subread-2.0.6-Windows-x86_64.tar.gz/downloadにアクセスし、subread-2.0.6-Windows-x86_64.tar.gzをダウンロードします。
ダウンロードディレクトリに保存されたsubread-2.0.6-Windows-x86_64.tar.gzを、作業ディレクトリであるC:\cygwin64\home\xxx\bulksemにドラッグ・アンド・ドロップでして下さい（xxxは各端末で設定されている任意のユーザー名です。）
下記コマンドを実行
```
cd ~/bulksem
tar -zxvf ./subread-2.0.6-Windows-x86_64.tar.gz
```
解凍後、下記コマンドを実行
```
./subread-2.0.6-Windows-x86_64/bin/featureCounts -v
```
以下のメッセージが表示されれば、featureCountsのインストール・コンパイルは完了です。

- PATHの追加

絶対PATHの入力を不要にするため、コマンドごとにPATHを追加しておきます。\ 下記のコマンドを実行します。（ただし、xxxは各端末で設定されている任意のユーザー名に変更して下さい。）

export PATH=/home/xxx/bulksem/sratoolkit.3.0.10-win64/bin:/home/xxx/bulksem/FastQC:/home/xxx/bulksem/hisat2-2.2.1:/home/xxx/bulksem/subread-2.0.6-Windows-x86_64/bin:$PATH

ただし、Cygwinを終了するとリセットされるため、再度実行する必要があります。\ （Cygwin終了の有無に関わらずPATHを追加することも可能ですが、今回は行いません。）

- Rのインストール

データ解析の再現性

『再現性』とは、同一の結果が同一の実験手法・解析手法によって得られるとき、それら結果の一致の度合いの高さを示す。\ 自分の解析結果を研究室内や会社内の他の人が同じ解析をする場合、エクセルなどの表計算ソフトにおけるメニュー操作やコピー＆ペーストを通して行ったデータの集計・加工・分析およびグラフ化（可視化）の作業は記録できず、結果の再現性が低い。\ 一方、RやPythonといったプログラミング言語を用いた解析では、スクリプトを作成することでデータの読み込みから結果の出力まで、「手作業」を極力排除してデータ解析を自動化することができ、結果の再現性が高い。

Rについて

Rはデータ分析や統計解析のために開発されたソフトウェアで、プログラミング言語としても十分な機能を備えている。
プログラミング言語といってもCやJavaなどの言語よりも比較的簡単。順番に処理を記述していけば一通りの分析が可能であり、プログラミング言語としうより「分析ツール」という感覚で使用している人も多い。
無料で入手できる統計解析に特化したプログラミング言語で、統計解析で最も広く使われている。
基本的な統計解析機能が標準パッケージに含まれており、初期状態で一通りの統計分析を行うことが可能。
様々な分野に適した拡張パッケージが提供されており、適宜インストールして使用することが可能。
Rは開発者のRoss Ihaka、Robert Clifford Gentlemanにより開発され、Rという名称は両者のイニシャルでもある。
現在は、R Development Core Teamによってメンテナンス・拡張がされている。

https://www.r-project.org/をクリック。
『download R』をクリック
任意のミラーサイトのリンクをクリック
『Download R for Windows』をクリック。※以降のインストール方法はOSがWindowsの場合です。その他のOSの場合は、適切なリンクをクリック。
『install R for the first time』をクリック。
『Download R-4.3.0 for Windows』をクリック。
ダウンロードフォルダにある『R-4.3.0-win.exe』をダブルクリック。
『はい』をクリック。
『日本語』が選択されていることを確認して、『OK』をクリック
『次へ』をクリック。
『次へ』をクリック。
３つともチェックされていることを確認して、『次へ』をクリック。
『いいえ』が選択されていることを確認して、『次へ』をクリック。
『次へ』をクリック。
『次へ』をクリック。
『完了』をクリック。

- RStudioのインストール

Rは単独で実行することも可能ですが、統合開発環境（Integrated Development Environment; IDE）というプログラムの記述や実行、デバックなどに必要なツールが一体化されているツールを利用して使用する方が便利です。\ IDEとして、最も一般的なのが『RStudio』です。

https://posit.co/download/rstudio-desktop/をクリック。
『DOWNLOAD RSTUDIO DESKTOP FOR WINDOWS』をクリック。※以降のインストール方法はOSがWindowsの場合です。その他のOSの場合は、適切なリンクをクリック。
ダウンロードフォルダにある『RStudio-2023.03.0-386.exe』をダブルクリック。
『はい』をクリック。
『次へ』をクリック。
『次へ』をクリック。
『インストール』をクリック。
『完了』をクリック。

RStudioによるRの起動

デスクトップの左下にある検索窓に「rstudio」と入力。
『タスクバーにピン留めする』をクリック。
『開く』をクリック。
『Use your machine's default 64-bit version of R』が選択されていることを確認して、『OK』をクリック。
クラッシュレポートを提供するかどうかを聞いている。どちらでもかまわない。
以下の画面が表示されればインストール成功。

- Rパッケージのインストール

『Console』に直接入力して実行してもかまわないが、スクリプトとして記録・保存できないため、\ タブから、File → New File → R Scriptを選択し、クリック。\ 下記コードを入力し、「Run」をクリック。

install.packages("BiocManager")
install.packages("ggplot2")
install.packages("reshape2")
BiocManager::install("GenomicFeatures")
BiocManager::install("clusterProfiler")
BiocManager::install("biomaRt")
BiocManager::install("org.Hs.eg.db")
BiocManager::install("DOSE")

macOS

macOSはUnixベースのOSですので、Cygwinのインストールは必要ありません。\ また、javaはデフォルトでインストールされているため、基本的にはインストール不要です。

ターミナルの起動

まずはターミナルを立ち上げます。\ Finder > アプリケーション > ユーティリティ > ターミナル\ 下記のアイコンのアプリケーションをダブルクリックして立ち上げます。\ スクリーンショット 2022-03-18 14 59 09 \ 以降は、ターミナルでコマンドラインインターフェース（Command Line Interface; CLI）により各種解析を実行します。

- 作業ディレクトリの作成

ホームディレクトリ下に作業ディレクトリを作成するため、下記のコマンドを実行

cd
mkdir bulksem

- NCBI SRA Toolkitのインストール

https://ftp-trace.ncbi.nlm.nih.gov/sra/sdk/3.0.10/sratoolkit.3.0.10-mac-x86_64.tar.gzにアクセスし、sratoolkit.3.0.10-mac-x86_64.tar.gzをダウンロードします。
Downloadsディレクトリにダウンロードされているため、下記コマンドを実行して作業ディレクトリに移動し、解凍します。
```
mv ~/Downloads/sratoolkit.3.0.10-mac-x86_64.tar.gz ~/bulksem
cd ~/bulksem
tar -zxvf ./sratoolkit.3.0.10-mac-x86_64.tar.gz
```

下記コマンドを実行

cd ~/bulksem
cd ./sratoolkit.3.0.10-mac-x86_64/bin/
./prefetch -V

実行後、以下のポップアップが表示された場合は、「キャンセル」をクリック

Finder > アプリケーション > システム環境設定 > セキュリティとプライバシーをクリックし、以下の画面が表示されているので、「このまま許可」をクリック

再度下記コマンドを実行
```
./prefetch -V
```
以下の画面が表示されているので、「開く」をクリック

以下の画面が表示されるので、「キャンセル」をクリック

「システム環境設定」→「セキュリティとプライバシー」をクリックし、以下の画面が表示されているので、「このまま許可」をクリック

再度下記コマンドを実行
```
./prefetch -V
```
以下の画面が表示されているので、「開く」をクリック

以下のメッセージが表示されれば、prefetchコマンドが実行可能になります。

- FastQCのインストール

下記コマンドを実行

cd ~/bulksem
curl -OL https://www.bioinformatics.babraham.ac.uk/projects/fastqc/fastqc_v0.12.1.zip
unzip fastqc_v0.12.1.zip
cd ./FastQC
./fastqc -v

FastQC v0.12.1のメッセージが表示されれば、FastQCのインストールのインストール完了

- Trimmomaticのインストール

下記コマンドを実行

cd ~/bulksem
curl -OL http://www.usadellab.org/cms/uploads/supplementary/Trimmomatic/Trimmomatic-0.39.zip
unzip Trimmomatic-0.39.zip
java -jar ./Trimmomatic-0.39/trimmomatic-0.39.jar -version

0.39のメッセージが表示されれば、Trimmomaticのインストール完了

- HISAT2のインストール

https://cloud.biohpc.swmed.edu/index.php/s/zMgEtnF6LjnjFrr/downloadにアクセスし、hisat2-2.2.1-OSX_x86_64.zipをダウンロードします。
Downloadsディレクトリにダウンロードされているため、下記コマンドを実行して作業ディレクトリに移動し、解凍します。
```
mv ~/Downloads/hisat2-2.2.1-OSX_x86_64.zip ~/bulksem
cd ~/bulksem
unzip ./hisat2-2.2.1-OSX_x86_64.zip
```

下記コマンドを実行

cd ~/bulksem
cd ./hisat2-2.2.1
./hisat2 --version

NCBI SRA Toolkitのprefetchと同じ要領でhisat2に実行許可を与えます。
再度下記コマンドを実行
```
./hisat2 --version
```
以下のメッセージが表示されれば、hisat2コマンドが実行可能になります。

下記コマンドを実行
```
./hisat2-build --version
```
同じ要領で実行許可を与えます。
再度下記コマンドを実行
```
./hisat2-build --version
```
以下のメッセージが表示されれば、hisat2-buildコマンドが実行可能になります。

- featureCountsのインストール

https://sourceforge.net/projects/subread/files/subread-2.0.6/subread-2.0.6-macOS-x86_64.tar.gz/downloadにアクセスし、subread-2.0.6-macOS-x86_64.tar.gzをダウンロードします。
Downloadsディレクトリにダウンロードされているため、下記コマンドを実行して作業ディレクトリに移動し、解凍します。
```
mv ~/Downloads/subread-2.0.6-macOS-x86_64.tar.gz ~/bulksem
cd ~/bulksem
tar -zxvf  ./subread-2.0.6-macOS-x86_64.tar.gz
```

下記コマンドを実行

cd ~/bulksem
cd ./subread-2.0.6-macOS-x86_64/bin
./featureCounts -v

NCBI SRA Toolkitのprefetchと同じ要領でfeatureCountsに実行許可を与えます。
再度下記コマンドを実行
```
./featureCounts -v
```
featureCounts v2.0.6のメッセージが表示されれば、hisat2コマンドが実行可能になります。

- PATHの追加

export PATH=/Users/xxx/bulksem/sratoolkit.3.0.10-mac-x86_64/bin:/Users/xxx/bulksem/FastQC:/Users/xxx/bulksem/hisat2-2.2.1:/Users/xxx/bulksem/subread-2.0.6-macOS-x86_64/bin:$PATH

ただし、ターミナルを終了するとリセットされるため、再度実行する必要があります。\ （ターミナル終了の有無に関わらずPATHを追加することも可能ですが、今回は行いません。）

2. 使用データ

下記のpaired-endでシーケンスされた２サンプルのデータを使用します。\ 今回のセミナー用に公共データから１サンプル１０万リードランダムサンプリングしたものです。\ 下記のコマンドを実行

cd ~/bulksem
curl -OL https://github.com/nojima-q/2021-12-13-15_PBL_analysis/raw/main/sample1_1_100K.fastq.gz
curl -OL https://github.com/nojima-q/2021-12-13-15_PBL_analysis/raw/main/sample1_2_100K.fastq.gz
curl -OL https://github.com/nojima-q/2021-12-13-15_PBL_analysis/raw/main/sample2_1_100K.fastq.gz
curl -OL https://github.com/nojima-q/2021-12-13-15_PBL_analysis/raw/main/sample2_2_100K.fastq.gz

3. 公共データベースの紹介

20190605_metacore

NCBI SRA (https://www.ncbi.nlm.nih.gov/sra)

DDBJ Sequence Read Archive (https://www.ddbj.nig.ac.jp/dra/index.html)

EBI ENA (https://www.ebi.ac.uk/ena/browser/home)

※データ検索は、スクリーン上で実演します。

データ取得方法

１．解析したデータセットのAccession numberをNCBI SRA Run Selectorに入力し、『Search』をクリック。

２．データセット内の全てのデータを解析する場合は、Total行の『Accession List』をクリックし、サンプルごとのAccession numberが記載されたテキストファイル（SRR_Acc_List.txt）をダウンロードする。一部のデータのみ解析する場合は、必要なデータに☑をいれSelected行の『Accession List』をクリックし、SRR_Acc_List.txtをダウンロードする。

SRR_Acc_List.txtの内容\

SRA Toolkitのprefetch、fastq-dumpを使ってデータを取得する。まず、prefetchでsraファイルがダウンロードされる。１つのsraファイルが20GBを超える場合は、-Xまたは--max-size 50Gなどのように最大数を変更する。--option-fileは使わずに直接Accession numberを入力して個別にダウンロードすることも可能です。

prefetch  --option-file ~/Downloads/SRR_Acc_list.txt

次にfastq-dumpでsraファイルからfastqファイルを取得する。PATHにfastqファイルを格納したいディレクトリのパスを記載。\

find . -name '*.sra' -exec fastq-dump --gzip --split-files --outdir ./PATH {} \;

--gzip：圧縮ファイルとして出力する
--split-files：レイアウトがpaired-endの際に指定する。single-endのときは不要。

4. FASTQファイルの形式、クオリティーコントロールについて

１．FASTQファイルとは

@SRR8615662.1 1 length=101\ TGATGGCCCTGCCTTCGTGGGAACAGAGGCTAAGGCCTTGAG\ +SRR8615662.1 1 length=101\ CCCFFFFFHHHHHJJJJIJJJJIJIJJJJJJJJJJJJJJJJJFIIGIJJIJIJJIJJJJJHHHHHFFFFF\ \ １行目：＠から始まり、リードIDが記載されている\ ２行目：シーケンスしたリードの塩基配列\ ３行目：＋を記載\ ４行目：２行目に記述した各塩基のクオリティ値

２．FASTQファイルのクオリティーコントロール

FastQCを使って、FASTQファイルの品質を確認します。\ 下記のコマンドを実行します。

cd ~/bulksem
mkdir fastqc_results
fastqc -t 4 -o ./fastqc_results/ ./sample1_1_100K.fastq.gz

-t：スレッド数（使用するPC環境に合わせて設定して下さい。）
-o：出力ディレクトリ

計算が終わるとAnalysis complete for sample1_1_100K.fastq.gzと表示れます。 fastqc_resultsディレクトリに、sample1_1_100K_fastqc.htmlというファイルが作成されています。\ これを適当なWeb Browserで開いて下さい。

3. FastQC解析結果

各QC結果の説明は、下記のブログでよくまとめられています。\ https://imamachi-n.hatenablog.com/entry/2017/03/27/234350

Per base sequence quality

リードの各塩基のクオリティスコアを示しています。 Phred quality scoreがだいたいグリーンの領域（Scoreが28以上）に収まっているかどうか確認します。結果として、クオリティが低いリードは含まれていないことが確認できます。\

Per tile sequence quality

フローセルの各タイルごとのクオリティスコアを示しています。\ Illumina社製の次世代シーケンサーでは、「フローセル」と呼ばれるガラス基板上でDNA合成反応を行います。このフローセルは「タイル」と呼ばれる単位に区切られており、各タイルごとに塩基配列を決定しています。\ シーケンスをかけたときに、例えば特定のタイルに気泡やゴミが入っているとクオリティの低下が見られることがあります。\ 特定のタイルで著しいクオリティの低下が見られる場合は、シークエンス時に上記のような問題があったと考えられます。

Per sequence quality scores

各リードの全長のクオリティスコアの分布を示しています。

Per base sequence content

各塩基のA, T, G, Cの含有量を示しています。 RNA-seqの場合、それぞれの含有量はほぼ25%ずつになりますが、PAR-CLIPのようにRNA結合タンパク質と結合しているRNA配列を抽出してきている場合、それぞれの含有率に偏りが見られます。

Per sequence GC content

リードのGC contentsの分布を示しています。

Per base N content

各塩基中に含まれるNの含有率（塩基を読めなかった箇所）を示しています。

Sequence Length Distribution

リード長の分布を示しています。

Sequence Duplication Levels

Duplidate readsの含まれている数を示しています。

Overrepresented sequences

頻出する特徴配列が示されています。リード中にアダプター配列などが混入している場合、その配列が示されます。

Adapter Content

各塩基ごとに見たときのリード中に含まれているアダプターの割合を示しています。あくまで、FastQCに登録されているアダプター配列しか確認していないので、登録されていないアダプター配列を使っていた場合、そのアダプター配列がリード中に混入していても確認できないことがあります。

4. アダプターの除去およびリードのトリミング

Trimmomaticを使って、アダプターの除去および低スコアな塩基のトリミングを同時に行います。アダプター配列を記載したFATSAファイルはここからダウンロード可能です。（今回は、illumina社のTruSeqシリーズの配列を記載しています。自前データで実行する際は、ライブラリー作製キットで使用している配列が記載されたFASTAファイルを用意して実行して下さい。）\ 下記はpaired-endでシーケンスしたFASTQファイルの場合です。

cd ~/bulksem
curl -OL https://github.com/nojima-q/2021-12-13-15_PBL_analysis/raw/main/Truseq_stranded_totalRNA_adapter.fa
java -jar ./Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 4 -phred33 ./sample1_1_100K.fastq.gz ./sample1_2_100K.fastq.gz ./sample1_1_100K_trim_paired.fastq.gz ./sample1_1_100K_trim_unpaired.fastq.gz ./sample1_2_100K_trim_paired.fastq.gz ./sample1_2_100K_trim_unpaired.fastq.gz ILLUMINACLIP:Truseq_stranded_totalRNA_adapter.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:25

PE：レイアウトがpaired-endシーケンスのときに指定。single-endの場合はSEを指定します。
-threads：スレッド数（使用するPC環境に合わせて設定して下さい。）
ILLUMINACLIP：Truseq_stranded_totalRNA_adapter.faはillumina社のTruSeqシリーズのアダプター配列をFASTA形式で記載してもの。後ろの数字は、許容ミスマッチ数:palindrome clip threshold:simple clip thresholdを表す。
LEADING：5'末端から設定したスコア値未満の塩基をトリム
TRAILING：3'末端から設定したスコア値未満の塩基をトリム
SLIDINGWINDOW：左の数字はウィンドウサイズ、右の数字は平均クオリティ値を表します。ウィンドウサイズの範囲内の塩基のスコア平均が設定値よりも低ければ、3'末端側の全ての塩基をトリム。\
MINLEN：設定値未満の塩基数になったリードを除去する。

awkコマンドを使ったバッチ処理

サンプル数が多い場合一つ一つコマンドを打って処理するのは面倒です。その場合、下記のようにawkコマンドを使って複数サンプルのコマンドを１つのシェルファイルに記載してバッチ処理することができます。\ 変数部分が%sに該当します。%sの数だけ後部にカンマ区切りの$1を記載します。

ls sample*_100K.fastq.gz | cut -f1 -d_ | sort | uniq | awk '{printf ("java -jar ./Trimmomatic-0.39/trimmomatic-0.39.jar PE -threads 4 -phred33 %s_1_100K.fastq.gz %s_2_100K.fastq.gz %s_1_100K_trim_paired.fastq.gz %s_1_100K_trim_unpaired.fastq.gz %s_2_100K_trim_paired.fastq.gz %s_2_100K_trim_unpaired.fastq.gz ILLUMINACLIP:Truseq_stranded_totalRNA_adapter.fa:2:30:10 LEADING:20 TRAILING:20 SLIDINGWINDOW:4:20 MINLEN:25 2> %s_trim.txt \n", $1, $1, $1, $1, $1, $1, $1)}' > trim.sh
sh trim.sh

リードトリミング後のFastQC結果

Per base sequence quality

Per tile sequence quality

Per sequence quality scores

各リードの全長のクオリティスコアの分布を示しています。

Per base sequence content

Per sequence GC content

リードのGC contentsの分布を示しています。

Per base N content

各塩基中に含まれるNの含有率（塩基を読めなかった箇所）を示しています。

Sequence Length Distribution

リード長の分布を示しています。

Sequence Duplication Levels

Duplidate readsの含まれている数を示しています。

Overrepresented sequences

頻出する特徴配列が示されています。リード中にアダプター配列などが混入している場合、その配列が示されます。

Adapter Content

各塩基ごとに見たときのリード中に含まれているアダプターの割合を示しています。あくまで、FastQCに登録されているアダプター配列しか確認していないので、登録されていないアダプター配列を使っていた場合、そのアダプター配列がリード中に混入していても確認できないことがあります。\

4 リファレンスゲノムファイル、アノテーションファイルの取得

ゲノムデータのダウンロードはかなり時間を要するため、本セミナーでは実行しません。下記コマンドは今後の参考にして下さい。

トリミングしたリードは全ゲノム配列が記載されたリファレンスゲノムファイルにマッピングします。\ Ensemblは、様々な生物種のゲノム配列情報が格納されているデータベースです。\ ここからヒトのリファレンスゲノムファイルとアノテーションファイル（遺伝子ごとにコードされている染色体と染色体内でのポジションが記載されたファイル）をダウンロードします。\

まずはEnsemblトップページの『Human』からをクリックします・
次のページの右側『Gene annotation』内の『Download FASTA』をクリック。
ftpサイトの『dna』をクリックすると複数のファイルが閲覧できます。このうち、『Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz』を使用します。
次はアノテーションファイルをダウンロードします。
Humanトップページの『Gene annotation』内の『Download GTF』をクリック。
ftpサイトの『Homo_sapiens.GRCh38.111.gtf.gz』をダウンロードします。
リファレンスゲノムファイルとアノテーションファイルを解凍します。

または、下記のコマンドを実行することでもダウンロード・解凍が可能です。

cd ~/bulksem
curl -OL https://ftp.ensembl.org/pub/release-101/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
curl -OL https://ftp.ensembl.org/pub/release-101/gtf/homo_sapiens/Homo_sapiens.GRCh38.101.gtf.gz
gunzip Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz
gunzip Homo_sapiens.GRCh38.101.gtf.gz

2024年3月1日現在の最新版はrelease-111になります。

5 リファレンスゲノムへのマッピング

5-1 リファレンスゲノムファイルのインデックス化

ゲノムのインデックス化はかなり時間がかかるため、本セミナーでは実行しません。下記コマンドは今後の参考にして下さい。

4でダウンロードしたリファレンスゲノムファイルをインデックス化します。一般的にゲノムデータのサイズは大きくなりがちでありそのままでは処理時間がかかってしまうため、インデックス（目次）ファイルを作成して高速にアクセスできるようにします。\ インデックス化に必要なプログラムは、多くの場合マッピングツールに含まれています。\ HISAT2に含まれるhisat2-buildを使用します。\ インデックス化する前段階として、まずは下記を実行します。それぞれスプライシングサイト、エキソンサイトを抽出するpythonスクリプトです。

cd ~/bulksem
./hisat2-2.2.1/extract_splice_sites.py ./Homo_sapiens.GRCh38.101.gtf > ./Homo_sapiens.GRCh38.101.ss
./hisat2-2.2.1/extract_exons.py ./Homo_sapiens.GRCh38.101.gtf > ./Homo_sapiens.GRCh38.101.exon

上記の２ファイルも使って、インデックス化します。

hisat2-build -p 4 --ss ./Homo_sapiens.GRCh38.101.ss --exon ./Homo_sapiens.GRCh38.101.exon ./Homo_sapiens.GRCh38.dna.primary_assembly.fa ./GRCh38.101

-p：スレッド数（使用するPC環境に合わせて設定して下さい。）
-ss：extract_splice_sites.pyで作成したファイルを指定
-exon：extract_exons.pyで作成したファイルを指定

本セミナーでは下記リンクからインデックス化済みリフェレンスゲノムファイルを取得して下さい。

https://www.dropbox.com/scl/fi/tjuf4u9s81rrg0wtya8s6/GRCh38.101.tar.gz?rlkey=63uz24zbvc7jjssez3tfdmpjl&dl=0\ ダウンロード後、GRCh38.101.tar.gzをWindows11の場合はC:\cygwin64\home\xxx\bulksemに、macOSの場合は~/bulksemに移動して下さい(xxxは任意のユーザー名に変更すること。)\ 下記コマンドを実行

tar -zxvf GRCh38.101.tar.gz

5-2 マッピング

ゲノムへのマッピングを行います。今回は時間短縮のため、１０万リードランダムサンプリングしたFASTQファイル（最初にダウンロードしたファイル）を5-1で作成したインデックス化したリファレンスゲノムにマッピングします。

hisat2 -p 4 --dta -x ./GRCh38.101/GRCh38.101 -1 ./sample1_1_100K_trim_paired.fastq.gz -2 ./sample1_2_100K_trim_paired.fastq.gz -S sample1_hisat2.sam 2> sample1_hisat2_log.txt

-p：スレッド数（使用するPC環境に合わせて設定して下さい。）
--dta：アセンブラーなど下流の解析ツールを使用する際のオプション。またメモリ使用量を改善する。
-x：インデックス化したリファレンスゲノムファイル
-1：FowardリードのFastqファイルを指定
-2：ReverseリードのFastqファイルを指定
-S：出力するSAMファイル名
2>：標準エラー出力。Mapping rate等をテキストファイルとして出力する。

※2024年3月6日追記 ※hisat2でバッチ処理をする場合は下記awkコマンドを参考にして下さい。

 ls sample*_1_100K.fastq.gz | cut -f1 -d_ | awk '{printf ("hisat2 -x ./GRCh38.101/GRCh38.101 -1 ./%s_1_100K.fastq.gz -2 ./%s_2_100K.fastq.gz -S ./%s_hisat2.sam\n", $1, $1, $1)}' > hisat2.sh

6 マッピングデータから遺伝子ごとにリードのカウントデータを取得する

マッピング結果であるSAMファイルからgeneごとまたはtranscriptごとにリードのカウント数を出力します。\ featureCountsを使用します。

curl -OL https://ftp.ensembl.org/pub/release-101/gtf/homo_sapiens/Homo_sapiens.GRCh38.101.gtf.gz
gunzip Homo_sapiens.GRCh38.101.gtf.gz
featureCounts -O -M -T 18 -p -t exon -g gene_id -a ./Homo_sapiens.GRCh38.101.gtf -o sample_count.txt sample*_hisat2.sam

-O：複数のfeatureで定義されている特定のゲノム領域にマッピングされたリードもカウントする
-M：マルチマッピングされたリードもカウントする
-T：スレッド数（使用するPC環境に合わせて設定して下さい。）
-p：paired-endの場合に指定（1リードペアを1カウントする。）
-t：GTFファイルのfeature領域（3カラム目、exon, gene, CDS等）を指定
-g：gene levelでカウントしたい場合はgene_idを、transcript levelでカウントしたい場合はtranscript_idを指定する
-a：アノテーションファイル（GTFファイル等）を指定
-o：出力ファイル名を指定

遺伝子IDとカウント値のみを抽出します。

cut -f1,7- ./sample_count.txt | grep -v ^\# > ./featureCounts_output.txt

7 カウントデータをTPM値に変換する

６で出力したファイルは本セミナー用のスモールデータです。ある疾患の公共RNA-Seqデータ（３９サンプル）のカウント値を出力したファイルfeatureCounts_all_output.txtを用意しました。下記コマンドをCygwin(Windows 11)またはターミナル(macOS)で実行してファイルをダウンロードします。

cd ~/bulksem
curl -OL https://github.com/nojima-q/2021-12-13-15_PBL_analysis/raw/main/featureCounts_all_output.txt

ここからはRStudioで作業します。\ RStudioを起動して下さい。\ まず初めに作業ディレクトリを移動します。\ 下記コードを実行

Windows 11の場合

setwd("C:/cygwin64/home/xxx/bulksem/") #xxxは任意のユーザー名に変更すること。

macOSの場合

setwd("/Users/xxx/bulksem/") #xxxは任意のユーザー名に変更すること。

featureCountsの出力ファイルからカウントデータを抽出したファイルを読み込みます。\

data <- read.table("./featureCounts_all_output.txt", header = TRUE, row.names = 1, sep = "\t")

続いて、GTFファイルからgene lengthを抽出します。featureCounts_all_output.txtを出力した際に使用したGTFファイルのバージョンがGRCh38.101だったため、そのバージョンのGTFファイルを指定しています。バージョンが違うとエラーになります。 下記スクリプトの実行は少し時間がかかるため、、本セミナーでは実行しません。今後の参考にして下さい。

library(GenomicFeatures)
txdb <- makeTxDbFromGFF("./Homo_sapiens.GRCh38.101.gtf", format = "gtf")
exons.list.per.gene <- exonsBy(txdb, by = 'gene', use.names = FALSE) #featureCountsでtranscriptレベルでカウントした場合は、by引数の'gene'を'tx'に、use.names引数をTRUEに変更します。
exonic.gene.sizes <- lapply(exons.list.per.gene, function(x){sum(width(reduce(x)))})
exonic.gene.sizes.2 <- as.numeric(exonic.gene.sizes)
names(exonic.gene.sizes.2) <- names(exonic.gene.sizes)
exonic.gene.sizes.2 <- data.frame(as.matrix(exonic.gene.sizes.2))
exonic.gene.sizes.2$ensembl_gene_id <- row.names(exonic.gene.sizes.2)

今回は事前に用意したファイルを読み込んで使用します。

下記コマンドをCygwin(Windows 11)またはターミナル(macOS)で実行してファイルをダウンロードします。

cd ~/bulksem
curl -OL https://github.com/nojima-q/2021-12-13-15_PBL_analysis/raw/main/exonic_gene_sizes_2_GRCh38.101.rds

RStudioに戻って

exonic.gene.sizes.2 <- readRDS("./exonic_gene_sizes_2_GRCh38.101.rds")
gene.len <- exonic.gene.sizes.2$as.matrix.exonic.gene.sizes.2.
names(gene.len) <- exonic.gene.sizes.2$ensembl_gene_id
gene.list.order <- row.names(data)
gene.len.sorted <- gene.len[gene.list.order]
tpm <- function(count, gene.len) {
  rate <- count / gene.len
  rate / sum(rate) * 1e6
}

カウントデータは、Transcript Per Million(TPM)値に変換します。TPMは、まずカウントデータを遺伝子長の長さで補正し、総リード数を100万リードにした値です。他にFPKMなどあり、こちらは先に総リード数を100万にしてから遺伝子長で補正します。したがって、最終的にサンプルごとに総リード数が異なることになるためサンプル間で比較する際はTPMの方が適しています。このような理由から最近ではTPM値を用いることが推奨されています。

tpms <- apply(data, 2, function(x) tpm(count = x, gene.len = gene.len.sorted))

サンプルの総リード数が１００万になっているか確認します。

colSums(tpms)

TPM値に1を足してlog2変換します。

tpms <- data.frame(log2(tpms + 1))

8 DEGsの検出

発現変動遺伝子（Differentially expressed genes; DEGs）の抽出を行います。\ 今回の公共データは患者20人、健常者19人からなるデータです。

disease <- as.matrix(tpms[,1:20])
control <- as.matrix(tpms[,21:39])
r1 <- matrix(nrow = nrow(disease))
for(i in 1:nrow(disease)){
  r1[i,] <- t.test(disease[i,], control[i,], var.equal=F, paired=F, alternative = "two.sided")$p.value
}
wel <- data.frame(tpms, r1)
wel <- na.omit(wel)
library(ggplot2)
ggplot(wel, aes(x = r1)) + geom_histogram(position = "dodge", alpha = 1, binwidth = 0.01)

P値のヒストグラムから、小さなP値の比率が高いことから、患者・健常者間で発現差異のある遺伝子は多いと考えられます。\ Rplot02 Benjamini and Hochberg法によるP値の補正

q.value <- data.frame(p.adjust(p = wel$r1, method = "BH")) #P値の補正
wel <- data.frame(wel, q.value) #元データと結合

Fold changeの算出

disease.mean <- rowMeans(wel[,1:20]) #患者サンプル
control.mean <- rowMeans(wel[,21:39]) #健常者サンプル
wel$FC <- disease.mean - control.mean #Fold changeの算出

基本的統計量が全て算出されましたので、DEGsを定義します。\ 今回は、|FC| ≥ 1かつadjusted P-value < 0.05を満たす遺伝子群をDEGsとします。

wel$Color = ifelse(abs(wel$FC) >= 1 & wel$p.adjust.p...wel.r1..method....BH.. < 0.05, "DEG", "non-DEG")

Volcano plotでDEGsを可視化します。

ggplot(wel, aes(x=FC, y=-log10(p.adjust.p...wel.r1..method....BH..), colour=Color)) + theme(panel.background = element_blank(), axis.line=element_line(colour = "black"), panel.grid = element_line(colour = "gray")) +
  geom_point(alpha=1, size=3) + xlab(expression(log[2]("Disease / Control"))) + ylab(expression(-log[10]("adjusted P-value"))) +
  theme(legend.position = "right", plot.title = element_text(size = rel(1.5), hjust = 0.5), axis.title = element_text(size = rel(1.25))) +
  geom_hline(yintercept=1.30102999566, linetype="dashed", colour="#00ff00", size = 1) +
  geom_vline(xintercept=1, linetype="dashed", colour="#00ff00", size = 1) +
  geom_vline(xintercept=-1, linetype="dashed", colour="#00ff00", size = 1) +
  theme(legend.title=element_text(size=30, face = "plain") , legend.text=element_text(size=30, face = "plain")) +
  theme(plot.title=element_text(size=30, colour="black", face = "bold",hjust = 1.0), axis.text=element_text(size=25, colour="black", face = "plain"), axis.title=element_text(size=30, face = "plain"))

Rplot01 赤が|FC| ≥ 1かつadjusted P-value < 0.05を満たす遺伝子群、青が満たさなかった遺伝子群です。\ 縦の緑破線がFCの閾値、横の緑破線がadjusted P-valueの閾値を示します。

9 エンリッチメント解析

疾患データでどのような生物学的イベントが起きているのか、どのようなパスウェイが動いているのか調査します。\ DEGsを入力データとしてGene Ontology(GO)解析とKEGGパスウェイ解析について紹介します。\ 今回の解析方法では、入力データとしてFC値とEntrez gene IDが必要です。\ そこで、biomaRtライブラリーを使ってEnsembl gene IDからEntrez gene IDに変換します。

library(biomaRt)
DEG <- data.frame(row.names(wel[wel$Color == "DEG",]), wel[wel$Color == "DEG",]$FC)
colnames(DEG) <- c("ensembl_gene_id", "FC")
db <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
hg <- useDataset("hsapiens_gene_ensembl", mart = db)
res.gene <- getBM(attributes = c("ensembl_gene_id", "entrezgene_id"), filters = "ensembl_gene_id", values = DEG, mart = hg, uniqueRows=T)
res.gene <- na.omit(res.gene)
DEG2 <- merge(res.gene, DEG, by = "ensembl_gene_id")

上記４、５行はbiomaRt側のサーバーが停止しているとエラーとなります。その場合は事前に用意したRDataを読み込んで使用して下さい。

load(file = "./biomaRt_104.RData")

9-1 GO解析

GO解析では、Biological Process(BP)、Molecular Function(MF)、Cellular Component(CC)の３種類があります。\ 下記はBPを選択していますが、ont引数を変更すればMF、CCを選択可能です。ALLを指定すると全て解析対象となります。棒グラフやネットワークグラフなど様々な方法で描写することが可能です。showCategory引数で表示するterm数を指定できます。\

library(clusterProfiler)
library(org.Hs.eg.db)
#GOエンリッチメント解析を実行
ego.result <- enrichGO(gene = as.character(DEG2$entrezgene_id),
                                        universe      = NULL,
                                        OrgDb         = org.Hs.eg.db,
                                        ont           = "BP", #"BP","CC","MF","ALL"から選択
                                        pAdjustMethod = "BH",
                                        pvalueCutoff  = 0.05,
                                        qvalueCutoff  = 0.05, 
                                        readable      = TRUE) #Gene IDを遺伝子名に変換

ego.result.simple <- simplify(ego.result) #GO termの冗長性を除去
barplot(ego.result.simple, drop=TRUE, showCategory=30)
dotplot(ego.result.simple)
d <- GOSemSim::godata("org.Hs.eg.db", ont = "BP")    
ego.result.simple <- enrichplot::pairwise_termsim(ego.result.simple, semData = d)
emapplot(ego.result.simple)
goplot(ego.result.simple)

9-2 KEGGパスウェイ解析

KEGG(Kyoto Encyclopedia of Genes and Genomes)パスウェイ解析を行います。\ GO解析と同様に様々な方法で描写可能です。showCategory引数で表示するterm数を指定できます。

ego.result.kegg <- enrichKEGG(gene = as.character(DEG2$entrezgene_id),
                                               universe      = NULL,
                                               organism = "hsa",
                                               keyType = "ncbi-geneid",
                                               pAdjustMethod = "BH",
                                               pvalueCutoff  = 0.05,
                                               qvalueCutoff  = 0.05)

barplot(ego.result.kegg, showCategory=10)
dotplot(ego.result.kegg, showCategory=10)   
ego.result.kegg <- enrichplot::pairwise_termsim(ego.result.kegg)
emapplot(ego.result.kegg, showCategory=10)

9-3 Disease enrichment解析

Disease enrichment解析では、入力した遺伝子リストがどのような疾患に関わっているか調査することができます。

library(DOSE)
DEG3 <- DEG2[abs(DEG2$FC) >= 2,]
geneList2 <- DEG3[,3]
names(geneList2) <- DEG3$entrezgene_id
edo2 <- enrichDGN(DEG3$entrezgene_id)
barplot(edo2, showCategory=20, font.size = 12)
edox2 <- setReadable(edo2, 'org.Hs.eg.db', 'ENTREZID')
cnetplot(edox2, foldChange=geneList2, showCategory = 11, layout = "kk", fixed = TRUE, categorySize = "geneNum", max.overlaps = 1000)

図の保存方法

右下に図が表示されたのち、下図赤枠のExportをクリック

Save a Imageをクリック

任意の保存先、フォーマット、ファイル名を指定して、Saveをクリック

コードをRスクリプトとして保存

タブからFile → Saveをクリックし、エンコーディングの方法を選択（デフォルトのUTF-8で構わない）しOKをクリック。\ 任意のファイル名を入力し、拡張子（ファイル名以降の部分）は.RとしてSaveをクリックし保存。（例:20230517.R）

今回使用したデータセット

元論文

特発性肺線維症(Idiopathic pulmonary fibrosis; IPF)のデータを使用しました。\ Disease enrichment解析ではPulmonary FibrosisやIdiopathic Pulmonary Fibrosisといったタームがエンリッチしています。

再解析論文

yosui-nojima / bulk_RNA-Seq_Kochi_Univ._2024-0305-0306

readme

高知大学 R5年度 DS/DXセミナー

統計基礎、bulk RNA-Seq解析ハンズオントレーニング

日時

場所

講義の概要と学習目標

用意するもの

スケジュール

１日目

２日目

統計基礎講義資料

bulk RNA-Seq解析ハンズオントレーニング

目次

Windows 11

macOS

使用するデータ(これ以降は一部を除いてWindows 11, macOS間共通)

1. 環境構築

Windows 11

- Cygwinのインストール

- 作業ディレクトリの作成

- NCBI SRA Toolkitのインストール

- FastQCのインストール

- javaのインストール

- Trimmomaticのインストール

- HISAT2のインストールとコンパイル

- featureCountsのインストール

- PATHの追加

- Rのインストール

データ解析の再現性

Rについて

- RStudioのインストール

RStudioによるRの起動

- Rパッケージのインストール

macOS

ターミナルの起動

- 作業ディレクトリの作成

- NCBI SRA Toolkitのインストール

- FastQCのインストール

- Trimmomaticのインストール

- HISAT2のインストール

- featureCountsのインストール

- PATHの追加

2. 使用データ

3. 公共データベースの紹介

NCBI SRA (https://www.ncbi.nlm.nih.gov/sra)

DDBJ Sequence Read Archive (https://www.ddbj.nig.ac.jp/dra/index.html)

EBI ENA (https://www.ebi.ac.uk/ena/browser/home)

※データ検索は、スクリーン上で実演します。

データ取得方法

4. FASTQファイルの形式、クオリティーコントロールについて

１．FASTQファイルとは

２．FASTQファイルのクオリティーコントロール

3. FastQC解析結果

Per base sequence quality

Per tile sequence quality

Per sequence quality scores

Per base sequence content

Per sequence GC content

Per base N content

Sequence Length Distribution

Sequence Duplication Levels

Overrepresented sequences

Adapter Content

4. アダプターの除去およびリードのトリミング

awkコマンドを使ったバッチ処理

リードトリミング後のFastQC結果

Per base sequence quality

Per tile sequence quality

Per sequence quality scores

Per base sequence content

Per sequence GC content

Per base N content

Sequence Length Distribution

Sequence Duplication Levels

Overrepresented sequences

Adapter Content

4 リファレンスゲノムファイル、アノテーションファイルの取得

5 リファレンスゲノムへのマッピング

5-1 リファレンスゲノムファイルのインデックス化

高知大学　R5年度　DS/DXセミナー　