【追記】
- tximport での作業箇所でコマンドラインの脱落を修正しました（2019.7.3）
- TCC-GUIの使い方を説明した動画が　生命科学系DB・ツール使い倒し系チャンネル　統合TVで公開されたのでリンクを貼りました。（2019.12.27）
- TCC-GUIをUbuntu18.04LTSでインストールする際のエラー対策はこちら（2020.2.5）

概要

salmonを用いることにより、alighnerを使用せず、転写物配列のインデックスを利用して、readsをカウントする。
TCC-GUIを用いて、PCA分析、クラスタリング、シルエット分析、これらの可視化、発現変動遺伝子(DEG)抽出を行う。
どの程度、コンピュータ資源と解析時間を節約して、DEG抽出が可能かを確認する。

環境

MacBook Pro 13-inch (Touch Bar) 2018 16GB RAM Core-i5 (4core, 8 threads)
内蔵用SSD 500GB　を USB3.1 Gen2 接続で外付けとして利用

macOS 10.14.4
anaconda3-2018.12
Python 3.7.1
R version 3.5.2
RStudio ver1.1.463

準備

1. Salmonのインストール

shell

conda create -n salmon salmon -c bioconda

今回ダウンロードされたパッケージのバージョン

package	build
jemalloc-5.1.0	h6de7cb9_1001 1.6 MB conda-forge
libboost-1.67.0	hebc422b_4 18.8 MB
salmon-0.13.1	heb0d2e1_0 13.5 MB bioconda
tbb-2019.4	h04f5b5a_0 156 KB conda-forge

起動

source activate salmon

2. 転写物の配列データの取得

１）Ensembleの場合
①マウス：Genome assembly: GRCm38.p6 (GCA_000001635.8)

shell

mkdir -p Volumes/Samsung500_J/mouse_salmon
cd /Volumes/Samsung500_J/mouse_salmon
wget ftp://ftp.ensembl.org/pub/release-95/fasta/mus_musculus/cdna/
gunzip Mus_musculus.GRCm38.cdna.all.fa.gz
grep ">" Mus_musculus.GRCm38.cdna.all.fa | cut -d' ' -f 1 | wc
    116067  116067 2561207

と出力されて、転写物：116,067個　であることがわかる。

②ヒト：Genome assembly: GRCh38.p12 (GCA_000001405.27)

shell

mkdir -p /Volumes/Samsung500_J/human_salmon
cd /Volumes/Samsung500_J/human_salmon
wget ftp://ftp.ensembl.org/pub/release-95/fasta/homo_sapiens/cdna/Homo_sapiens.GRCh38.cdna.all.fa.gz
gunzip Homo_sapiens.GRCh38.cdna.all.fa.gz
grep ">" Homo_sapiens.GRCh38.cdna.all.fa | cut -d' ' -f 1 | wc
　  187626  187626 3568377

と出力されて、転写物が　187,626個　であることがわかる。

2) GENCODEの場合
①マウス: gencode.vM20.transcripts.fa

shell

mkdir -p /Volumes/Samsung500_J/mouse_GENCODE/salmon/
cd /Volumes/Samsung500_J/mouse_GENCODE/salmon/
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_mouse/release_M20/gencode.vM20.transcripts.fa.gz
gunzip gencode.vM20.transcripts.fa.gz
grep ">" gencode.vM20.transcripts.fa | cut -d' ' -f 1 | wc
   138835  138835 16998089

転写物：138,835個　(Ensemble 116,067個より多い)

②ヒト

shell

mkdir -p /Volumes/Samsung500_J/human_GENCODE/salmon
cd /Volumes/Samsung500_J/human_GENCODE/salmon
wget ftp://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_29/gencode.v29.transcripts.fa.gz
gunzip gencode.v29.transcripts.fa.gz
grep ">" gencode.v29.transcripts.fa | cut -d' ' -f 1 | wc
　  206694  206694 24332877

転写物：206,694個 (Ensemble 187,626個より多い）

3. salmon indexの作成

構文

shell

salmon index -t /path/to/hoge.cdna.all.fa -i hoge_salmon_index -k 31

salmonのドキュメントによれば、-k オプションについては

We find that a k of 31 seems to work well for reads of 75bp or longer, but you might consider a smaller k if you plan to deal with shorter reads. Also, a shoter value of k may improve sensitivity even more when using the –validateMappings flag. So, if you are seeing a smaller mapping rate than you might expect, consider building the index with a slightly smaller k.

となっているが、salmon index -hで確認すると 31 がデフォルト値なので、現状のリード長なら省略で可。
実行時間は以下のようにそれほどかからない。

1) Emsemble の場合
①マウス

shell


cd /Volumes/Samsung500_J/mouse_salmon
salmon index -t Mus_musculus.GRCm38.cdna.all.fa -i Mus_musculus.GRCm38_salmon_index -k 31

② ヒト

shell

cd /Volumes/Samsung500_J/human_salmon
time salmon index -t Homo_sapiens.GRCh38.cdna.all.fa -i Homo_sapiens.GRCh38_salmon_index -k 31

real    3m6.119s
user    2m54.204s
sys    0m8.833s

2）GENCODE の場合
① マウス

shell

cd /Volumes/Samsung500_J/mouse_GENCODE/salmon/
time salmon index -t gencode.vM20.transcripts.fa -i gencode.vM20_salmon_index -k 31

real    2m47.658s
user    2m39.345s
sys    0m7.678s

②ヒト

shell

cd /Volumes/Samsung500_J/human_GENCODE/salmon
time salmon index -t gencode.v29.transcripts.fa -i gencode.v29_salmon_index -k 31

real    3m31.654s
user    3m17.814s
sys    0m13.055s

4. mapping-based modeでのカウント

alighnerを用いずに直接定量することを売り物としているが、実際にはsalmonはBAMを利用することもできる。
ここでは、starとの比較のためfastqからの直接定量の方法を示す。

構文 (PE/unstranded reads の場合）

shell

./bin/salmon quant -i transcripts_index -l "LIBTYPE" -1 reads1.fq -2 reads2.fq --validateMappings -o transcripts_quant

注１：anacondaでインストールした場合、./bin/ は省略
注２： -l -1 -2 ( -r :SEの場合）の順番が守られていないといけない。

Order of command-line parameters
The library type -l should be specified on the command line before the read files (i.e. the parameters to -1 and -2, or -r). This is because the contents of the library type flag is used to determine how the reads should be interpreted.

また、"LIBTYPE"の指定も重要。readが各ストランド上にどのような向きで張り付くかを考慮してカウントするので、ライブラリー調整法に合わせて設定する必要がある。今回はIU (an unstranded paired-end library where the reads face each other)を選択。

注３：-0 で指定されたoutput dirを作成し結果を保存する

実行（例：マウスGENCODE と　hoge_1.fastq/hoge_2.fastq）

shell

mkdir -p /Volumes/Samsung500_J/salmon_hoge
cd /Volumes/Samsung500_J/salmon_hoge

salmon quant -i /Volumes/Samsung500_J/mouse_GENCODE/salmon/gencode.vM20_salmon_index \
-l IU \
-1 /Volumes/Samsung500_J/fastq/hoge_1.fastq \
-2 /Volumes/Samsung500_J/fastq/hoge_2.fastq \
--validateMappings \
-o /Volumes/Samsung500_J/salmon_hoge/hoge_fastq_salmon_quant

結果　所要時間：約3分　(star-RSEMの10倍以上の速さ）

出力されたカウントデータのファイルは、quant.sf である。以下、作成されたファイルとサブディレクトリです。

shell

(salmon) $ ls -lh  /Volumes/Samsung500_J/salmon_hoge/hoge_fastq_salmon_quant
total 40576
drwxr-xr-x  0 piyo  staff     0B  4  1 16:49 aux_info
-rw-r--r--  1 piyo  staff   397B  4  1 16:48 cmd_info.json
drwxr-xr-x  0 piyo  staff     0B  4  1 16:49 libParams
-rw-r--r--  1 piyo  staff   648B  4  1 16:49 lib_format_counts.json
drwxr-xr-x  0 piyo  staff     0B  4  1 16:46 logs
-rw-r--r--  1 piyo  staff    20M  4  1 16:49 quant.sf

shell

(salmon)$ ls -lh /Volumes/Samsung500_J/salmon_hoge/hoge_fastq_salmon_quant/aux_info/
total 2048
-rw-r--r--  1 piyo  staff   695K  4  1 16:49 ambig_info.tsv
-rw-r--r--  1 piyo  staff    89B  4  1 16:49 expected_bias.gz
-rw-r--r--  1 piyo  staff   498B  4  1 16:49 fld.gz
-rw-r--r--  1 piyo  staff   1.4K  4  1 16:49 meta_info.json
-rw-r--r--  1 piyo  staff    54B  4  1 16:49 observed_bias.gz
-rw-r--r--  1 piyo  staff    54B  4  1 16:49 observed_bias_3p.gz

shell

(salmon) $ ls -lh /Volumes/Samsung500_J/salmon_hoge/hoge_fastq_salmon_quant/logs
total 128
-rw-r--r--  1 piyo  staff   3.6K  4  1 16:49 salmon_quant.log`

(salmon) $ ls -lh /Volumes/Samsung500_J/salmon_hoge/hoge_fastq_salmon_quant/libParams/
total 128
-rw-r--r--  1 piyo  staff    11K  4  1 16:49 flenDist.txt

一括処理　　　< bash スクリプトで全fastqファイルをまとめて処理>
①configファイルの作成
list_sample.txt にサンプル名に相当するprefixのリストを記す

shell

pwd
/Volumes/Samsung500_J/salmon_hoge

touch list_sample.txt
nano list_sample.txt
# nanoで以下をlist_sample.txt に記入して保存する

hoge
hoge2
hoge3
fuga
fuga2
fuga3

② script ファイルの作成

shell

touch salmon_count.sh
nano salmon_count.sh
# nanoで以下をsalmon_count.shに記入して保存する


#!/bin/bash

for sample in \`cat ./list_sample.txt\`
do
　salmon quant -i /Volumes/Samsung500_J/mouse_GENCODE/salmon/gencode.vM20_salmon_index -l IU -1 /Volumes/Samsung500_J/fastq/\${sample}\_1.fastq -2 /Volumes/Samsung500_J/fastq/\${sample}\_2.fastq --validateMappings -o /Volumes/Samsung500_J/salmon_hoge/${sample}.fastq_salmon_quant
done

③ シェルスクリプトの実行

shell

chmod +x /Volumes/Samsung500_J/salmon_hogesalmon_count.sh

sh /Volumes/Samsung500_J/salmon_hogesalmon_count.sh

5. 発現量の群間比較　ー quant.sf をtximport でファイル整形し、TCC-GUI へ流し込む

TCCそのものはコンソールで動かして使用しても構わない。しかしTCC-GUIを使用する利点は単にGUIで簡易というより、ブラウザ画面上で「グリグリ」動かせるグラフで結果を表示してくれるところにあると思う。
*おそらくPlotlyを使用しているのだろうと思う。これはJavaScriptでグラフを作成するライブラリで、PlotlyはPython、R、Julia、matlab、JavaScriptなどの言語に対応し同じ見た目のグラフを作れる。

①準備　　 RStudioでの環境設定
定量したデータ(quant.sf ファイル)を使って、edgeR や DESeq2 などで発現量の群間比較を行うことができる。この際に、Bioconductor の tximport パッケージ（Soneson et al, 2015）を利用することで、salmon の定量結果を edgeR/DESEq2 に取り込める形式にする。
　ここでは、edgeR/DESeq2を内包するTCC-GUI へ流し込み、クラスター分析、averaded shilluete score の算出、DEGの抽出、PCA分析等を行う。raw count dataが必要になることに注意。

R-3.5.2 であることを確認しておく。依存バッケージの一部に3.5.2で作成されたものがあった。
私は3.5.1から　3.5.2へ切り替えました
以下のコマンドで一旦Ｒをアンインストール

shell

sudo rm -rf /Library/Frameworks/R.framework /Applications/R.app \> /usr/local/bin/R /usr/local/bin/Rscript

R-3.5.2のMac用インストーラーでインストール
jsonlite, readrをCRANでインストール

② RStudioでの環境設定
RStudioの起動してコンソールで以下を実行していく

> getwd()
[1] "/Volumes/Samsung500_J/salmon_hoge"

tximportをインストール
ただしCRANではインストールできない。Bioconductor で検索しインストールコマンドを確認,以下をコピペし実行

if (!requireNamespace("BiocManager", quietly = TRUE))
   install.packages("BiocManager")
BiocManager::install("tximport", version = "3.8")

③ tximportでquant.sf を整形

１）directoryの構造を確認

shell

(salmon) $ ls -lh
total 168
drwxr-xr-x  0 chao  staff     0B  4  1 17:33 fuga.fastq_salmon_quant
drwxr-xr-x  0 chao  staff     0B  4  1 16:49 fuga2.fastq_salmon_quant
drwxr-xr-x  0 chao  staff   4.0K  4  1 17:24 fuga3.fastq_salmon_quant
drwxr-xr-x  0 chao  staff   4.0K  4  1 17:26 hoge.fastq_salmon_quant
drwxr-xr-x  0 chao  staff   4.0K  4  1 17:29 hoge2.fastq_salmon_quant
drwxr-xr-x  0 chao  staff   4.0K  4  1 18:08 hoge3.fastq_salmon_quant
-rw-r--r--  1 chao  staff    68B  4  1 17:14 list_sample.txt
-rwxr-xr-x  1 chao  staff   343B  4  1 17:22 salmon_count.sh

２）実行　　　コマンドライン　　directoryのPrefix を参照して以下のようにして、Salmon の定量結果を読み込んで発現量行列を取得する。

library(tximport)
library(jsonlite)
library(readr)

salmon.files <- file.path(list.files('.', pattern = 'fastq_salmon_quant'), 'quant.sf')
names(salmon.files) <- c('hoge', 'hoge2', 'hoge3', 'fuga', 'fuga2', 'fuga3')
tx.exp <- tximport(salmon.files, type = "salmon", txOut = TRUE)

【追記】Salmonの定量結果は転写産物発現量。そこで、転写産物の名前と遺伝子の名前の対応表（データフレーム型）を作成したあと、これとtximportの summarizeToGene 関数を使って、遺伝子発現量に換算する。

tx2gene <- data.frame(
    TXNAME = rownames(tx.exp$counts),
    GENEID = sapply(strsplit(rownames(tx.exp$counts), '\\.'), '[', 1)
)

3) txname|expでのcount tabelをgene id | expに変換する

salmonForTCC.gene.exp <- summarizeToGene(tx.exp, tx2gene)
summarizing abundance
summarizing counts
summarizing length

head(salmonForTCC.gene.exp$counts)
                       hoge    hoge2    hoge3     fuga    fuga2    fuga3
ENSMUST00000000001  832.000  722.000  764.000  719.000 1057.000  621.000
ENSMUST00000000003    0.000    0.000    0.000    0.000    0.000    0.000
ENSMUST00000000010    0.000    0.000    0.000    0.000    0.000    0.000
ENSMUST00000000028   36.918   13.449   33.598   17.705   41.000   21.000
ENSMUST00000000033 2081.155 1911.423 1889.926 1683.395 2561.719 1521.427
ENSMUST00000000049    2.000    3.000    7.000    6.000    3.000    4.000

salmonForTCC.gene.exp$counts をTCC-GUIへの入力とすれば良い。この段階がstar-RSEMでの出力結果に相当する。
ここで

head(salmonForTCC.gene.exp)

とすると、他にどんなデータがobjectに格納されているか確認できる、

$abundance
                           hoge    hoge2    hoge3     fuga    fuga2    fuga3
ENSMUST00000000001    19.429877    18.611211    19.115514    17.293916    20.951762    16.330796
ENSMUST00000000003     0.000000     0.000000     0.000000     0.000000     0.000000     0.000000
ENSMUST00000000010     0.000000     0.000000     0.000000     0.000000     0.000000     0.000000
・
・
・
 [ reached getOption("max.print") -- 137105 行を無視しました ] 

$counts
                         hoge    hoge2    hoge3     fuga    fuga2    fuga3
ENSMUST00000000001    832.000    722.000    764.000    719.000   1057.000    621.000
ENSMUST00000000003      0.000      0.000      0.000      0.000      0.000      0.000
ENSMUST00000000010      0.000      0.000      0.000      0.000      0.000      0.000
・
・
・
 [ reached getOption("max.print") -- 137105 行を無視しました ] 

$length
                           hoge    hoge2    hoge3     fuga    fuga2    fuga3
ENSMUST00000000001 3.084807e+03 3.086667e+03 3.069008e+03 3.071679e+03 3.086726e+03 3.070576e+03
ENSMUST00000000003 7.182677e+02 7.182677e+02 7.182677e+02 7.182677e+02 7.182677e+02 7.182677e+02
ENSMUST00000000010 2.390244e+03 2.390244e+03 2.390244e+03 2.390244e+03 2.390244e+03 2.390244e+03
・
・
・
 [ reached getOption("max.print") -- 137105 行を無視しました ] 
\$countsFromAbundance
[1] "no"

TCC-GUIに用いるのは salmonForTCC.gene.exp$counts

③ TCC-GUIへの入力のためにファイルへ書き出す

1) as.data.frame()関数を用いる

書き出す内容：　salmonForTCC.gene.exp$counts

２）実行：

table <- as.data.frame(salmonForTCC.gene.exp$counts)
write.table(table, file = "salmonForTCC.gene.txt", col.names = T, row.names = T, sep = "\t")

３）結果：

shell

$ ls -lh
total 78528
-rw-r--r--@ 1 chao  staff    10M  4  2 14:51 190402_salmon_edgeR_result.xlsx
drwxr-xr-x  0 chao  staff     0B  4  1 17:33 hoge.fastq_salmon_quant
drwxr-xr-x  0 chao  staff     0B  4  1 16:49 hoge2.fastq_salmon_quant
drwxr-xr-x  0 chao  staff   4.0K  4  1 17:24 hoge3.fastq_salmon_quant
drwxr-xr-x  0 chao  staff   4.0K  4  1 17:26 fuga.fastq_salmon_quant
drwxr-xr-x  0 chao  staff   4.0K  4  1 17:29 fuga2.fastq_salmon_quant
drwxr-xr-x  0 chao  staff   4.0K  4  1 18:08 fuga3.fastq_salmon_quant
-rw-r--r--  1 chao  staff    68B  4  1 17:14 list_sample.txt
-rw-r--r--  1 chao  staff   6.0M  4  2 15:47 salmonForTCC.gene.txt  ＃ 作成された6MBのdata
-rw-r--r--  1 chao  staff   205B  4  2 13:30 salmon_count.Rproj
-rwxr-xr-x  1 chao  staff   343B  4  1 17:22 salmon_count.sh
-rw-r--r--@ 1 chao  staff    22M  4  2 14:12 salmon_edgeR_result.txt

４）内容確認

shell

$ wc -l salmonForTCC.gene.txt 
  137272 salmonForTCC.gene.txt

$ head  salmonForTCC.gene.txt 
"hoge"    "hoge2"    "hoge3"    "fuga"    "fuga2"    "fuga3"
"ENSMUST00000000001"    832    722    764    719    1057    621
"ENSMUST00000000003"    0    0    0    0    0    0
"ENSMUST00000000010"    0    0    0    0    0    0
"ENSMUST00000000028"    36.918    13.449    33.598    17.705    41    21
"ENSMUST00000000033"    2081.155    1911.423    1889.926    1683.395    2561.719    1521.427
"ENSMUST00000000049"    2    3    7    6    3    4
"ENSMUST00000000058"    831.867    789    703.472    609.765    983.875    630.786
"ENSMUST00000000080"    1209.643    1207.307    988.747    848.888    1328.56    1063.687
"ENSMUST00000000087"    571.847    697.089    685.971    383.457    714.186    335.981

５）TCC-GUI用にヘッダー行の修正
ヘッダー行がいきなり、サンプル名で始まっている。
これはTCC-GUIでは読み込みエラーとなる。行頭にタブを打ち込んで保存。

６）現在のprojectは　salmon_count なので一旦終了

④ TCC-GUI で解析

詳細は以下を参照
https://github.com/swsoyee/TCC-GUI
生命科学系DB・ツール使い倒し系チャンネル「統合TV」でも使い方の流れを動画で確認できます。

RStudioを起動してインストール

# Check "BiocManager"
if (!requireNamespace("BiocManager", quietly = TRUE))
    install.packages("BiocManager")

# Package list
libs <- c("shiny", "shinydashboard", "shinyWidgets", "plotly", "dplyr", "DT", "heatmaply", "tidyr","utils","rmarkdown","data.table","RColorBrewer", "knitr", "cluster", "shinycssloaders", "shinyBS", "MASS", "TCC")


# Install packages if missing
for (i in libs){
  if( !is.element(i, .packages(all.available = TRUE)) ) {
     BiocManager::install(i)
  }
}

起動

  shiny::runGitHub("TCC-GUI", "swsoyee", subdir = "TCC-GUI", launch.browser = TRUE)

この起動方法はGitHubのドキュメントにあるMethod 1。ネット接続がいらないMethod ２を用いてもよい

ブラウザで起動
salmonForTCC.gene.txtを入力解析後、htmlで保存

RNA-seq 発現変動遺伝子解析：salmonとTCC-GUIを使って

概要

環境

準備

1. Salmonのインストール

2. 転写物の配列データの取得

3. salmon indexの作成

4. mapping-based modeでのカウント

5. 発現量の群間比較 ー quant.sf をtximport でファイル整形し、TCC-GUI へ流し込む

5. 発現量の群間比較　ー quant.sf をtximport でファイル整形し、TCC-GUI へ流し込む