macでインフォマティクス

macでインフォマティクス

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

illuminaのショートリードシミュレータ Sandy(RNA seqにも対応)

 

 Sandyは、与えられたfastaファイルからシングルエンド/ペアエンドのリードを生成するシンプルなバイオインフォマティックツールである。多くの次世代シーケンシング分析は、実際には正確には満足されていない仮説モデルおよび原理に依存している。ポジティブコントロールを提供するシミュレーションデータは、これらの困難を克服するための完璧な方法である。ここでは、Sandyを紹介する。Sandyは、第2世代および第3世代シーケンスの合成リードをを生成する簡単で使いやすく、高速かつ完全なツールセットである。 Sandyは全ゲノムシーケンス、全エクソームシーケンス、RNAseqリードをシミュレートする。 Sandyはまた、ユーザーがデータを操作するためのいくつかの機能を提供するとともに、リードのゲノム上のポジション、遺伝子および転写産物の発現、シーケンシングエラー、「真の」情報(生成データに基づく)を含む体系化されたデータベースも提供する。Sandyの最も印象的な機能の1つは、シーケンシングエラーとともに、多型をsnv、indel、および構造変異(例えば、遺伝子重複、逆位、遺伝子融合)としてシミュレートする能力である - さらなる処理ステップを必要としない。そのため、Sandyは、ゲノミクスやトランスクリプトミクスにおけるさまざまなパイプラインの結果のベンチマーク、および新しい仮説の生成やシーケンスプロジェクトの最善の設計、おそらく時間とコストの最適化に役立つだろう。

 

HP

Welcome to Sandy simulator! ## | sandy

manual

https://galantelab.github.io/sandy/v0.22/main.html#docker-usage

 

インストール

ubuntu16.04でテストした。

依存

#step1 Debian/Ubuntu
apt update && apt install perl zlib1g-dev gcc make

#step1 Centos/Fedora
yum install perl zlib gcc make

#step2 install cpanminus package 
cpan -i App::cpanminus

本体 Github

cpanm App::Sandy

 

またはオーサーらが準備してくれているdockerイメージを使う

docker pull galantelab/sandy

 >sandy

# sandy

Usage:

     sandy [options]

     sandy help <command>

     sandy <command> ...

 

     Options:

      -h, --help               brief help message

      -u, --man                full documentation

 

     Help commands:

      help                     show application or command-specific help

      man                      show application or command-specific documentation

 

     Misc commands:

      version                  print the current version

      citation                 export citation in BibTeX format

 

     Database commands:

      quality                  manage quality profile database

      expression               manage expression-matrix database

      variation                manage structural variation database

 

     Main commands:

      genome                   simulate genome sequencing

      transcriptome            simulate transcriptome sequencing

 

sandy genome --help

# sandy genome --help

Usage:

     sandy genome [options] <fasta-file>

 

     Arguments:

      a fasta-file

 

     Options:

      -h, --help                         brief help message

      -u, --man                          full documentation

      -v, --verbose                      print log messages

      -p, --prefix                       prefix output [default:"out"]

      -o, --output-dir                   output directory [default:"."]

      -O, --output-format                bam, sam, fastq.gz, fastq [default:"fastq.gz"]

      -1, --join-paired-ends             merge R1 and R2 outputs in one file

      -x, --compression-level            speed compression: "1" - compress faster,

                                         "9" - compress better [default:"6"; Integer]

      -i, --append-id                    append to the defined template id [Format]

      -I, --id                           overlap the default template id [Format]

      -j, --jobs                         number of jobs [default:"1"; Integer]

      -s, --seed                         set the seed of the base generator

                                         [default:"time()"; Integer]

      -c, --coverage                     genome coverage [default:"8", Number]

      -t, --sequencing-type              single-end or paired-end reads

                                         [default:"paired-end"]

      -q, --quality-profile              sequencing system profiles from quality

                                         database [default:"poisson"]

      -e, --sequencing-error             sequencing error rate for poisson

                                         [default:"0.001"; Number]

      -m, --read-mean                    read mean size for poisson

                                         [default:"100"; Integer]

      -d, --read-stdd                    read standard deviation size for poisson

                                         [default:"0"; Integer]

      -M, --fragment-mean                the fragment mean size for paired-end reads

                                         [default:"300"; Integer]

      -D, --fragment-stdd                the fragment standard deviation size for

                                         paired-end reads [default:"50"; Integer]

      -a, --genomic-variation            a list of genomic variation entries from

                                         variation database. This option may be passed

                                         multiple times [default:"none"]

      -A, --genomic-variation-regex      a list of perl-like regex to match genomic

                                         variation entries in variation database.

                                         This option may be passed multiple times

                                         [default:"none"]

 

> sandy transcriptome -h

# sandy transcriptome -h

Usage:

     sandy transcriptome [options] <fasta-file>

 

     Arguments:

      a fasta-file

 

     Options:

      -h, --help                     brief help message

      -u, --man                      full documentation

      -v, --verbose                  print log messages

      -p, --prefix                   prefix output [default:"out"]  

      -o, --output-dir               output directory [default:"."]

      -O, --output-format            bam, sam, fastq.gz, fastq [default:"fastq.gz"]

      -1, --join-paired-ends         merge R1 and R2 outputs in one file

      -x, --compression-level        speed compression: "1" - compress faster,

                                     "9" - compress better [default:"6"; Integer]

      -i, --append-id                append to the defined template id [Format]

      -I, --id                       overlap the default template id [Format]

      -j, --jobs                     number of jobs [default:"1"; Integer]

      -s, --seed                     set the seed of the base generator

                                     [default:"time()"; Integer]

      -n, --number-of-reads          set the number of reads

                                     [default:"1000000", Integer]

      -t, --sequencing-type          single-end or paired-end reads

                                     [default:"paired-end"]

      -q, --quality-profile          sequencing system profiles from quality

                                     database [default:"poisson"]

      -e, --sequencing-error         sequencing error rate for poisson

                                     [default:"0.001"; Number]

      -m, --read-mean                read mean size for poisson

                                     [default:"100"; Integer]

      -d, --read-stdd                read standard deviation size for poisson

                                     [default:"0"; Integer]

      -M, --fragment-mean            the fragment mean size for paired-end reads

                                     [default:"300"; Integer]

      -D, --fragment-stdd            the fragment standard deviation size for

                                     paired-end reads [default:"50"; Integer]

      -f, --expression-matrix        an expression-matrix entry from database

 

root@84830ec770b4:/data# 

> sandy expression #RNAseq向けデータベース組み込み済み発現パターン

# sandy expression

.----------------------------------------------------------------------------------.

| expression-matrix                   | source             | provider | date       |

+-------------------------------------+--------------------+----------+------------+

| adipose_subcutaneous                | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| adipose_visceral                    | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| adrenal_gland                       | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| artery_aorta                        | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| artery_coronary                     | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| artery_tibial                       | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| bladder                             | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_amygdala                      | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_anterior_cingulate_cortex     | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_caudate                       | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_cerebellar_hemisphere         | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_cerebellum                    | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_cortex                        | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_frontal_cortex                | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_hippocampus                   | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_hypothalamus                  | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_nucleus_accumbens             | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_putamen                       | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_spinal_cord                   | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| brain_substantia_nigra              | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| breast_mammary_tissue               | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| cells_ebv_transformed_lymphocytes   | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| cells_leukemia_cell_line            | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| cells_transformed_fibroblasts       | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| cervix_ectocervix                   | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| cervix_endocervix                   | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| colon_sigmoid                       | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| colon_transverse                    | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| esophagus_gastroesophageal_junction | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| esophagus_mucosa                    | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| esophagus_muscularis                | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| fallopian_tube                      | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| heart_atrial_appendage              | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| heart_left_ventricle                | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| kidney_cortex                       | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| liver                               | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| lung                                | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| minor_salivary_gland                | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| muscle_skeletal                     | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| nerve_tibial                        | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| ovary                               | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| pancreas                            | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| pituitary                           | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| prostate                            | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| skin_not_sun_exposed                | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| skin_sun_exposed                    | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| small_intestine_terminal_ileum      | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| spleen                              | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| stomach                             | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| testis                              | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| thyroid                             | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| uterus                              | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| vagina                              | Xena GTEx Kallisto | vendor   | 2018-08-08 |

| whole_blood                         | Xena GTEx Kallisto | vendor   | 2018-08-08 |

'-------------------------------------+--------------------+----------+------------'

詳細なマニュアル

> sandy genome --man

 

dockerでヘルプを見るなら

> docker run galantelab/sandy --help

$ docker run galantelab/sandy --help

Usage:

     sandy [options]

     sandy help <command>

     sandy <command> ...

 

     Options:

      -h, --help               brief help message

      -u, --man                full documentation

 

     Help commands:

      help                     show application or command-specific help

      man                      show application or command-specific documentation

 

     Misc commands:

      version                  print the current version

      citation                 export citation in BibTeX format

 

     Database commands:

      quality                  manage quality profile database

      expression               manage expression-matrix database

      variation                manage structural variation database

 

     Main commands:

      genome                   simulate genome sequencing

      transcriptome            simulate transcriptome sequencing

 

 


実行方法

1、genome

ペアエンドfastqをゲノムの20x発生させる。

sandy genome --verbose --sequencing-type=paired-end --coverage=20 input_genome.fa -o outdir

#dockerなら
docker run -v $PWD:/data/ galantelab/sandy genome --verbose --sequencing-type=paired-end --coverage=20 /data/input_genome.fa -o /data/outdir

default settings(manualより)

  • The strand is randomly chosen;
  • The number of reads is calculated by the coverage;
  • The chromossomes are raffled following a weighted raffle with the sequence length as the bias;

-s <integer>で乱数発生を固定し(defaltはtimeを使う)、再現性ある結果をreproduceできる。リードのidentifier(fastqの1、3行目)は"--id="で指定できる。フォーマットはCのprintf関数に似ている。詳細はdocumentation参照(リンク)。quaity profileのdefaultはhiseq。

 出力

f:id:kazumaxneo:20190304220750p:plain

 

2、transcriptome

sandy transcriptome \
--quality-profile=hiseq_101 \
--sequencing-type=paired-end \
--number-of-reads=1000000
--fragment-mean=350 \
--fragment-stdd=100 \
--prefix=pancreas_sim \
--output-dir=simulation_dir \
--id="%i.%U read=%c:%t-%n mate=%c:%T-%N length=%r" \
--verbose \
--seed=123 \
--jobs=20 \
input_cDNAs.fa

default settings(manualより)

  • Choose the Minus strand;
  • The number of reads is directly passed;
  • The genes/transcripts are raffled following the expression matrix;

上の条件で実行して得たfastqのstatistics

file                 format  type   num_seqs      sum_len  min_len  avg_len  max_len

out_R1_001.fastq.gz  FASTQ   DNA   1,000,000  100,000,000      100      100      100

out_R2_001.fastq.gz  FASTQ   DNA   1,000,000  100,000,000      100      100      100

$ head -n 4 out_R1_001.fastq

$ head -n 4 out_R1_001.fastq 

@SR.1:syn:sll0006  aspC, aat; putative aminotransferase; K00837 1

agggctgtttctgctgctttttgtacgactcgaaaaattccataatccaagttagttttcagcgttctcaaaccttgaataatttcttgatttcccacga

+

@IAFFAHFACIDCAEH@FFIIFBEAHEFGEGEBIFGDBCBFABBHDEEHI@@EIFHF?@FFIDCAH78==658><8<6=9=99.,11-2+/+#&(%##%$

 

公式Documentでは以下の実行例がある。

sandy transcriptome \
--quality-profile=hiseq_101 \
--sequencing-type=paired-end \
--fragment-mean=350 \
--fragment-stdd=100 \
--prefix=pancreas_sim \
--output-dir=sim_dir \
--id="%i.%U read=%c:%t-%n mate=%c:%T-%N length=%r" \
--verbose \
--seed=123 \
--jobs=30 \
--no-gzip \
--expression-matrix=brain_cortex \
gencode_pc_v26.fa.gz

expression-matrixはリファレンスと紐づいており、ヘッダーが違うと機能しないので注意。

 

自分の発現データとリファレンスを組み込み、"-expression-matrix"フラグで使えるようにする。

sandy expression add my_expression.txt
sandy expression -f my_expression.txt my_fasta.fa

 

他にもquality profileを追加するコマンドなどがある。

引用

galantelab/sandy: A straightforward and complete next-generation sequencing read simulator

Thiago Miller. (2018, May 6).
Zenodo. http://doi.org/10.5281/zenodo.1241587