2018-06-13

バクテリアのRNA seq定量ツール EDGE-pro

　バクテリアゲノム中の遺伝子の発現を測定することは、感染の治療法の開発から合成ゲノムの作成まで、非常に幅広い用途を有する。バクテリアにおける遺伝子発現研究は、代謝経路を研究し、変異株の特性を同定し、他の点ではそれらのゲノムにおける分子過程をよりよく理解するために用いられてきた（論文より ref1.2)。

　10年以上にわたって、マイクロアレイは遺伝子発現を研究するための主な手段であった。しかしながら、マイクロアレイ技術は、プローブが設計された転写産物のみを捕捉するので、したがって、よく研究された種の既知の株における既知の遺伝子へ適用を制限する。あるいは、定量的PCR（qPCR）の使用は、ゲノム中のすべての遺伝子ではなく、特定の遺伝子を定量することを可能にするが、この技術はゲノム全体のスケールではるかに高価である。第2世代シーケンシングの品質、効率、およびコストの近年の改良によって、マイクロアレイ解析を置き換え、RNAを直接捕捉してシーケンスするいわゆるRNA-seq実験の爆発がもたらされた。マイクロアレイ実験では、参照株と新規の株との違いによって、マイクロアレイ上のいくつかのプローブとのハイブリダイゼーションが妨げられる可能性がある。対照的に、RNA-seq実験では、転写された遺伝子が配列決定され、ゲノムにアライメントされる。アラインメントアルゴリズムは、多くのミスマッチを許容することができ、それによって標的ゲノムが参照から分岐した場合でも遺伝子発現の高感度測定が可能になる。さらに、全転写物が配列決定されるので、RNA-seqデータはゲノムのオペロン構造も明らかにする。

　RNA-seq技術の導入以来、サンプル中の遺伝子発現を定量化し、複数の試料間の遺伝子発現の差異を見出すためにソフトウェア方法が開発されている(ref.5-9）。しかし、遺伝子発現を推定するためのツールは、真核生物遺伝子の構造を同定する目的で開発された。これらのツールは、遺伝子内のイントロン領域を見いだし、高等真核生物で一般的な選択的スプライスバリアントを見出す努力の多くに焦点を当てている。逆に、バクテリア遺伝子はイントロンを有さず、選択的にスプライシングされない。したがって、それらの転写産物を分析する際にスプライスバリアントを探す必要はない。

　バクテリアゲノムは非常に密集しており、その多くは互いに重なっている。これまでのRNA-seqソフトウェア法は、一般に、ヒトや他の哺乳動物（ほとんどのRNA-seq実験の主な標的）においては重複が極めてまれであるため、重複する遺伝子を扱う手段を提供していない。これとは対照的に、典型的なバクテリアゲノムの約90％はタンパク質をコードしている（ref.10）。広い進化範囲にまたがる220の原核生物ゲノムの研究（ref.10）から、これらの種の全遺伝子の29％は5 'または3'で互いにオーバーラップがあることが示されている。これらのオーバーラップは、わずか数塩基対（bp）から100bpをはるかに超える範囲にわたる。重複する遺伝子は、同じ鎖上または反対の鎖上に存在し得る; 従って、ストランド特異的RNA-seqは、せいぜい部分的な解を提供するだけである。オーバーラップ領域内にマップされたRNAシーケンスリードの場合、2つの遺伝子のどちらがリードをもたらしたかを決定することが不可能であり、したがって、各原核生物遺伝子の遺伝子発現を決定する際に課題が生じる。

　さらに、分析のための要求を複雑にするバクテリアのRNA-seq技術は、ポリアデニル化不在のために、真核生物RNA seqプロトコルと少なくとも1つの大きな相違点を有することである。真核生物の転写産物上の長いポリAテイルは捕捉するプローブとして使用できるが、バクテリアのRNA-seqではランダムプライミング法が必要となる。もう一つの課題は、捕捉された細菌転写産物の80％がrRNAになることである。 rRNAを除去するための方法が開発されているが、RNA-seq実験ではまだ多くのrRNAが依然として出現する。

　真核生物ゲノムとバクテリアゲノムとの間の相違、およびRNA-seqプロトコル間の相違のために、発現解析のための既存のプログラムは、バクテリアゲノムからRNA-seqデータに適用されると、ほとんど機能しないか完全にbreak down する。従って、バクテリア RNA-seqデータにおける遺伝子発現のレベルを推定するためには、新しいバイオインフォマティクス法が必要である。現在のところ（論文執筆時点）、この目的のための独立型ツールは存在しない。複数のバクテリア RNA-seqプロジェクトが公開されているが、これらは遺伝子発現を定量化するためにad hoc(その場限り)な方法を用いており、Bowtie、MAQ、SOAP、BWA、ELAND、Novoalign、などの次世代配列アライナーを使用して、最初に入力ゲノムにリードをアライメントする。これらのアライメントから、各遺伝子にマッピングされたリード数をカウントし、通常は各遺伝子の長さでカウントを正規化する。これらのad hocなアプローチのいくつかのレビューは、Guell et al (ref.27) やVan Vliet.(ref.28)に見つかる。

　標準的なアラインメント手法が直面する課題の1つは、すべてのデータセットにおいて、いくつかのリードがゲノムの複数の場所にアライメントされうることである。これらの複数の場所のどれが真のシーケンス源であるかを決定することは困難であり、時には不可能である場合がある。特に、リピートがある場合はそうである。この問題を回避するために、いくつかの以前の方法では、マルチアラインされたリードを単に破棄する。この戦略は、リピートを含む遺伝子の見かけの発現レベルを有意に（そして間違って）減少させる可能性がある。より深刻な問題は、ファミリー内のどのリードもファミリー遺伝子のすべてのコピーに等しくマッピングされるような遺伝子ファミリーで生じる。マルチアライメントリードをカウントする他の方法では、リードがマップされている各場所に分数を割り当てたり、リードをマルチマップされた場所のいずれかにランダムにアサインしたりする。これらのどれも完全な解決法を提供していないが、我々（この論文の著者ら）が示すように、分数のリードカウントの使用は実際にはうまく機能する。

　バクテリアの遺伝子発現レベルを推定するために現在使用されている臨時の方法は、合理的にうまくいくが、しばしば使用するのは容易ではない。一部のユーザーは、複数のソフトウェアツールを連続して実行する必要があり、あるプログラムの出力は、次のツールの入力が正しくない場合があり、データを再フォーマットするために付随するプログラムが必要になる。本論文では、原核生物ゲノムの遺伝子発現レベルを推定するために特別に設計された最初のスタンドアローン方式であるEDGE-pro（Estimated Degree of Gene Expression in PROkaryots）を紹介する。 EDGE-proは、上記の課題に対する解決策を提供する効率的なソフトウェアシステムである。

　EDGE-proパイプラインは、リファレンスゲノム、そのゲノム内のタンパク質コード遺伝子の座標を含むタンパク質翻訳テーブル（ptt）、tRNAおよびrRNA遺伝子の座標を含む別のテーブル（rnt）、およびRNA-seq 自分自身を読み込む。 pttとrntのテーブルファイルが利用できない場合は、Glimmerを実行してタンパク質コード遺伝子を検索し、tRNAscan-SEとRNAmmerをRNA遺伝子に変換することで、ゲノム配列から別々に生成することもできる。EDGE-proパイプラインは、以下の4つの手順で構成されている。

メインマッピングステップ。すべてのリードがリファレンスゲノムにアライメントされる。
マルチアライメントリードのフィルタリングステップ。
リファレンスゲノムにおける各塩基のカバレッジデプスの計算
カバレッジデプスを、各遺伝子のRPKM（マッピングされた百万単位当たりの遺伝子キロベース当たりのリード量）に変換する。

各ステップの詳細は論文で確認してください（リンク）。このツールはedgeRと直接関係があるわけではありません。このツールは、内部でbowtieを動かしマッピングして、フィルタリングした後にリードカウントを出すプロセスを自動で行うツールです。

公式ページ

https://ccb.jhu.edu/software/EDGE-pro/index.shtml

マニュアル

https://ccb.jhu.edu/software/EDGE-pro/MANUAL

インストール

cent os6でテストした。

依存

Bowtie2

本体中に、linux 64bit向けのコンパイル済みのBowtie2があるので、Bowtie2にパスが通ってなければそれを使ってください。

EDGE-pro本体のダウンロード (linux向けバイナリリンク)。

wget https://ccb.jhu.edu/software/EDGE-pro/EDGE_pro_v1.3.1.tar.gz
tar -vxf EDGE_pro_v1.3.1.tar.gz

perl EDGE_pro_v1.3.1/edge.pl -h

Usage

edge.pl <-g genome> <-p ptt> <-r rnt> <-u reads>

# [OMP_NUM_THREADS=n] PATH/edge.pl <-g genome> <-p ptt> <-r rnt> <-u reads> [-v reads2, if paired-end] [-m minInsert] [-M maxInsert] [-t bowtie threads] [-s path to Source code] [-o prefix of Output files] [-w window] [-i initial transcription site window] [-x similarity] [-l read length] [-c minimum coverage] [-n count type: 0-partial count (default), 1-random count] [-h]

(anaconda3-5.0.0) [uesaka@cyano EDGE_pro_v1.3.1]$ perl edge.pl -h

USAGE:

------

[OMP_NUM_THREADS=n] PATH/edge.perl <-g genome> <-p ptt> <-r rnt> <-u reads> [options]

MANDATORY FILES:

----------------

-g genome: fasta file containing bacterial genome. If multiple chromosomes/plasmids exist, they must be combined into one file before running EDGE (see MANUAL)

-p ptt: ptt file with coordinates of coding genes, in Genbank format. If multiple chromosomes/plasmids exist, they must be combined into one file (see MANUAL).

-r rnt: rnt file with coordinates of rRNAs and tRNAs, in Genbank format. If multiple chromosomes/plasmids exist, they must be combined into one file (see MANUAL).

-u reads: fastq file of reads. If multiple fastq files exists, see MANUAL to combine them, if possible.

OPTIONAL FILES/PARAMETERS:

--------------------------

-v reads2: fastq file of mates in paired-end data. If this file is not entered, single-end reads are assumed.

-m min: min is an integer value. It is minimum insert size in paired-end library. Default: 0.

-M max: max is an integer value. It is maximum insert size in paired-end library. Default: 500.

-t threads: threads is an integer value. It is the number of threads to be used by Bowtie2. Default: 1.

OMP_NUM_THREADS is an integer environmental optional parameter that specifies the number of threads to be used to count per base coverage. Note that it is entered before the command ./edge.perl. Default: 16.

-s source_dir: It is a string specifying the absolute of relative path to the directory that contains all scripts. Default: working directory.

-o prefix: It is a string specifying the prefix of all output files. Default: out.

-w window: It is an integer specifying the window size close to overlapping region used to approximate the coverage of a gene close to the overlapping region in order to distrbute the coverage of the overlapping region between two overlapping genes. Default: 100.

-i untranslated region: It is an integer specifying the window size of the untranslated region bewteen the initial transcription site and the start codon. Default: 40.

-x similarity: It is a decimal number spcifying the percentage used to determine when two coverage values are considered similar. For example, if the similarity is x, and coverage of a region is C, then another region is considered similarly expressed if its coverage is in the interval [(1-x)*C,(1+x)*C]. Default: 0.15.

-l read length: It is an integer specifying the read length. If read length is not specofied, the first 1000 reads are used to approximate the read length.

-c minimum coverage: It is an integer specifying the minimum average coverage of gene for gene to be considered expressed. Coverage lower than specified is assumed to be noise and gene is considered to not be expressed. Default: 3.

-n count type: It is 0 or 1 specifying how to count reads that map to multiple places. 0 denotes giving a partial count to each place where the read maps. 1 denotes picking randomly one of the places where the read maps and assigning full count to that one place. Default: 0.

-h: display this help

ラン

Role of CRP K100 positive charge on Escherichia coli global transcriptomeのWTのtripicates（SRR5416993 - SRR5416995）と K100 variant （SRR5416999 - SRR5417001）を使い、２群間比較を行う。SRAにDepositされたデータは３グループあるので、３グループの多群間比較も行える。目的により変わるが、ここではクオリティトリミング後もシーケンスデータは十分量あるものと仮定し、解析は盲目的に最後まで実行する（検証ステップはOptional stepに記載）。

１、シーケンスデータのダウンロード

以下のURLのSRAデータを使用する。E.coliのWTと変異株のRNA seq解析データとなる。

https://www.ncbi.nlm.nih.gov/bioproject/PRJNA381697

#ヒアドキュメントでダウンロードリスト作成
cat >sra_run_ids.txt <<EOF 
SRR5416993 
SRR5416994 
SRR5416995 
SRR5416999 
SRR5417000 
SRR5417001
EOF

#ダウンロード(srapyを使う（Github）)
get-run.py -s -f sra_run_ids.txt -d sra_files/

#またはsra-toolkitを使う
prefetch --option-file sra_run_ids.txt
mv ~/ncbi/public/sra/SRR541* sra_files/

２、ペアエンドのfastqへの変換

forループで処理。pfastq-dump（紹介）を使えばより早く終わります。

mkdir fastq

for srr in $(cat sra_run_ids.txt) 
do
 fastq-dump --split-files --skip-technical --readids \
 --defline-seq '@$sn/$ri' --defline-qual '+' sra_files/${srr}.sra \
 -O fastq/
done

３、クオリティフィルタリング

高速なfastp（紹介）を使い、アダプターとクオリティトリミングしてfastq.gz出力。forループで処理。

for srr in $(cat sra_run_ids.txt)
do 
 fastp -i fastq/${srr}_1.fastq -I fastq/${srr}_2.fastq \
 -o fastq/${srr}_clean1.fastq.gz -O fastq/${srr}_clean2.fastq.gz \
 -h fastq/${srr}_report.html -j fastq/${srr}_report.json
done

htmlとjsonのレポートもそれぞれ出している。統合した１つのレポート出したいなら、afterQC（紹介）などMultiQC（紹介）が対応したツールに切り替えて、MulitiQCに渡すのが楽ですfastp実行時に"-w 8"フラグ立てるとエラーが起きたのでdefaultの"-w 3"で実行した。

４、コンカテネート

EDGE-proはペアエンド情報を使わないので、EDGE-pro入力前にペアエンドをコンカテネートする（間違いが発生しないよう丁寧に書いたが、今回のケースではfastq-dumpのコマンドを--split-filesを付けずに実行してOK）。EDGE-proはfq.gzを扱えないのでgz解凍も行う。

for srr in $(cat sra_run_ids.txt)
do 
 cat fastq/${srr}_clean2.fastq.gz fastq/${srr}_clean2.fastq.gz > fastq/${srr}_merge.fq.gz
 sleep 1s
 gzip -dv fastq/${srr}_merge.fq.gz
done

論文より　 However, because bacterial RNA-seq analysis does not need to link together exons across splice junctions, paired-end sequencing does not provide as much of an advantage.

Optional1

BBsketch（紹介）でシーケンスデータを１つ分析。

sendsketch.sh in=SRR5416995.fastq.gz reads=100k

sendsketch.sh in=SRR5416995.fastq.gz reads=100k

Adding /home/uesaka/.linuxbrew/Cellar/bbtools/37.77/resources/blacklist_refseq_species_300.sketch to blacklist.

Loaded 1 sketch in 0.738 seconds.

Query: 1/1 DB: RefSeq SketchLen: 10032 Seqs: 200000 Bases: 27673848 gSize: 3391752 Quality: 0.9365 File: SRR5416995.fastq.gz

WKID KID ANI Complt Contam Matches Unique noHit TaxID gSize gSeqs taxName

47.53% 36.64% 97.51% 72.88% 5.45% 4768 0 4717 1403831 4414395 1 Escherichia coli str. K-12 substr. MC4100

43.88% 30.91% 97.24% 64.02% 9.10% 4402 0 4717 457401 4800053 8 Escherichia sp. 4_1_40B

43.98% 30.00% 97.25% 62.08% 9.00% 4412 0 4717 2082622 4975840 193 Escherichia sp. R7

43.99% 29.79% 97.25% 61.64% 8.99% 4413 0 4717 2082619 5014011 219 Escherichia sp. R4

42.40% 30.81% 97.13% 64.98% 10.58% 4254 0 4717 457400 4659828 47 Escherichia sp. 1_1_43

44.06% 29.39% 97.26% 60.75% 8.92% 4420 0 4717 1806490 5074112 170 Achromobacter sp. ATCC35328

K12のゲノム配列をリファレンスとして使って問題なさそう。

Optional2

統合レポート作成（multiqc紹介）。

fastqc fastq/*merge.fq
multiqc fastq/

html出力。極端におかしなサンプルはない。

f:id:kazumaxneo:20180613131848j:plain

step１−４でfastqのpreprocessingは完了。続いてmappingとリードカウントを行う。

５、リファレンス配列とコード領域情報の入手

Genome2D（HP）のbacteira（リンク）のMG1655（リンク）から、リファレンスゲノム（.fna）、プロテインテーブル(.ptt)、tRNA＆rRNAテーブル（.rnt）をダウンロード（リンク）。

wget http://genome2d.molgenrug.nl/Bacteria/Escherichia_coli_K_12_substr__MG1655_uid57779/NC_000913.fna 
wget http://genome2d.molgenrug.nl/Bacteria/Escherichia_coli_K_12_substr__MG1655_uid57779/NC_000913.ptt 
wget http://genome2d.molgenrug.nl/Bacteria/Escherichia_coli_K_12_substr__MG1655_uid57779/NC_000913.rnt

> head NC_000913.rnt

head NC_000913.rnt

Escherichia coli str. K-12 substr. MG1655, complete genome. - 1..4641652

178 RNAs

Location Strand Length PID Gene Synonym Code COG Product

16952..17006 + 55 556503834 sokC b4413 - - -

77367..77593 + 227 556503834 sgrS b4577 - - -

189712..189847 + 136 556503834 tff b4414 - - -

223771..225312 + 1542 556503834 rrsH b0201 - - 16S ribosomal RNA of rrnH operon

225381..225457 + 77 556503834 ileV b0202 - - Ile tRNA

225500..225575 + 76 556503834 alaV b0203 - - Ala tRNA

225759..228662 + 2904 556503834 rrlH b0204 - - 23S ribosomal RNA of rrnH operon

> head NC_000913.ptt

$ head NC_000913.ptt

Escherichia coli str. K-12 substr. MG1655, complete genome. - 1..4641652

4140 proteins

Location Strand Length PID Gene Synonym Code COG Product

190..255 + 21 16127995 thrL b0001 - - thr operon leader peptide

337..2799 + 820 16127996 thrA b0002 - COG0527E,COG0527E Bifunctional aspartokinase/homoserine dehydrogenase 1

2801..3733 + 310 16127997 thrB b0003 - COG0083E,COG0083E homoserine kinase

3734..5020 + 428 16127998 thrC b0004 - COG0498E,COG0498E L-threonine synthase

5234..5530 + 98 16127999 yaaX b0005 - - DUF2502 family putative periplasmic protein

5683..6459 - 258 16128000 yaaA b0006 - COG3022S,COG3022S peroxide resistance protein, lowers intracellular iron

6529..7959 - 476 16128001 yaaJ b0007 - COG1115E,COG1115E putative transporter

６、EDGE-proによるマッピングと定量

EDGE-proでマッピングを実行。forループで処理。不安なら１つずつ。

#最初にダウンロードしたEDGE-proのプログラムディレクトリ(EDGE_pro_v1.3.1)がカレントにあるものとする。
#出力ディレクトリ
mkdir read_count

#edge.plを実行。
for srr in $(cat sra_run_ids.txt)
do 
 
 perl EDGE_pro_v1.3.1/edge.pl -g NC_000913.fna -p NC_000913.ptt -r NC_000913.rnt \
 -u fastq/${srr}_merge.fq -o read_count/${srr}.out \
 -s EDGE_pro_v1.3.1/ -t 20 -c 3 
done

確認できないようなハイスピードで多量のログが標準出力されるのでトラブルではないかと心配になりますが、エラーコードが出てこなければ大丈夫です。

１サンプルの出力を確認。

$ head SRR5416993.out.rpkm_0

gene_ID start_coord end_coord average_cov #reads RPKM

b0001 190 255 40.5 19 84

b0002 337 2799 3.5 61 7

b0003 2801 3733 16.3 109 34

b0004 3734 5020 13.0 120 27

b0005 5234 5530 11.3 24 23

b0006 5683 6459 35.6 198 74

b0007 6529 7959 19.8 203 41

b0008 8238 9191 1053.5 7179 2185

空行のほか、大量にスペースができてしまう。次のステップで消す。

７、フォーマット変更

DESeq2のカウントテーブル形式に変換するスクリプトを走らせる。

cd read_count/

perl ../EDGE_pro_v1.3.1/additionalScripts/edgeToDeseq.perl \
SRR5416993.out.rpkm_0 SRR5416994.out.rpkm_0 SRR5416995.out.rpkm_0 \
SRR5416999.out.rpkm_0 SRR5417000.out.rpkm_0 SRR5417001.out.rpkm_0

perl スクリプト処理でエラーが起きたので、手作業で修正する。DESeqやedgeRで利用するリードカウントを抽出する。rpkmはいらない。

大量のスペース、余分な列（２列目）、そして空行があるのでsedとgrepで消す。そしてリードカウントのカラムだけ抽出し、各サンプルのカラムを結合する。rpkmを使うなら、5カラムでなく6カラム目を抽出する。ループにし忘れましたが、できればループ処理してください。

sed -e 's/s+/	/g' SRR5416993.out.rpkm_0 |cut -f 1 | grep -v '^$' > name
sed -e 's/s+/	/g' SRR5416993.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5416993_read_count
sed -e 's/s+/	/g' SRR5416994.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5416994_read_count
sed -e 's/s+/	/g' SRR5416995.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5416995_read_count
sed -e 's/s+/	/g' SRR5416999.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5416999_read_count
sed -e 's/s+/	/g' SRR5417000.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5417000_read_count
sed -e 's/s+/	/g' SRR5417001.out.rpkm_0 |cut -f 5 | grep -v '^$' > SRR5417001_read_count

#出力ファイルのカラムを結合
paste name SRR5416993_read_count SRR5416994_read_count SRR5416995_read_count SRR5416999_read_count SRR5417000_read_count SRR5417001_read_count > count.txt

#vi count.txt #１列目を修正。以下のようにした。
> head -n count.txt

original count.txt

$ head -n 1 count.txt

gene_ID #reads #reads #reads #reads #reads #reads

modified count.txt

$ head count.txt

gene_ID WT1 WT2 WT3 K100R1 K100R2 K100R3

b0001 19 21 23 9 109 31

b0002 61 69 80 47 89 127

b0003 109 120 152 132 172 181

b0004 120 129 116 157 125 207

b0005 24 9 11 13 22 26

b0006 198 195 234 241 157 303

b0007 203 197 202 243 194 203

b0008 7179 7028 8557 7635 7808 8414

b0009 345 276 340 337 403 351

８、edgeRやDEseq2で正規化、DEGを抽出。Rで可視化

本ツールではrpkmを出しますが、現在はrpkmは推奨されていません。最近の論文（pubmed）のように、リードカウントの方を使いedgeRやDEseq2で正規化してDEGを出してください。edgeRの流れについては、手前味噌ですが、2017年に勉強会用の記事をまとめましたので、参考にしてください（シロイヌナズナのRNA seq解析）。DEseq2を使った流れについては、北大の中川先生がデータのダウンロードから分かりやすく説明されています。

https://ncrna.jp/blog/item/388-deseq2-ggplot2

また次のリンクは英語になりますが、こちらはEDGE-proのチュートリアルです。

https://bioinformatics.uconn.edu/edge-pro-tutorial/

データの取り込みからDESeq2の検定まで通して説明されています。ただしDEseq2のアップデートで一部使えないコマンドが出てきている？ようです（未確認）。

バクテリアのRNA seqを行なっている人は、情報が少なく苦労されてると思います。今更かもしれませんが、参考にしてください。

このツールにこだわらなければ、Kallistoでもいいかもしれません。

https://groups.google.com/forum/#!topic/kallisto-sleuth-users/U8KdYSCjPyE

https://groups.google.com/forum/#!topic/kallisto-sleuth-users/Yu_T4jkgSuE

https://www.biostars.org/p/238650/

Nature Biotechnology volume 34, pages 525–527 (2016)

https://www.nature.com/articles/nbt.3519

引用