macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

自動化されたProkaryotesのRNA seq解析パイプライン ProkSeq

2020 6/29 補足説明追加。 

 

RNA-seq技術は、導入以来、病原性細菌の研究において異なる条件にさらされた細菌からの複数のサンプルにわたる遺伝子発現の違いを同定し、定量化するために広く利用されてきた。一部の例外を除いて、遺伝子発現を評価するための現在のツールは、真核生物の遺伝子の構造を中心に設計されている。原核生物用に設計されたスタンドアローンのツールは少なく、改善が必要である。原核生物の品質管理、差分遺伝子発現の決定、ダウンストリームパスウェイ解析、極限の生物学的条件で収集されたデータの正規化のために必要なツールをすべて含む原核生物のための明確なパイプラインはまだ不足している。ここでは、原核生物用に設計されたユーザーフレンドリーで完全自動化されたRNA-seqデータ解析パイプラインであるProkSeqについて説明する。ProkSeqは、 differential gene expressionの解析、発現データの正規化、データと結果の可視化のための様々なオプションを提供し、publication-qualityの数値を生成する。
 ProkSeqはPythonで実装されており、ISCオープンソースライセンスで公開されている。ツールと詳細なユーザーマニュアルは、Docker: https://hub.docker.com/repository/docker/snandids/prokseq-v2.1、Anaconda: https://anaconda.org/snandiDS/prokseqGithub: https://github.com/snandiDS/prokseq でホストされている。

 

 原核生物RNA-seq解析は、利用可能なほとんどのRNA-seqパッケージは、入力データが原核生物のものとは多くの面で異なる真核生物の遺伝子構造を反映していることを前提としているため、困難である(Johnson, et al., 2016)。細菌の転写物はイントロンを持たず、交互にスプライスされていないため、スプライスジャンクションを考慮して開発されたアライナーを使用するとしばしばゲノム中の誤って割り当てられたリードが増加することがある(Magoc, et al.)。さらに、真核生物とは異なり、特定のストレス下では、ほぼすべての原核生物遺伝子の発現が変化する可能性がある(Creecy and Conway, 2015)。さらに、実験セットアップのばらつき、プラスミド遺伝子の存在および過剰発現、ならびにRNA-seqプロトコルの違いのために、原核生物のデータには、品質トリミング、アダプター除去、および歪んだデータの正規化がしばしば必要とされる(Magoc, et al., 2013; McClure, et al)。

 原核生物には、SPARTA(Johnson, et al., 2016)、EdgePRO(Magoc, et al., 2013)、およびRockHopper(McClure, et al., 2013)のような、RNA-seqデータの解析を容易にすることができるいくつかのソフトウェアパッケージがあるが、いずれも、データ操作に関する実質的な知識を必要とする。そこで、原核生物RNA-seqデータ解析を行う際の人間の介入を減らすために、完全自動化されたコマンドラインベースのワークフローであるProkSeqを開発した。ProkSeqは、利用可能な様々なツールとPythonで書かれた組み込み関数を統合している。ProkSeqは、RNA-seqデータの品質管理ステップから、異なる発現を持つ遺伝子のパスウェイエンリッチメント解析までを処理する。差分発現、正規化発現、可視化のための多様なオプションを提供し、出版品質の数値を生成する。人間の介入が少ないため、ProkSeqの使用は、データの再フォーマットが必要となることが多い個別のツールを逐次適用するよりも時間がかからない。利便性に加えて、ProkSeqのマルチスレッド機能により、パイプラインにかかる時間が短縮される。

 ユーザーは任意の分析を行うためにパラメータファイルの設定を柔軟に変更することができる。サンプルファイルには、解析に含まれる fastq ファイルの名前が記載されており、処理サンプルやコントロールサンプルなどの実験クラスも定義されている。ProkSeqは、まずリードの品質をFastQC(http://www.bioinformatics.babraham.ac.uk/projects/fastqc)とafterQC(Chen, et al)でチェックし、bowtie2(Langmead and Salzberg, 2012)を使用してリードをリファレンスゲノムにマッピングする。ProkSeqはその後、各ライブラリのアライメント品質に関するレポートを図とテキストの両方で生成し、カバレッジの均一性、タンパク質コーディング配列に沿った分布、5'および3'UTR領域、RseQC(Wang, et al., 2012)およびその他の組み込み機能によって生成されたリード重複率とストランド特異性に関する情報を提供する。遺伝子あたりの総リード数はfeatureCounts(Liao, et al., 2014)を用いて計算される。またProkSeqは、各遺伝子の正規化された遺伝子発現値を、100万あたりの転写物(TPM)(Wagner, et al. 2012)や100万あたりのカウント(CPM)(Wagner, et al. 2012)の形で計算する。これらを算出した計算式については、補足(S1)で説明する。さらに、バイアスを考慮したアルゴリズムを使用して転写物の発現を定量化するために、salmon(Patro et al、2017)を統合しており、これはその後のdifferential expressionの解析の精度および信頼性を実質的に向上させる。
 ProkSeqは、DEseq2(Love, et al., 2014)、edgeR(Robinson, et al., 2010)、およびNOISeq(Tarazona, et al., 2015)などのdifferential expression analysisのためのいくつかのツールを統合している。differential expression genesの下流解析には、ProkSeqではGOエンリッチメント、clusterProfilerを統合したパスウェイエンリッチメントを用いている(Yu, et al. アライメント前後の品質統計やグラフィカルな可視化に関するレポートはpdfとHTML形式で作成される。ProkSeqの重要なユニークな機能の1つは、歪んだデータに対する確立された正規化手法の統合である(Creecy and Conway, 2015; Zhu, et al)。 さらに、このパッケージは、任意のゲノムブラウザで可視化するための1ヌクレオチド分解能のwiggleファイルを生成する。ProkSeqは、データ解析の各ステップで鮮やかなグラフィックと出版可能な数値を生成し、ユーザーのデータへの信頼性と理解度を高める。その方法については、補足方法(S1)で詳しく説明している。

 

 

インストール

推奨されているdockerイメージを使ってテストした(ホストはmacos)。

本体 Github

#docker hub(link)
docker pull snandids/prokseq-v2.1:v1

#Bioconda(link)
conda install -c bioconda prokseq -y

python scripts/pipeline-v2.8.py -h

# python scripts/pipeline-v2.8.py -h

Usage: pipeline-v2.8.py [options] arg

 

Options:

  -h, --help            show this help message and exit

  -s SAMPLE_FILE_NAME, --sample=SAMPLE_FILE_NAME

                        provide the sample file

  -p PARAMETER_FILE_NAME, --param=PARAMETER_FILE_NAME

                        provide the parameter file

  -n NUMBER OF PROCESSORS, --numproc=NUMBER OF PROCESSORS

                        provide the number of processors

 

 

テストラン

ここではdockerイメージを使う。

docker run -it snandids/prokseq-v2.1:v1
source /etc/profile.d/conda.sh
conda activate py36
cd prokseq

カレントのファイル構成

ls -lh /root/prokseq

# ls -lh

total 111M

drwxr-xr-x 10 root root 4.0K May 29 06:14 Output

-rw-rw-r--  1 1000 1000  21K Jun  7 07:08 README.txt

-rw-r--r--  1 1000 1000 4.6M May  1 06:29 SequenceChromosome.fasta

drwxr-xr-x  2 root root 4.0K May 29 05:33 bowtie2_genome

-rw-r--r--  1 root root   18 May 29 06:14 chr.sizes

-rw-r--r--  1 root root 5.0K May 29 06:16 class.log

drwxr-xr-x  2 root root 4.0K May 29 06:09 data

drwxr-xr-x  1 root root 4.0K May 29 06:16 depend

-rw-rw-r--  1 1000 1000 273K May  1 06:28 oldAnnotationGFF.bed

-rw-rw-r--  1 1000 1000 367K May  1 06:28 oldAnnotationGFF.gtf

-rw-rw-r--  1 1000 1000 2.9K Jun  8 07:00 param.input.bowtie

-rw-rw-r--  1 1000 1000 3.1K Jun  8 07:00 param.input.salmon

-rw-r--r--  1 1000 1000 1.3M Mar 30 07:12 reads_1.fastq

-rw-r--r--  1 1000 1000 1.3M Mar 30 07:13 reads_2.fastq

-rw-r--r--  1 1000 1000 1.3M Mar 30 08:28 reads_3.fastq

-rw-r--r--  1 1000 1000 1.3M Mar 30 08:28 reads_4.fastq

-rw-rw-r--  1 1000 1000 8.2M May 27 13:41 sampleCtrl_1.R1.fq

-rw-rw-r--  1 1000 1000 8.6M Jun  2 05:39 sampleCtrl_1.R2.fq

-rw-rw-r--  1 1000 1000 8.2M May 27 13:41 sampleCtrl_2.R1.fq

-rw-rw-r--  1 1000 1000 8.6M May 27 13:42 sampleCtrl_2.R2.fq

-rw-rw-r--  1 1000 1000 8.2M May 27 13:41 sampleCtrl_3.R1.fq

-rw-rw-r--  1 1000 1000 8.6M May 27 13:42 sampleCtrl_3.R2.fq

-rw-rw-r--  1 1000 1000 8.2M May 27 13:40 sampleTreat_1.R1.fq

-rw-rw-r--  1 1000 1000 8.6M May 27 13:41 sampleTreat_1.R2.fq

-rw-rw-r--  1 1000 1000 8.2M May 27 13:40 sampleTreat_2.R1.fq

-rw-rw-r--  1 1000 1000 8.6M May 27 13:40 sampleTreat_2.R2.fq

-rw-rw-r--  1 1000 1000 8.2M May 27 13:40 sampleTreat_3.R1.fq

-rw-rw-r--  1 1000 1000 8.6M May 27 13:40 sampleTreat_3.R2.fq

-rw-rw-r--  1 1000 1000 1.2K May 28 12:28 samples.bowtie.PEsample

-rw-rw-r--  1 1000 1000 1.1K May 27 12:35 samples.bowtie.SEsample

-rw-rw-r--  1 1000 1000  868 May 28 12:49 samples.salmon.PEsample

-rw-rw-r--  1 1000 1000  840 May 28 14:13 samples.salmon.SEsample

drwxrwxr-x  1 1000 1000 4.0K Jun  8 06:59 scripts

-rw-r--r--  1 root root   18 May 29 06:14 test.fasta

 

ランにはサンプル情報を記載したファイルを使用する。

cat samples.bowtie.PEsample

# cat samples.bowtie.PEsample

###################################################################################

#       File "sample" - sample description file

#       Specify the names of the sample files and tag them as "treat" and "control".

###################################################################################

#       Specify the genome file, and specify the path where the

#       indexed file will be stored, and the prifex of the indexed genome.

#       Default is 'bowtie2_genome'.

GENOME SequenceChromosome.fasta bowtie2_genome/sequenceChr

#       Specify the fastq files

#       Specify the output name of the sam files.   

#       Followed by the tag/class/condition of the sample (treated or control)

#       List all the fastq files as bellow.

FASTQ sampleTreat_1.R1.fq sampleTreat_1.R2.fq sampleTreat_1.sam treat

FASTQ sampleTreat_2.R1.fq sampleTreat_2.R2.fq sampleTreat_2.sam treat

FASTQ sampleTreat_3.R1.fq sampleTreat_3.R2.fq sampleTreat_3.sam treat

FASTQ sampleCtrl_1.R1.fq sampleCtrl_1.R2.fq sampleCtrl_1.sam control

FASTQ sampleCtrl_2.R1.fq sampleCtrl_2.R2.fq sampleCtrl_2.sam control

FASTQ sampleCtrl_3.R1.fq sampleCtrl_3.R2.fq sampleCtrl_3.sam control

#

#       End of file "sample"

#

下の方のコメントアウトされていない6行(白)がfastq情報になる。ペアエンドのR1<space>R2<space>output.sam<space>case/control を記載する。(3群以上の多変量解析は行えない)。

実際のランでは、赤字の行のGENOME SequenceChromosome.fasta bowtie2_genome/sequenceChrを自分の使用しているgenomeのFASTAファイルとbowtie indexのディレクトリに修正する。indexはラン中に作成されるので、名前だけ指示しておけば良い(上ではbowtie2_genome/にsequenceChr~ができる)。

 

もう1つ必要なのはパラメータを記載したファイル。 exampleではparam.input.bowtie。テストはデフォルト設定でランできるが、実際のランでは白字のgtfとbedファイル名を修正する必要がある。また、GO enrichment解析はTERM2GENE.csvをもとに実行される(GO databaseはdata/TERM2NAME.csv)。GO enrichment解析を実行したいなら、事前に対応表を作って置いておく必要がある(*1)。

> cat param.input.bowtie

 cat param.input.bowtie 

####################################################################

#       File "param.input" - defination file

#       Define the paths and parameters to run the external packages.

####################################################################

#       Describe the Bowtie options below as:

BOWTIE -I 0

BOWTIE -X 500

BOWTIE -k 1

BOWTIE -p 1

# In case the package salmon is to be run, uncomment the options.

#       These are the default values.

#SALMONINDEX -k 29

#SALMONQUANT -l A

#SALMONQUANT -p 2

#SALMONQUANT --validateMappings TRUE

#       Define the Featurecounts options as bellow:

FEATURECOUNTS a oldAnnotationGFF.gtf

#       Define the Featurecounts output file name.

FEATURECOUNTS o FeatCount

#       Specify a count file name

COUNTFILE countFile.csv

#       geneBody coverage.r require a bed file. Specify the name of bed file as bellow:

geneBody_coverage r oldAnnotationGFF.bed

#       Specify if batch effect removal is required. FALSE if not required, else TRUE.

BATCH_EFFECT_REMOVE FALSE

#       For pathway analysis, define the log fold change cutoff as bellow:

PATHWAY cutoffPositive 2.0

PATHWAY cutoffNegative -2.0

#       For pathway analysis, define the organism in three alphabets as bellow.

#       ypy = Yersinia pseudotuberculosis

#       User need to change the keg abbreviation of their genome which can be

#       found in https://www.genome.jp/kegg/catalog/org_list.html. Here ypy is

#       the Yersinia pseudotuberculosis YPIII

PATHWAY Organism ypy

#       For Gene Ontology of the pathway analysis, define GO term and gene name file.

PATHWAY TERM2GENE data/TERM2GENE.csv

PATHWAY TERM2NAME data/TERM2NAME.csv

#       Specify the path to samtools

PATH SAMTOOLS /root/prokseq/depend/samtools/bin

#       Specify the root path. That means where the ProkSeq bundle is unpacked.

#       The location should have the following folders

#       1. depend - contains all the binaries of the external packages

#       2. scripts - contains all the modules required for ProkSeq, and pipeline-vx.x.sh

#       3. data - Contains Gene ontology files for pathway analysis.

PATH ROOT /root/prokseq

#       If the above environment (depend, scripts, data) is true, the following

#       line maye uncommented.

#PATH DEFAULT

#       Specify the path to geneBody_coverage

PATH geneBody_coverage /root/prokseq/depend/RSeQC-2.6.2/scripts/

#       Specify the path to FEATURECOUNTS

PATH FEATURECOUNTS /root/prokseq/depend/subread-1.4.6-p5-Linux-i386/bin/

#       Specify the path to fastqc

PATH FASTQC /root/prokseq/depend/FastQC

#       Specify the path to bowtie

PATH BOWTIE /root/prokseq/depend/bowtie2/bowtie2-2.3.5.1-linux-x86_64

#       Specify the path to pypy required for running afterqc

PATH PYPY /root/prokseq/depend/pypy2.7-v7.2.0-linux64/bin

#       Specify the path to readfasta

PATH READFASTA /root/prokseq/depend

#

#       End of file "param.input"

#

(py36) sh-5.0# 

注意;エラーが起こるので実際のランでも/root/proseq/を使うこと。

 

実行

python scripts/pipeline-v2.8.py -s samples.bowtie.PEsample -p param.input.bowtie -n 8
  • -s   provide the sample file
  • -p   provide the parameter file
  • -n   provide the number of processors

 

Checking required tools or packages:

 

       samtools : Found

         salmon : Not found

         fastqc : Found

         bowtie : Found

        afterqc : Found

      readFasta : Found

              R : Found

        GnBdCov : Found

        FeatCnt : Found

 

Checking R packages:

 

DESeq2 : TRUE

ggplot2 : TRUE

edgeR : TRUE

NOISeq : TRUE

limma : TRUE

clusterProfiler : TRUE

apeglm : TRUE

 

 

 

SALMON:

        salmon IS NOT FOUND!

 

 

     

Done with package checks. Seems all the required packages are available.

Do you want to continue? (Y/N) : Y

Yで実行。

 

PATHWAY ANALYSIS

 

Rscript --vanilla /root/prokseq/depend/clusterProf.r 2.0 -2.0 ypy data/TERM2GENE.csv data/TERM2NAME.csv

Output: 

 

Registered S3 method overwritten by 'enrichplot':

  method               from

  fortify.enrichResult DOSE

clusterProfiler v3.14.3  For help: https://guangchuangyu.github.io/software/clusterProfiler

 

If you use clusterProfiler in published research, please cite:

Guangchuang Yu, Li-Gen Wang, Yanyan Han, Qing-Yu He. clusterProfiler: an R package for comparing biological themes among gene clusters. OMICS: A Journal of Integrative Biology. 2012, 16(5):284-287.

[1] "ypy"

[1] "ypy"

[1] "Found organism"

`universe` is not in character and will be ignored...

 

 Success.

終了

 

出力

QC_preFilter/  #fastqc出力

f:id:kazumaxneo:20200628193238p:plain

QC_afterFilter/   #afterQC出力

f:id:kazumaxneo:20200628193630p:plain

sam/

f:id:kazumaxneo:20200628193436p:plain

bam/

f:id:kazumaxneo:20200628193352p:plain
wig/

f:id:kazumaxneo:20200628193452p:plain

AlignmentStat/

f:id:kazumaxneo:20200628193536p:plain

Results/

f:id:kazumaxneo:20200628193511p:plain

edgeR_results.txt

f:id:kazumaxneo:20200628214447p:plain

>DESeq2lfcShrink_results.txt

f:id:kazumaxneo:20200628214508p:plain

> DESeq2_results.txt

f:id:kazumaxneo:20200628214532p:plain

> afterNoiseq.txt

f:id:kazumaxneo:20200628214559p:plain

>KEGGpathway.txt

f:id:kazumaxneo:20200628214213p:plain

> KEGGenricher.txt

f:id:kazumaxneo:20200628214232p:plain

> GOpathways.txt

f:id:kazumaxneo:20200628214301p:plain

> GOenricher.txt

f:id:kazumaxneo:20200628214401p:plain

 

plots/

f:id:kazumaxneo:20200628193326p:plain

correlation_violin.pdf

f:id:kazumaxneo:20200628192851p:plain

correlation.pdf

f:id:kazumaxneo:20200628192911p:plain

DESeq-Volcano.pdf

f:id:kazumaxneo:20200628192926p:plain

edgeRvolcano.pdf

f:id:kazumaxneo:20200628192950p:plain

MA-plot.pdf

f:id:kazumaxneo:20200628193009p:plain

pca.pdf

f:id:kazumaxneo:20200628193024p:plain

pValue-histogram.pdf

f:id:kazumaxneo:20200628193040p:plain

RSeQC.geneBodyCoverage.curves.pdf

f:id:kazumaxneo:20200628193102p:plain


sampleCtrl_1.sortedBAM.bam_RSeQC.geneBodyCoverage.curves.pdf(サンプルごとに個別出力されるので合計6つ)

f:id:kazumaxneo:20200628193201p:plain

 

注意;ノートパソコなどで実行する場合は空き容量に注意して下さい。データによりますが、SAM出力するので、dockerの利用できる容量と空き容量を増やしておいて下さい。目安は100GB以上です。

引用

ProkSeq for complete analysis of RNA-seq data from prokaryotes
Authors: A K M Firoj Mahmud1, Soumyadeep Nandi#2, Maria Fällman

bioRxiv, Posted June 09, 2020.

 

*1

TERM2GENE.csvの中身

f:id:kazumaxneo:20200629210453p:plain

遺伝子名は他のアノテーションファイルと合致していること。 

 

TERM2NAME.csvはGO termを定義するファイル。

f:id:kazumaxneo:20200629210709p:plain