macでインフォマティクス

macでインフォマティクス

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

eukaryotesのアノテーションパイプライン BRAKER2

2020 8/14 補足追記

 

 遺伝子予測の完全自動化は、次世代シーケンシングの出現以来、バイオインフォマティクスの重要な課題となっている。真核生物ゲノムアノテーションパイプラインBRAKER1では、自己学習型のGeneMark ETとAUGUSTUSを組み合わせて、トランスクリプトームデータをサポートした遺伝子座標(genes coordinates)を生成していた。ここでは、GeneMark EP+とAUGUSTUSによる外部からのタンパク質配列のゲノムアラインメントをサポートしたパイプラインBRAKER2を紹介する。このパイプラインの開発において課題となったのは、相同性はあるが進化的に遠いタンパク質から、タンパク質をコードするエクソン境界の位置を信頼できるヒントを生成することであった。同等の条件で、BRAKER2の遺伝子予測精度は、他のゲノムアノテーションパイプラインであるMAKER2よりも高いことが示された。また、大量のトランスクリプトデータに支えられたBRAKER1と比較して、タンパク質データベースの参照種との進化距離が小さい場合には、BRAKER2の方がより良い遺伝子予測精度が得られることが示された。また、その結果、BRAKER2は真核生物ゲノムの構造アノテーションを高速かつ正確に行うことができることが示された。

 

下に引用したMethods in Molecular Biologyの論文ではbrakerを使う流れがコマンドも含めて詳細に記述されています。確認して下さい。

 

インストール

dockerhubで公開されているイメージを使ってテストした。

Github

#condaで導入可能(試してません)
#bioconda(link) (v1は別link)
conda create -n braker2 -y
conda activate braker2
conda install -c bioconda braker2 -y

 

ここでは公開されているdockerイメージを使う。nanozooさんが公開しているイメージを選んだ。

docker pull nanozoo/braker2
docker run -it nanozoo/braker2

#cdbfastaがないので入れる(condaでも入る)。
apt update && apt-get install cdbfasta

#cpanmを使ってperl libraryを導入,makeも入れる
apt update && apt install make cpanminus
cpanm YAML Hash::Merge Logger::Simple Parallel::ForkManager MCE::Mutex Thread::Queue threads

#GeneMarkはライセンスの関係でHPより直接ダウンロードする必要がある (link)。プログラム本体とキーをダウンロードして、
cd GeneMark_dir
#Install the key
zcat gm_key_64.gz ~/.gm_key
#test run
bash check_install.bash
#パスを通す
export GENEMARK_PATH=/your_path_to_GeneMark-ET/

#diamond
conda install -c bioconda -y diamond

#以下のperlライブラリも入れた
cpanm Exception/Class.pm Object::InsideOut Clone::Choose Class::Data::Inheritable Devel::StackTrace Moo::Role Module::Runtime List::MoreUtils

braker.pl -v  #version 2.1.4

# braker.pl -v

Option v is ambiguous (verbosity, version)

Use of uninitialized value $epath in concatenation (.) or string at /opt/conda/bin/braker.pl line 1986.

#**********************************************************************************

#                               BRAKER CONFIGURATION                               

#**********************************************************************************

# BRAKER CALL: /opt/conda/bin/braker.pl -v

# Thu Aug 13 06:07:03 2020: braker.pl version 2.1.4

# Thu Aug 13 06:07:03 2020: Configuring of BRAKER for using external tools...

# Thu Aug 13 06:07:03 2020: Found environment variable $AUGUSTUS_CONFIG_PATH. Setting $AUGUSTUS_CONFIG_PATH to /opt/conda/config/

# Thu Aug 13 06:07:03 2020: Found environment variable $AUGUSTUS_BIN_PATH. Setting $AUGUSTUS_BIN_PATH to /opt/conda/bin/

# Thu Aug 13 06:07:03 2020: Found environment variable $AUGUSTUS_SCRIPTS_PATH. Setting $AUGUSTUS_SCRIPTS_PATH to /opt/conda/bin/

# Thu Aug 13 06:07:03 2020: Did not find environment variable $GENEMARK_PATH  (either variable does not exist, or the path given in variable does not exist). Will try to set this variable in a different way, later.

# Thu Aug 13 06:07:03 2020: Trying to guess $GENEMARK_PATH from location of gmes_petap.pl executable that is available in your $PATH

#*********

# WARNING: Guessing the location of $GENEMARK_PATH failed.  is not a directory!

#*********

# Thu Aug 13 06:07:03 2020: ERROR: in file /opt/conda/bin/braker.pl at line 2015

$GENEMARK_PATH not set!

There are 3 alternative ways to set this variable for

braker.pl:

   a) provide command-line argument --GENEMARK_PATH=/your/path

   b) use an existing environment variable $GENEMARK_PATH

      for setting the environment variable, run

           export GENEMARK_PATH=/your/path

      in your shell. You may append this to your .bashrc or 

      .profile file in order to make the variable available to

      all your bash sessions.

   c) braker.pl can try guessing the location of 

      $GENEMARK_PATH from the location of a gmes_petap.pl

      executable that is available in your $PATH variable.

      If you try to rely on this option, you can check by

      typing

           which gmes_petap.pl

      in your shell, whether there is a bamtools executable in

      your $PATH

braker.pl

# braker.pl 

 

braker.pl     Pipeline for predicting genes with GeneMark-ET and AUGUSTUS with

            RNA-Seq

 

SYNOPSIS

 

braker.pl [OPTIONS] --genome=genome.fa --bam=rnaseq.bam

 

INPUT FILE OPTIONS

 

--genome=genome.fa                  fasta file with DNA sequences

--bam=rnaseq.bam                    bam file with spliced alignments from

                                    RNA-Seq

--hints=hints.gff                   Alternatively to calling braker.pl with a

                                    bam file, it is possible to call it with a

                                    file that contains introns extracted from

                                    RNA-Seq (or other data) in gff format.

                                    This flag also allows the usage of hints

                                    from additional extrinsic sources for gene

                                    prediction with AUGUSTUS. To consider such

                                    additional extrinsic information, you need

                                    to use the flag --extrinsicCfgFiles to

                                    specify parameters for all sources in the

                                    hints file (including the source "E" for

                                    intron hints from RNA-Seq).

--prot_seq=prot.fa                  A protein sequence file in multiple fasta

                                    format. This file will be used to generate

                                    protein hints for AUGUSTUS by running one

                                    of the three alignment tools Exonerate

                                    (--prg=exonerate), Spaln (--prg=spaln) or

                                    GenomeThreader (--prg=gth). Default is

                                    GenomeThreader if the tool is not

                                    specified.  Currently, hints from protein

                                    sequences are only used in the prediction

                                    step with AUGUSTUS.

--prot_aln=prot.aln                 Alignment file generated from aligning

                                    protein sequences against the genome with

                                    either Exonerate (--prg=exonerate), or

                                    Spaln (--prg=spaln), or GenomeThreader

                                    (--prg=gth).

                                    To prepare alignment file, run Spaln2 with

                                    the following command:

                                    spaln -O0 ... > spalnfile

                                    To prepare alignment file, run Exonerate

                                    with the following command:

                                    exonerate --model protein2genome \

                                        --showtargetgff T ... > exfile

                                    To prepare alignment file, run

                                    GenomeThreader with the following command:

                                    gth -genomic genome.fa  -protein \

                                        protein.fa -gff3out \

                                        -skipalignmentout ... -o gthfile

                                    A valid option prg=... must be specified

                                    in combination with --prot_aln. Generating

                                    tool will not be guessed.

                                    Currently, hints from protein alignment

                                    files are only used in the prediction step

                                    with AUGUSTUS.

--AUGUSTUS_ab_initio                output ab initio predictions by AUGUSTUS

                                    in addition to predictions with hints by

                                    AUGUSTUS

 

FREQUENTLY USED OPTIONS

 

--species=sname                     Species name. Existing species will not be

                                    overwritten. Uses Sp_1 etc., if no species

                                    is assigned

--softmasking                       Softmasking option for soft masked genome

                                    files. (Disabled by default.)

--esmode                            Run GeneMark-ES (genome sequence only) and 

                                    train AUGUSTUS on long genes predicted by 

                                    GeneMark-ES. Final predictions are ab initio

--epmode                            Run GeneMark-EP with intron hints provided

                                    from protein data. This mode is not 

                                    comptabile with using the aligners

                                    GenomeThreader, Exonerate and Spaln within

                                    braker.pl because etpmode and epmode require

                                    a large database of proteins and such

                                    mapping should be done outside of braker.pl 

                                    e.g. on a cluster.

--etpmode                           Run GeneMark-ETP with hints provided from

                                    proteins and RNA-Seq data. This mode is not

                                    compatible with using the aligners

                                    GenomeThreader, Exonerate and Spaln within

                                    braker.pl because etpmode and epmode require

                                    a large database of proteins and such

                                    mapping should be done outside of braker.pl 

                                    e.g. on a cluster.

--gff3                              Output in GFF3 format (default is gtf

                                    format)

--cores                             Specifies the maximum number of cores that

                                    can be used during computation. Be aware:

                                    optimize_augustus.pl will use max. 8

                                    cores; augustus will use max. nContigs in

                                    --genome=file cores.

--workingdir=/path/to/wd/           Set path to working directory. In the

                                    working directory results and temporary

                                    files are stored

--nice                              Execute all system calls within braker.pl

                                    and its submodules with bash "nice"

                                    (default nice value)

 

--alternatives-from-evidence=true   Output alternative transcripts based on

                                    explicit evidence from hints (default is

                                    true).

--crf                               Execute CRF training for AUGUSTUS;

                                    resulting parameters are only kept for

                                    final predictions if they show higher

                                    accuracy than HMM parameters.

--keepCrf                           keep CRF parameters even if they are not

                                    better than HMM parameters

--UTR=on                            create UTR training examples from RNA-Seq

                                    coverage data; requires options

                                    --bam=rnaseq.bam and --softmasking.

                                    Alternatively, if UTR parameters already

                                    exist, training step will be skipped and

                                    those pre-existing parameters are used.

--prg=gth|exonerate|spaln           Alignment tool gth (GenomeThreader),

                                    exonerate (Exonerate) or Spaln2

                                    (spaln) that will be used to generate

                                    protein alignments that will be the

                                    basis for hints generation for gene

                                    prediction with AUGUSTUS (if specified

                                    in combination with --prot_seq) or that

                                    was used to externally generate an

                                    alignment file with the commands listed in

                                    description of --prot_aln (if used in

                                    combination with --prot_aln).

--gth2traingenes                    Generate training gene structures for

                                    AUGUSTUS from GenomeThreader alignments.

                                    (These genes can either be used for

                                    training AUGUSTUS alone with

                                    --trainFromGth; or in addition to

                                    GeneMark-ET training genes if also a

                                    bam-file is supplied.)

--trainFromGth                      No GeneMark-Training, train AUGUSTUS from

                                    GenomeThreader alignments

--makehub                           Create track data hub with make_hub.py 

                                    for visualizing BRAKER results with the

                                    UCSC GenomeBrowser

--email                             E-mail address for creating track data hub

--version                           Print version number of braker.pl

--help                              Print this help message

 

CONFIGURATION OPTIONS (TOOLS CALLED BY BRAKER)

 

--AUGUSTUS_CONFIG_PATH=/path/       Set path to config directory of AUGUSTUS

                                    (if not specified as environment

                                    variable). BRAKER1 will assume that the

                                    directories ../bin and ../scripts of

                                    AUGUSTUS are located relative to the

                                    AUGUSTUS_CONFIG_PATH. If this is not the

                                    case, please specify AUGUSTUS_BIN_PATH

                                    (and AUGUSTUS_SCRIPTS_PATH if required).

                                    The braker.pl commandline argument

                                    --AUGUSTUS_CONFIG_PATH has higher priority

                                    than the environment variable with the

                                    same name.

--AUGUSTUS_BIN_PATH=/path/          Set path to the AUGUSTUS directory that

                                    contains binaries, i.e. augustus and

                                    etraining. This variable must only be set

                                    if AUGUSTUS_CONFIG_PATH does not have

                                    ../bin and ../scripts of AUGUSTUS relative

                                     to its location i.e. for global AUGUSTUS

                                    installations. BRAKER1 will assume that

                                    the directory ../scripts of AUGUSTUS is

                                    located relative to the AUGUSTUS_BIN_PATH.

                                    If this is not the case, please specify

                                    --AUGUSTUS_SCRIPTS_PATH.

--AUGUSTUS_SCRIPTS_PATH=/path/      Set path to AUGUSTUS directory that

                                    contains scripts, i.e. splitMfasta.pl.

                                    This variable must only be set if

                                    AUGUSTUS_CONFIG_PATH or AUGUSTUS_BIN_PATH

                                    do not contains the ../scripts directory

                                    of AUGUSTUS relative to their location,

                                    i.e. for special cases of a global

                                    AUGUSTUS installation.

--BAMTOOLS_PATH=/path/to/           Set path to bamtools (if not specified as

                                    environment BAMTOOLS_PATH variable). Has

                                    higher priority than the environment

                                    variable.

--GENEMARK_PATH=/path/to/           Set path to GeneMark-ET (if not specified

                                    as environment GENEMARK_PATH variable).

                                    Has higher priority than environment

                                    variable.

--SAMTOOLS_PATH=/path/to/           Optionally set path to samtools (if not

                                    specified as environment SAMTOOLS_PATH

                                    variable) to fix BAM files automatically,

                                    if necessary. Has higher priority than

                                    environment variable.

--ALIGNMENT_TOOL_PATH=/path/to/tool Set path to alignment tool

                                    (GenomeThreader, Spaln, or Exonerate) if

                                    not specified as environment

                                    ALIGNMENT_TOOL_PATH variable. Has higher

                                    priority than environment variable.

--DIAMOND_PATH=/path/to/diamond     Set path to diamond, this is an alternative

                                    to NCIB blast; you only need to specify one 

                                    out of DIAMOND_PATH or BLAST_PATH, not both.

                                    DIAMOND is a lot faster that BLAST and yields 

                                    highly similar results for BRAKER.

--BLAST_PATH=/path/to/blastall      Set path to NCBI blastall and formatdb

                                    executables if not specified as

                                    environment variable. Has higher priority

                                    than environment variable.

--PYTHON3_PATH=/path/to             Set path to python3 executable (if not 

                                    specified as envirnonment variable and if

                                    executable is not in your $PATH).

--MAKEHUB_PATH=/path/to             Set path to make_hub.py (if option --makehub

                                    is used).

--CDBTOOLS_PATH=/path/to            cdbfasta/cdbyank are required for running

                                    fix_in_frame_stop_codon_genes.py. Usage of

                                    that script can be skipped with option 

                                    '--skip_fixing_broken_genes'.

 

 

EXPERT OPTIONS

 

--augustus_args="--some_arg=bla"    One or several command line arguments to

                                    be passed to AUGUSTUS, if several

                                    arguments are given, separated by

                                    whitespace, i.e.

                                    "--first_arg=sth --second_arg=sth".

--overwrite                         Overwrite existing files (except for

                                    species parameter files)

--skipGeneMark-ES                   Skip GeneMark-ES and use provided

                                    GeneMark-ES output (e.g. provided with 

                                    --geneMarkGtf=genemark.gtf)

--skipGeneMark-ET                   Skip GeneMark-ET and use provided

                                    GeneMark-ET output (e.g. provided with 

                                    --geneMarkGtf=genemark.gtf)

--skipGeneMark-EP                   Skip GeneMark-EP and use provided

                                    GeneMark-EP output (e.g. provided with

                                    --geneMarkGtf=genemark.gtf)

--skipGeneMark-ETP                  Skip GeneMark-ETP and use provided

                                    GeneMark-ETP output (e.g. provided with

                                    --geneMarkGtf=genemark.gtf)

--geneMarkGtf=file.gtf              If skipGeneMark-ET is used, braker will by

                                    default look in the working directory in

                                    folder GeneMarkET for an already existing

                                    gtf file. Instead, you may provide such a

                                    file from another location. If geneMarkGtf

                                    option is set, skipGeneMark-ES/ET/EP/ETP is

                                    automatically also set.

--rounds                            The number of optimization rounds used in

                                    optimize_augustus.pl (default 5)

--skipAllTraining                   Skip GeneMark-EX (training and

                                    prediction), skip AUGUSTUS training, only

                                    runs AUGUSTUS with pre-trained and already

                                    existing parameters (not recommended).

                                    Hints from input are still generated.

                                    This option automatically sets

                                    --useexisting to true.

--useexisting                       Use the present config and parameter files

                                    if they exist for 'species'

--filterOutShort                    It may happen that a "good" training gene,

                                    i.e. one that has intron support from

                                    RNA-Seq in all introns predicted by

                                    GeneMark, is in fact too short. This flag

                                    will discard such genes that have

                                    supported introns and a neighboring

                                    RNA-Seq supported intron upstream of the

                                    start codon within the range of the

                                    maximum CDS size of that gene and with a

                                    multiplicity that is at least as high as

                                    20% of the average intron multiplicity of

                                    that gene.

--skipOptimize                      Skip optimize parameter step (not

                                    recommended).

--skipGetAnnoFromFasta              Skip calling the python3 script

                                    getAnnoFastaFromJoingenes.py from the

                                    AUGUSTUS tool suite. This script requires

                                    python3, biopython and re (regular 

                                    expressions) to be installed. It produces

                                    coding sequence and protein FASTA files 

                                    from AUGUSTUS gene predictions and provides

                                    information about genes with in-frame stop 

                                    codons. If you enable this flag, these files 

                                    will not be produced and python3 and 

                                    the required modules will not be necessary

                                    for running braker.pl.

--skip_fixing_broken_genes          If you do not have python3, you can choose

                                    to skip the fixing of stop codon including

                                    genes (not recommended).

--fungus                            GeneMark-ET option: run algorithm with

                                    branch point model (most useful for fungal

                                    genomes)

--rnaseq2utr_args=params            Expert option: pass alternative parameters

                                    to rnaseq2utr as string, default parameters:

                                    -r 76 -v 100 -n 15 -i 0.7 -m 0.3 -w 70 

                                    -c 100 -p 0.5 

--eval=reference.gtf                Reference set to evaluate predictions

                                    against (using the eval package)

--AUGUSTUS_hints_preds=s            File with AUGUSTUS hints predictions; will

                                    use this file as basis for UTR training;

                                    only UTR training and prediction is

                                    performed if this option is given.

--flanking_DNA=n                    Size of flanking region, must only be

                                    specified if --AUGUSTUS_hints_preds is given

                                    (for UTR training in a separate braker.pl 

                                    run that builds on top of an existing run)

--verbosity=n                       0 -> run braker.pl quiet (no log)

                                    1 -> only log warnings

                                    2 -> also log configuration

                                    3 -> log all major steps

                                    4 -> very verbose, log also small steps

--downsampling_lambda=d             The distribution of introns in training

                                    gene structures generated by GeneMark-EX

                                    has a huge weight on single-exon and

                                    few-exon genes. Specifying the lambda

                                    parameter of a poisson distribution will

                                    make braker call a script for downsampling

                                    of training gene structures according to

                                    their number of introns distribution, i.e.

                                    genes with none or few exons will be

                                    downsampled, genes with many exons will be

                                    kept. Default value is 2. 

                                    If you want to avoid downsampling, you have 

                                    to specify 0. 

--checkSoftware                     Only check whether all required software

                                    is installed, no execution of BRAKER

--nocleanup                         Skip deletion of all files that are typically not 

                                    used in an annotation project after 

                                    running braker.pl. (For tracking any 

                                    problems with a braker.pl run, you 

                                    might want to keep these files, therefore

                                    nocleanup can be activated.)

 

 

DEVELOPMENT OPTIONS (PROBABLY STILL DYSFUNCTIONAL)

 

--splice_sites=patterns             list of splice site patterns for UTR

                                    prediction; default: GTAG, extend like this:

                                    --splice_sites=GTAG,ATAC,...

                                    this option only affects UTR training

                                    example generation, not gene prediction

                                    by AUGUSTUS

--extrinsicCfgFiles=file1,file2,... Depending on the mode in which braker.pl

                                    is executed, it may require one ore several

                                    extrinsicCfgFiles. Don't use this option

                                    unless you know what you are doing!

--stranded=+,-,+,-,...              If UTRs are trained, i.e.~strand-specific

                                    bam-files are supplied and coverage 

                                    information is extracted for gene prediction, 

                                    create stranded ep hints. The order of 

                                    strand specifications must correspond to the

                                    order of bam files. Possible values are

                                    +, -, .

                                    If stranded data is provided, ONLY coverage

                                    data from the stranded data is used to 

                                    generate UTR examples! Coverage data from 

                                    unstranded data is used in the prediction

                                    step, only.

                                    The stranded label is applied to coverage

                                    data, only. Intron hints are generated

                                    from all libraries treated as "unstranded"

                                    (because splice site filtering eliminates

                                    intron hints from the wrong strand, anyway).

--optCfgFile=ppx.cfg                Optional custom config file for AUGUSTUS

                                    for running PPX (currently not

                                    implemented)

--grass                             Switch this flag on if you are using braker.pl 

                                    for predicting genes in grasses with 

                                    GeneMark-ES/ET. The flag will enable

                                    GeneMark-ES/ET to handle GC-heterogenicity

                                    within genes more properly.

                                    NOTHING IMPLEMENTED FOR GRASS YET!

--transmasked_fasta=file.fa         Transmasked genome FASTA file for GeneMark

                                    (to be used instead of the regular genome

                                    FASTA file).  

--min_contig=INT                    Minimal contig length for GeneMark, could

                                    for example be set to 10000 if transmasked_fasta

                                    option is used because transmasking might

                                    introduce many very short contigs.

--translation_table=INT             Change translation table from non-standard

                                    to something else. 

                                    DOES NOT WORK YET BECAUSE BRAKER DOESNT

                                    SWITCH TRANSLATION TABLE FOR GENEMARK-EX, YET!

 

 

DESCRIPTION

 

Example:

 

 

braker.pl [OPTIONS] --genome=genome.fa --species=speciesname \

    --bam=accepted_hits.bam

braker.pl [OPTIONS] --genome=genome.fa --species=speciesname \

    --hints=rnaseq.gff

 

To run with protein data from remote species and GeneMark-EP:

 

braker.pl [OPTIONS] --genome=genome.fa --hints=proteinintrons.gff --epmode=1

 

To run with protein data from a very closely related species:

 

braker.pl [OPTIONS] --genome=genome.fa --prot_seq=proteins.fa --prg=gth \

    --gth2traingenes --trainFromGth

 

 

 

テストラン

git clone https://github.com/Gaius-Augustus/BRAKER.git
cd BRAKER/example/
wget http://topaz.gatech.edu/GeneMark/Braker/RNAseq.bam

f:id:kazumaxneo:20200813165357p:plain

brakerを実行する。

#Testing BRAKER with proteins of close homoogy and RNA-Seq data
braker.pl --genome genome.fa --bam RNAseq.bam --softmasking --cores 8

#Testing BRAKER with pre-trained parameters
braker.pl --genome=genome.fa --bam RNAseq.bam --species=arabidopsis \
--skipAllTraining --softmasking --cores 8

#Testing BRAKER with genome sequence
braker.pl --genome=genome.fa --esmode --softmasking --cores 8

#Testing BRAKER with proteins of close homology
braker.pl --genome genome.fa --prot_seq proteins.fa --prg gth \
--trainFromGth --softmasking --cores 8
  • --genome     fasta file with DNA sequences
  • --bam    bam file with spliced alignments from RNA-Seq
  • --softmasking    Softmasking option for soft masked genome files. (Disabled by default.)
  • --cores    Specifies the maximum number of cores that can be used during computation. Be aware: optimize_augustus.pl will use max. 8 cores; augustus will use max. nContigs in --genome=file cores.

ls /opt/conda/config//species/

f:id:kazumaxneo:20200813184848p:plain

ls -lh /opt/conda/config//species/arabidopsis/

f:id:kazumaxneo:20200813184954p:plain

 

 

 

実行方法

genomeとbamファイルを指定する。

braker.pl --species=yourSpecies --genome=genome.fasta \
--bam=file1.bam,file2.bam --cores 8

#リピートがソフトマスクされているゲノム(リピート配列が小文字に置き換えられている)
braker.pl --species=yourSpecies --genome=genome.fasta \
--bam=file.bam --softmasking --UTR=on --cores 8
  • --softmasking    Softmasking option for soft masked genome files. (Disabled by default.)
  • --UTR=on    create UTR training examples from RNA-Seq coverage data; requires options --bam=rnaseq.bam and --softmasking.  Alternatively, if UTR parameters already exist, training step will be skipped and  those pre-existing parameters are used.

UTR parametersをつけるときはStranded RNA-Seqかどうかも重要になる。詳しくはGithub参照。

 

すでにbamや他のソースからgtfを得ている場合、hintフラグを立ててgtfファイルを指定する。

braker.pl --species=yourSpecies --genome=genome.fasta \
--hints=hints1.gff,hints2.gff --cores 8

 

RNA seqなどの信頼性の高いエビデンスが利用できない場合、近縁なモデル生物のprotein情報を使用する。

braker.pl --genome=genome.fa --prot_seq=proteins.fa --softmasking --cores 8

 

引用

BRAKER2: Automatic Eukaryotic Genome Annotation with GeneMark-EP+ and AUGUSTUS Supported by a Protein Database

Tomas Bruna, Katharina Hoff, Mario Stanke, Alexandre Lomsadze, Mark Borodovsky

bioRxiv, Posted August 11, 2020

 

Whole-Genome Annotation with BRAKER
Katharina Hoff, Alexandre Lomsadze, Mark Borodovsky, Mario Stanke

Methods Mol Biol. 2019;1962:65-95

 

 

補足

GeneMarkは最新版とキーを両方ダウンロードし、キーと合わせて使うこと。GeneMarkのランを検証するには、例えばダウンロードしたperl gmes_linux_64の外でテストランする。

perl gmes_linux_64/gmes_petap.pl --sequence BRAKER/example/genome.fa --ES --cores 8 --max_intron 100 --min_gene_in_predict 100

エラーが起きるならperlのライブラリが配列@INCに含まれるパスに存在しない可能性がある。perlモジュールはcpanm(apt-getで導入できる)を使ってインストールすれば楽。

cpanmはーLでモジュールの導入パスを指定できる。このとき、ーLで指定したディレクトリのさらにサブディレクトリであるlib/perlに実際には入る。例えば@INCに/root/perl5/lib/perl5/が含まれ、ここにcpanmでperlモジュールをインストールしたいなら、ーLを”-L /root/perl5”と指定する。 Hash::Mergeの導入なら、

cpanm -L /root/perl5 Hash::Merge

と打つ。これで/root/perl5/lib/perl5/にHash/ができ、中にMerge.pmがビルドされモジュールが見えるようになる(ビルドにはmakeが必要)。

 

export PERL_CPANM_OPT="--local-lib=~/perl5"
export PATH=$HOME/perl5/bin:$PATH;
export PERL5LIB=$HOME/perl5/lib/perl5:$PERL5LIB;

 

 

関連

複数のロングリードドラフトアセンブリを使って連続性の高いアセンブリを得る GALA

 

 高品質のゲノムアセンブリは、遺伝学や医学研究の分野で幅広く応用されている。しかし、現在のワークフローでは、ギャップのない染色体スケールのアセンブリを実現することは非常に困難である。ここでは、preliminaryなアセンブリやキメラを含む生データからミスアセンブリを識別し、染色体スケールのリンケージグループにデータを分割する多層コンピュータグラフを用いた、chromosome-by-chromosomeなアセンブリ戦略を提案する。各リンケージグループの後続の独立したアセンブリは、通常、既存のワークフローを悩ませるミスアセンブリエラーから解放されたギャップフリーアセンブリを生成する。この柔軟なフレームワークは、Pacbio、Nanopore、Hi-C、遺伝地図などの様々な技術からのデータを統合して、ギャップフリーな染色体スケールのアセンブリを生成することも可能である。GALAを用いて、公開されているデータセットからPacbioとNanoporeのシーケンスデータを組み合わせてC.elegansとA.thalianaのゲノムをde novoでアセンブルした。また、ヒトゲノムの2本の染色体をギャップフリーで組み立てることで、GALAの適用性を実証した。さらに、GALAはPacbioのhigh-fidelityなロングリードに対しても有望な性能を示した。この方法は、複数のデータソースと複数の計算ツールを用いてゲノムを簡単にアセンブルすることができ、de novoゲノムアセンブル技術の適用の障壁を克服することができる。 

 

インストール

condaを使ってpython2.7の仮想環境を作ってテストした(ホストOS; ubuntu18.04LTS)

依存

  • Minimap2
  • bwa
  • samtools
  • python2.7
  • canu

Github

conda create -n gala python=2.7 -y
conda activate gala
conda install -c bioconda minimap2 bwa samtools canu -y

git clone https://github.com/ganlab/gala.git
cd gala/
#必要ならパスを通す
sudo ./install

> gala -h

$ gala -h

usage: gala -h  [options] <draft_names & paths> <fa/fq> <reads> <platform>

 

GALA Gap-free Long-reads Assembler

 

positional arguments:

  draft_names           Draft names and paths [required]

  input_file            input type (fq/fa) [required]

  reads                 raw/corrected reads [required]

  sequencing_platform   pacbio-raw pacbio-corrected nanopore-raw nanopore-

                        corrected [required]

 

optional arguments:

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

  -a [ASSEMBLER [ASSEMBLER ...]]

                        Chr-by_Chr assembler (canu flye miniasm) [default

                        canu]

  -b Alignment block length [default 5000]

  -p Alignment identity percentage [default 70%]

  -l lowest number of misassemblies indecator [default 1]

  -c Shortest contig length [default 5000]

  -k Mis-assembly block [default 175]

It is better to extend the misassembly block in case of

unpolished assemblies or expected mis-assemblies

in highly repetative regions (5000-10000)

  -q Mapping quality [default 20]

  -f Output files name [default gathering]

  -t cut on a threshold passed by -u [default False]

  -u threshold cut value [default 3]

  -o output files path [default current directory]

  -v, --version         show program's version number and exit

 

 

実行方法

一括して行うモードと、1プロセスごとに進めるモードがある。一括して行うにはgalaコマンドを使う。

ドラフトアセンブリのリストファイルを指定する。リストは

draft_1=path/to/draft_fasta_file
draft_2=path/to/draft_fasta_file
draft_3=path/to/draft_fasta_file

のような形式になっている必要がある。リストに加え、ロングリードを指定する。ここではont-rawのロングリードを指定。

gala ./list fq ont_reads.fq nanopore-raw
  • sequencing_platform   pacbio-raw pacbio-correctednanopore-raw | nanopore-corrected [required]

エラーになる。ランできるようになったら追記します。パフォーマンスについてはプレプリント表1や図3で少し触れられています。

引用

GALA: gap-free chromosome-scale assembly with long reads

Mohamed Awad, Xiangchao Gan

bioRxiv, Posted May 16, 2020

 

関連


 

 

 

 

 

LongAGE

 

 1塩基のブレークポイント分解能で構造変異(SV)の正確な位置を定義することは、アライメントのギャップが大きいために困難な問題である。これまでは、AGE(Alignment with Gap Excision)を用いることで、1塩基分解能でSVのブレークポイントを定義することができたが、長い配列のペアをアライメントする場合、AGEは膨大な量のメモリを必要とする。この問題を解決するために、著者らは古典的なHirschbergアルゴリズムをベースにしたメモリ効率の良い実装であるLongAGEを開発した。LongAGEを使ってPacBioの10Kbpを超える長さのリードのセグメント重複に埋め込まれたSVのブレークポイントを解決し、LongAGEを実証した。さらに、同じ遺伝子座の欠失と重複について異なるブレークポイントが観察され、このようなマルチアレルコピーナンバーバリアント(mCNV)が2つ以上の独立した先祖代々の突然変異から生じることを直接的に証明している。LongAGEはC ++で実装されており、Githubhttps://github.com/Coaxecva/LongAGE)で利用できる。

 

インストール

Github

git clone https://github.com/Coaxecva/LongAGE.git
cd LongAGE
make

> ./long_age_align -h

# ./long_age_align -h

Usage:

long_age_align

[-version] [-allpos]

[-indel|-tdup|-inv|-invl|-invr]

[-match=value:1] [-mismatch=value:-3]

[-go=value:-7] [-ge=value:-1]

[-both] [-revcom1] [-revcom2]

[-coor1=start-end] [-coor2=start-end] file1.fa file2.fa

 

テストラン

リファレンスとロングリードを指定する。デフォルトはindel検出モード(-indel)になっているが、フラグを立てることで変更可能。

#indel mode
long_age_align -indel part_chr19.fa seq_first.fa
# TDUP mode
long_age_align -tdup part_chr19.fa seq_second.fa
  •  -indel|-tdup|-inv|-invl|-invr

出力

 

MATCH = 1, MISMATCH = -3, GAP OPEN = -7, GAP EXTEND = -1, INDEL

 

First  seq [1,33001] =>     33001 nucs 'chr19:54210000-54243000'

Second seq [1,17004] =>     17004 nucs 'm54015_171028_124452/40502023/4621_21625'

 

Score:        4819

Aligned:     13127        nucs

Identic:     11683 ( 89%) nucs =>      9173 ( 89%)      2510 ( 88%)

Gaps:          855 (  7%) nucs

 

Alignment:

 first  seq =>  [   1, 9914] EXCISED REGION [29680,32426]

 second seq =>  [4106,14095] EXCISED REGION [14112,16859]

 

EXCISED REGION(S):

 first  seq =>     19765 nucs [9915,29679]

 second seq =>        16 nucs [14096,14111]

 

Identity at breakpoints: 

 first  seq =>         3 nucs [9914,9916] to [29679,29681]

 second seq =>         0 nucs

Identity outside breakpoints: 

 first  seq =>         1 nucs [9914,9914] to [29680,29680]

 second seq =>         1 nucs [14095,14095] to [14112,14112]

Identity inside breakpoints: 

 first  seq =>        21 nucs [9915,9935] to [29659,29679]

 second seq =>         0 nucs

 

        1 GTTCTAATTCTAGAGGCAAGCCCCTTCCATGT-CCTGAGCTCTGTAAT-GCATCTTTT-C 57

          |||||||||||||||||||||||||||||||| |||.||||||||||| ||.|||||| |

     4106 GTTCTAATTCTAGAGGCAAGCCCCTTCCATGTTCCTTAGCTCTGTAATCGCCTCTTTTTC 4165

 

       58 TTTCATGAGCCTTGCGAT-CAGGCGATGTTTATTCAGTGGTTACCACATCCAGGCATGCT 116

          |||||||||| ||||||| ||||...||||| |||||||||||||||||||||||.  ||

     4166 TTTCATGAGC-TTGCGATTCAGGGCGTGTTT-TTCAGTGGTTACCACATCCAGGCT--CT 4221

 

      117 GCCAGGAGGAGGGGAGTCGTGGGTGAAGCTGATAGGATTCCTGCTGGACTCACAGAGCCT 176

          ||||.||||||.| ||||||||| ||||||||||.|||||||||||| |||||.||||||

     4222 GCCATGAGGAGTG-AGTCGTGGG-GAAGCTGATATGATTCCTGCTGG-CTCACCGAGCCT 4278

 

      177 GGGTTAATGACACATTACCCATGTTTAGATAGGAGGTAATTCTGCTCCGGTTTCGACAAG 236

          |||||| ||||||||||||||| ||||||||||  ||||||||||..|||||||||||||

     4279 GGGTTA-TGACACATTACCCAT-TTTAGATAGG--GTAATTCTGCCTCGGTTTCGACAAG 4334

 (以下略)

 

 

 

引用

LongAGE: defining breakpoints of genomic structural variants through optimal and memory efficient alignments of long reads
Quang Tran, Alexej Abyzov
Bioinformatics, Published: 10 August 2020

完全性、正確性、連続性を考量してゲノムアセンブリを評価する PDR

 

 既存のゲノムアセンブリ評価指標は、ゲノムアセンブリの品質の特定の側面についての限られた知見しか提供しておらず、時にはお互いに意見が合わないこともある。アセンブリ間の統合的な比較をより良くするために、著者らはここで新しいゲノムアセンブリ評価指標PDRを提案する。この評価指標は、遺伝学研究における共通の関心事に由来しており、完全性、連続性、正確性を考慮に入れている。また、PDRの計算を高速化するための近似実装を提案する。
 公開されているデータセットに対する評価結果は、ゲノムアセンブリの品質を統合的に評価するPDRの能力を肯定している。実際、このことはその定義によって保証されている。また、近似によって生じる誤差は極めて小さく、無視できる程度であることが示された。

 

1.1 Continuity

 ゲノムアセンブリー評価のための標準的な指標は原則として存在しない。それにもかかわらず、いくつかのメトリックは、そのような仕事のほとんどすべてで共通して使用されている。典型的な例はN50であり、これは、アセンブリの少なくとも50%の長さが、それと等しいかそれよりも長いscaffolds(またはコンティグ)によって寄与されている長さとして定義されている。この統計量は、広く使われているだけでなく、批判もされている。初期のチャレンジでは、サイズの異なるアセンブリ、特にこれらのアセンブリが同一のサンプルから得られた場合には、比較が不公平になる可能性があると主張があった。そこで提案されたのがNG50で、これはアセンブリサイズを推定ゲノムサイズに置き換えることを除いてはN50と似ている(Earl et al., 2011)。しかし、NG50であっても、場合によっては誤解を招く可能性がある。論文図1は、青のアセンブリが赤のアセンブリよりも全体的に優れているが、NG50とN50が低い例を示している。QUAST(Gurevichら、2013)は、50%の支配を避けるためにNGx(Nx)プロットを提供することによって、この問題に対処した。
 しかし、このプロットは、NG50における別の問題を解決することができない:より大きなコンティグをもたらし、したがってNG50を膨らませるミスアセンブリ(a.k.a misjoinまたは構造的エラー)。これを解決するために、リファレンスゲノムが利用可能な場合、多くのベンチマークや評価ツールでは、NG50のカウントにおいてコンティグ長を置き換えるためにアライメントブロック長を使用していた。これらの指標には、コンティグパスNG50(Earlら、2011)、補正Nx(Salzbergら、2012)、正規化N50(Mäkinenら、2012)、およびNGA50(Gurevichら、2013)がある。これらは定義と実装が若干異なるが、原理は同じである:誤ったコンティグを誤ってアセンブリされた点で切断することである。この一連のメトリクスは、通常、連続性の評価に使用される。
 これら以外にも、限られた研究でしか使用されていないにもかかわらず、いくつかのコンティギュイティ指標がある。例えば、E-sizeはGAGEで提案されたもので、ランダムな位置に位置するコンティグの長さの期待値として定義されている。N50と同様に、この指標はコンティグ長の分布を反映しているに過ぎず、ミスアセンブリによって上昇する可能性がある。リファレンスゲノムを用いて、別のメトリックU50(CastroおよびNg、2017)は、N50メトリックに固有のいくつかの制限を回避することを目的として、重複する配列を除去することによって、ユニークで標的特異的なコンティグを同定した。アセンブラソン(Earlら、2011;Bradnamら、2013)では、CC50は明確な考えとして定義された。染色体内の位置のペアは、正しい順序と同じコンティグで識別されていれば、正しく連続している(CC)と言える。CCペアの割合は分離距離とともに減少する。CC50は、この長さで区切られた50%のペアがCCになる長さである。計算に時間がかかるため、サンプリングベースの方法でしか推定できまなかった。

1.2 Completeness
 完全性はゲノムアセンブリの評価におけるもう一つの次元である。いくつかのベンチマークでは、連続性のみが評価され、完全性はその一部として考慮されていた。実際には、完全性はアセンブリによって引き起こされる損失に焦点を当てているのに対し、連続性は局所的なコンテキストの再構築を反映している。完全性の普遍的なメトリックは、アライメントカバレッジである。これは非常にわかりやすい指標だが、時として識別力が低いことがある。リファレンスゲノムが利用できない場合、完全性はCEGMA(Parraら、2007)およびBUSCO(Simãoら、2015)によって評価することができる。これらは、保存されたシングルコピーオルソログのセットを収集し、アセンブリがそれらを含むかどうかをテストする。言い換えれば、これらが反映しているのは遺伝子空間の完全性であり、アセンブルされたゲノムの正確な完全性ではない。つまり、これらは実際にはサンプリングテストである。
1.3 Correctness
 ゲノムアセンブリの評価におけるもう一つの側面は、正しさである。一般的に使用されているメトリクスは、連続性のメトリクスよりもわかりやすく、時間の経過とともに改善されることはあまりない。これらのメトリクスには、通常、1塩基エラー(ミスマッチ)、indel、ミスアセンブリが含まれる。ミスアセンブルは最も有害なタイプとみなされ、正しさを評価するためのメトリクスで広く使われている。このタイプをより良くプロファイルするために、QUASTではさらに以下のように分類されている。(a) relocation、すなわち、隣接する配列が同じ染色体に配置されているが互いに離れている位置にある、(b)逆位、すなわち、隣接する配列が同じ染色体に配列しているが反対側のストランドに配置されている、(c)転座、すなわち、隣接する配列が異なる染色体に配置されている、に分類される。これらは正確性を評価するには十分だが、デノボアセンブリでは通常利用できないリファレンスゲノムが必要となるため、利用には制限がある。さらに、これらのカウントは、アセンブリの数に依存しているため、アセンブリの大きさや効果を反映していない。そこで、REAPR(Huntら、2013)およびLAP(Ghodsiら、2013)は、リファレンスゲノムが存在しない場合に、アセンブリとサンプルのシーケンシングリードとの間の整合性を確認することによってアセンブリを評価するために提案された。REAPRは、長いインサートサイズ(1000bp以上)のペアエンドリードをマッピングしたものを必要とし、そのペアリングとマッピング情報を評価エビデンスとして利用した。そのため、REAPRの性能はリードの品質やインサートサイズに左右される。LAPでは、品質をアセンブルされた配列のリードを観察する条件付き確率と定義している。しかし、それは同じリードセットに由来するアセンブリを比較する場合にのみ適用される。REAPRやLAPと同様に、FRCurve (Narzisi and Mishra, 2011)は、ミスアセンブリの特徴を検出するためにリードのレイアウト情報を利用する。そして、それは、ミスアセンブリ特徴数(X軸)の与えられた数内のコンティグの最大総長(Y軸)を示す特徴-応答曲線をプロットする。実際、FRCurveは正しさだけでなく、コンティグ性も評価する。定量的なメトリックを与えるために、プロットは、補正されたN50を計算するために使用される。NGA50とは異なり、補正されたN50はN50の欠点を完全に回避しているわけではない。大きなコンティグの中にいくつかのミスアセンブリがあっても、補正後のN50は大きく膨らむ可能性がある。FRCurveのもう一つの問題点は、すべてのエラー特徴(アライメントブレークポイント、low depth、異常なリード方向など)が、下流の解析での影響に関係なく等しく重み付けされていることである。
1.4 Need for an overall metric
 連続性と正しさの間にはトレードオフがある。具体的には、アグレッシブなアセンブリの中には、正確性を犠牲にしても連続性に優れているものもあれば、エラーが少なくコンティグが短いものもある。ダイバージェンスは、ほとんどの場合、繰り返し領域でコンティグが伸びている間に発生する。リピート領域が長すぎて1回のリードではカバーできない場合、アセンブラはコンティグを拡張するために複数の選択肢を持っているが、どれが正しい選択肢なのかはわからない。保守的なアセンブラは、決定を下すのに十分な証拠がない場合に伸長を停止し、結果として断片化したアセンブラリを生成する。対照的に、アグレッシブなアセンブラは、ミスアセンブルのリスクを冒しながらも、微妙な手がかりを得て、シーケンスを拡張し続ける。異なるアセンブラは、拡張を継続するかどうかを判断するための独自の戦略を持っており、継続する場合はどちらの選択が良いかを判断する。
 ベンチマークや評価ツールでは、通常、アセンブリの多次元的な比較やグローバルなプロファイルを提供するために、連続性、完全性、正しさ(C3)のそれぞれをいくつかのメトリクスで調査している。しかし、Haiminenら(2011)が批判しているように、異なるアセンブリメトリクスの表が提供されると、各メトリクスが独自のフロントランナーを持っているため、アセンブリ間の比較が複雑化する。このような過度に詳細な情報は、実際のexperienceの指針にはならない。どのアセンブリが全体的に優れているかを判断するのはまだ難しい。これに対処するためには、連続性、完全性、正確性の3つの側面をすべて統合した総合的なメトリックまたはスコアが望ましい。そこで、このギャップを埋めるために、ゲノムアセンブリの評価のための新しい指標を提案する。(以下略)

 

論文の式(8)の近似がPDR(論文参照

 

インストール

ubuntu18.04LTSでテストした。

依存

  • PDRi only needs Java(1.8 or above)

ここではJDK1.8を導入

apt update && apt install openjdk-8-jdk

GIthub

リリースからPDRi.jarをダウンロードする。

> java -jar PDRi.jar

$ java -jar PDRi.jar 

Usage: pdmega [OPTIONS] REFERENCE ASSEMBLY

 

Options:

  --threads INT  Threads to use [default: CPU core]

  -k INT         Block size [default: 1000]

  -d PATH        Temporary folder for intermediate files [default: PDRTmp]

  -a TEXT        Executable path of aligner (BWA or minimap2) [default: bwa]

  -e INT         Maximum offset for two alignment segment to be jointed

                 [default: 0]

  -m INT         Minimum chromosome length (in bp) to summarize and report

                 alignment statistics. This doesn't change PDR result.

                 [default: 1% genome]

  -h, --help     Show this message and exit

 

Arguments:

  REFERENCE  Reference genome

  ASSEMBLY   Assembly to evaluate

 

 

実行方法

リファレンスゲノムとアセンブリ配列を指定する。

java -jar PDRi.jar reference_geonme.fasta input_assembly.fasta

 出力

====== Finished ======

Genome payload: 3954589.0

PDR Total:      4924756.2547394745

PDR Ratio:      3.1490679542361645E-7

 

引用

PDR: a new genome assembly evaluation metric based on genetics concerns
Luyu Xie, Limsoon Wong

Bioinformatics, Published: 06 August 2020

 

関連

 

publication readyなggplot2 プロット出力を行う ggpubr

2020 8/9 誤字修正

 

ggplot2' パッケージは R でのエレガントなデータ可視化のための優れた柔軟性を持っているが、デフォルトで生成されるプロットは、出版前にいくつかのformatingを必要とする。さらに、'ggplot'をカスタマイズするための構文は不透明であり、高度なRプログラミングスキルを持たない研究者にとっては難易度が高くなる。ggpubr' は、'ggplot2'ベースのpublication readyプロットを作成し、カスタマイズするための使いやすい関数を提供する。 

 

インストール

Rstudio最新版を使ってテストした。

Github 

本体 Github 

install.packages("ggpubr")

#最新版
if(!require(devtools)) install.packages("devtools")
devtools::install_github("kassambara/ggpubr")

 

実行方法

1、ライブラリのロード

library("ggpubr")

 

2、データの読み込み

次のような整然構造のフォーマットを読み込む(excelライクな形式は*1参照)。

  genome GC
1 genome1 77
2 genome1 64
3 genome1 70
4 genome1 55
5 genome1 77
6 genome1 76
7 genome1 76
8 genome1 79
9 genome1 71
10 genome2 60
11 genome2 61
12 genome2 59
13 genome2 59
14 genome2 63
15 genome2 64
16 genome2 61
17 genome2 57
18 genome3 91
19 genome3 90
20 genome3 90
21 genome3 97
22 genome3 85
23 genome3 91
24 genome3 83
25 genome3 88

excelライクなフォーマットと違って欠損値(欠損セル)問題が起きない。

 

コピペ読み込み。上の表をコピーして下のコマンドを実行する。

#mac
input <- read.table(pipe("pbpaste"), header = TRUE)

#windows
input <- read.table("clipboard"), header = TRUE)

またはファイルinput.tsvに保存し、データフレームに読み込む。  

#カレント
input <- read.table("input.tsv", header=T, sep="\t")

#ファイルがカレントにないならフルパス指定か. セパレータがコンマならsep=","
input <- read.table("/home/kazu/data/input.csv", header=T, sep=",")

 

3、視覚化

violin plot

ggviolin(input, x = "genome", y = "GC", fill = "genome",
palette = c("#00AFBB", "#E7B800", "#FC4E07"),
add = "boxplot", add.params = list(fill = "white"))+
stat_compare_means(comparisons = my_comparisons, label = "p.signif")+
stat_compare_means(label.y = 100)

f:id:kazumaxneo:20200301142355p:plain

 

box plot

input <- read.table(pipe("pbpaste"), header = TRUE)

ggboxplot(input, x = "genome", y = "GC", color = "genome", palette =c("#00AFBB", "#E7B800", "#FC4E07"), add = "jitter", shape = "genome") 

f:id:kazumaxneo:20200807122419p:plain

 

p-valueを追加

p <- ggboxplot(input, x = "genome", y = "GC", color = "genome", palette =c("#00AFBB", "#E7B800", "#FC4E07"), add = "jitter", shape = "genome")

#一度描画
p

#比較対象を定義
my_comparisons <- list( c("genome1", "genome2"), c("genome2", "genome3"), c("genome1", "genome3") )

#Piarwise比較してp-valueを追加
p + stat_compare_means(comparisons = my_comparisons)

 

f:id:kazumaxneo:20200807123049p:plain

 

#比較対象を定義
my_comparisons <- list( c("genome1", "genome2"), c("genome2", "genome3"), c("genome1", "genome3") )

#global p-valueを求めた結果も追加.label.y = 120でy=120の位置に追加
p + stat_compare_means(comparisons = my_comparisons) + stat_compare_means(label.y = 120)  

f:id:kazumaxneo:20200807123409p:plain

 

追記

7つ以上データ(ここではgenome)がある場合、手動で色を指定する必要がある。

ggviolin(input, x = "genome", y = "GC", fill = "genome",
palette =c("#CC0000", "#CC6600", "#CCCC00", "#66CC00", "#00CCCC", "#0066CC", "#CC00CC", "#606060"),
add = "boxplot", add.params = list(fill = "white"))

f:id:kazumaxneo:20200809192607p:plain

 

Githubやマニュアルには他にも様々な作図例が掲載されています。確認して下さい。

引用

Ggpubr: ‘Ggplot2’ Based Publication Ready Plots

Alboukadel Kassambara

2018 R Package Version 0.2. Available online: https://CRAN.R-project.org/package=ggpubr

 

 *1

下の様なexcelライクなフォーマットの方が見やすいが、全ての列で行数が同じなどの制約がある。さらに、excelライクなフォーマットだとRに読み込んでから一手間かかる。

 

genome1 genome2
77 67
64 54
70 60
55 45
77

67

のようなデータのbar plotを出力するなら

#読み込み
x <- read.table(pipe("pbpaste"), header = TRUE)

#名前、平均、SD
colnames(x) <- c("genome1", "genome2") #カラム名を再アサイン(スキップ可)
xmean <- apply(x, 2, mean) #平均
xsd <- apply(x, 2, sd) #S.D

#プロット
plot <- barplot(bindmean, xlab = "name", ylab = "logFPKM", ylim = c(0, max(xmean + xsd)))

f:id:kazumaxneo:20200301141736p:plain

さらにエラーバーをつける。

arrows(plot, xmean - xsd, b, xmean + xsd, code = 3, lwd = 1, angle = 90, length = 0.1)

 f:id:kazumaxneo:20200301141750p:plain

このように、視覚化まで時間がかかる。整然フォーマットのデータを準備し、さらにggpubrを使うことで煩雑なコードを削減できる。

 

 

参考

 

関連

 

(ヒトゲノム)高速かつ精度の高いロングリードのSVコーラー cuteSV

 

 構造変化(SV)とは、欠失、挿入、逆位、重複、転座などのゲノムリアレンジメントで、その大きさが50 bpを超えるものを指す。ヒトゲノム上で最大のdivergencesとして、SV はヒトの疾患(遺伝性疾患やガンなど)、進化(遺伝子欠損やトランスポゾン活性など)、遺伝子制御、その他の表現型(交配や本質的な生殖隔離など)と密接に関連している。

 ショートリードベースのSVコールアプローチを開発するための努力がなされてきた。それらの多くは、リードデプス、不一致リードペア、スプリットリードアラインメント、ローカルアセンブリ、またはそれらの組み合わせなどの方法を用いており、1000 Genomes Projectのような大規模なゲノム研究において重要な役割を果たしてきた。しかし、これらのツールは比較的短いリード長であるため、高感度なSV検出には限界があり、偽陽性も存在する。

 Pacific Bioscience (PacBio)や Oxford Nanopore Technology (ONT)などのロングリードシークエンシング技術の急速な発展に伴い、ロングレンジスパニング情報は、より高分解能でより包括的に SV を検出する機会を提供している。しかし、高いシーケンシングエラー率(通常5~20%)とリード長(平均10kbp以上)に対応するためには、新たな計算アプローチが必要である。これまでの研究では、主に2種類のアプローチが採用されてきた。

 デノボアセンブリーベースのアプローチ は、より長いゲノム配列(コンティグやスキャホールド)にリードをアセンブリし、アセンブリした配列とリファレンス配列とのアラインメントからSVを発見することを目的としている。このようなアプローチは、リードアライメントに基づくアプローチに比べて、リファレンスの影響を受けにくく、特にリードアライメントのアーチファクトがない。しかし、これらのアプローチは通常計算量が多く、大規模なゲノムのハプロタイプ配列の再構成にはまだいくつかの困難がある。

 リードアライメントに基づくアプローチは、リードをリファレンスに対して直接アライメントし、アライメント結果を解析することでSVを検出する。このようなアプローチは、感度を欠くことなく計算リソースに対してより費用対効果が高く、ロングリードベースのSVコーラーで広く使用されている。いくつかのリードアライメントベースのSVコーラーが提案されており、PB-Honey [ref.31]、SMRT-SV [ref.32]、Sniffles [ref.33]、PBSV (https://github.com/PacificBiosciences/pbsv)、およびSVIM [ref.34]などがある。これらの手法では、発散性の高いアラインメントを持つ局所的なゲノム領域の同定、クリッピングされたリード部分の局所的なアセンブリと再アラインメント、SV-スパンニングシグネチャクラスタリングなど、リードアラインメントによって暗示されるSVの証拠を見つけるために様々な方法を使用している。さらに、通常、BLASR [ref.36]、NGMLR [ref.33]、Minimap2 [ref.37]、PBMM2 (https://github.com/PacificBiosciences/pbmm2) などの最先端のロングリードアライナーがリードアライメントに使用されていた。

 しかし、リードアライメントに基づくSVコールは、依然として非自明である。シーケンシングエラーが多く、複雑なSVが存在する状況下では、SVブレークポイント周辺のリードのアライメントはキメラ的で不均質であり、通常は感度や精度が低い。そのため、リードアラインメントから得られるSVシグネチャは非常に複雑であり、様々な種類のSVに対して感度の高い検出を行うためには、それらを収集・解析することが困難である。主に、最先端のツールでは、以下のような技術的な課題がある。(1) 全体的に感度がまだ満足のいくものではない。(2) いくつかのアプローチ(rMETL [ref.38]、rCANID [ref.39]、npInv [ref.40]など)は、特定の設計のために、サブセットまたは特定のクラスのSVしか検出できない。(3) いくつかのアプローチ(PBSVやSMRT-SVなど)は、まだ時間がかかり、スケーリング性能が良くないため、多くの大規模データセットには適していない可能性がある;そして(4) いくつかのアプローチ(SMRT-SVやPB-Honeyなど)は、1種類のシーケンシングデータしかサポートしていない(例えばPacBioのリードのみ)ものもある。このような欠点が、ロングリードシーケンシングデータを広く利用する上でのボトルネックとなっている。 

 ここでは、いくつかの利点を持つ多用途なSV検出アプローチであるcuteSVを紹介する。(1) cuteSVは、最先端のSVコーラーよりも高いSV検出率を持っている。特に、カバレッジの低いデータに対しても、精度を落とさずに高感度に検出することができる。(2) cuteSVは、主流のロングリードシーケンシングプラットフォームで作成された様々なエラー率のデータセットをサポートしており、様々なタイプのSV(欠失、挿入、重複、逆位、転座を含む)を検出することができる。(3) cuteSVは、RAMの使用量が少なく、最先端のアプローチと同等の速度を持っている。さらに重要なことは、CPUのスレッド数に応じてほぼ直線的に高速化できるスケーラビリティを持っている。これらの特徴により、cuteSVは大規模なデータ解析に適しており、最先端のゲノミクス研究への可能性を秘めている。

cuteSVの概要
ソートされたBAMファイルを入力として、cuteSVは、アライメントされたリード中の大きな挿入/欠失や分割アライメントをSVシグネチャクラスターとして抽出し、解析してSVをコールする。このアプローチには大きく分けて以下の3つのステップがある。

ステップ1:cuteSVは複数のシグネチャ抽出手法を用いて、様々なタイプのSVのシグネチャを網羅的に収集する。さらに、挿入と欠失をヒューリスティックに組み合わせて、脆弱なアライメントから本物のSVの証拠を回収する。

ステップ2: cuteSVは、特別に設計されたクラスタリングと精緻化のアプローチを用いて、キメラアラインメントされたリードを局所領域にクラスタリングし、さらにクラスタを精緻化して、SVのシグネチャヘテロ接合型SVと正確に区別する。

ステップ3:cuteSVは、いくつかのオーダーメイドのルールを使用して、洗練されたSVシグネチャクラスターに基づいてSVコールとジェノタイピングを実行する。

 

 

インストール

依存

  • python3
  • pysam
  • Biopython
  • cigar
  • numpy

Github

#bioconda (link)
conda install -c bioconda cutesv -y

#pip
pip install cuteSV

cuteSV -h

$ cuteSV -h

usage: cuteSV [-h] [--version] [-t THREADS] [-b BATCHES] [-S SAMPLE]

              [--retain_work_dir] [-p MAX_SPLIT_PARTS] [-q MIN_MAPQ]

              [-r MIN_READ_LEN] [-md MERGE_DEL_THRESHOLD]

              [-mi MERGE_INS_THRESHOLD] [-s MIN_SUPPORT] [-l MIN_SIZE]

              [-L MAX_SIZE] [-sl MIN_SIGLENGTH] [--genotype]

              [--gt_round GT_ROUND]

              [--max_cluster_bias_INS MAX_CLUSTER_BIAS_INS]

              [--diff_ratio_merging_INS DIFF_RATIO_MERGING_INS]

              [--diff_ratio_filtering_INS DIFF_RATIO_FILTERING_INS]

              [--max_cluster_bias_DEL MAX_CLUSTER_BIAS_DEL]

              [--diff_ratio_merging_DEL DIFF_RATIO_MERGING_DEL]

              [--diff_ratio_filtering_DEL DIFF_RATIO_FILTERING_DEL]

              [--max_cluster_bias_INV MAX_CLUSTER_BIAS_INV]

              [--max_cluster_bias_DUP MAX_CLUSTER_BIAS_DUP]

              [--max_cluster_bias_TRA MAX_CLUSTER_BIAS_TRA]

              [--diff_ratio_filtering_TRA DIFF_RATIO_FILTERING_TRA]

              [BAM] output work_dir

 

Long read based fast and accurate SV detection with cuteSV.

 

Current version: v1.0.6

Author: Tao Jiang

Contact: tjiang@hit.edu.cn

 

Suggestions:

 

For PacBio CLR/ONT data:

--max_cluster_bias_INS 100

--diff_ratio_merging_INS 0.2

--diff_ratio_filtering_INS 0.6

--diff_ratio_filtering_DEL 0.7

For PacBio CCS(HIFI) data:

--max_cluster_bias_INS 200

--diff_ratio_merging_INS 0.65

--diff_ratio_filtering_INS 0.65

--diff_ratio_filtering_DEL 0.35

 

 

 

positional arguments:

  [BAM]                 Sorted .bam file form NGMLR or Minimap2.

  output                Output VCF format file.

  work_dir              Work-directory for distributed jobs

 

optional arguments:

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

  --version, -v         show program's version number and exit

  -t THREADS, --threads THREADS

                        Number of threads to use.[16]

  -b BATCHES, --batches BATCHES

                        Batch of genome segmentation interval.[10000000]

  -S SAMPLE, --sample SAMPLE

                        Sample name/id

  --retain_work_dir     Enable to retain temporary folder and files.

 

Collection of SV signatures:

  -p MAX_SPLIT_PARTS, --max_split_parts MAX_SPLIT_PARTS

                        Maximum number of split segments a read may be aligned

                        before it is ignored.[7]

  -q MIN_MAPQ, --min_mapq MIN_MAPQ

                        Minimum mapping quality value of alignment to be taken

                        into account.[20]

  -r MIN_READ_LEN, --min_read_len MIN_READ_LEN

                        Ignores reads that only report alignments with not

                        longer than bp.[500]

  -md MERGE_DEL_THRESHOLD, --merge_del_threshold MERGE_DEL_THRESHOLD

                        Maximum distance of deletion signals to be merged.[0]

  -mi MERGE_INS_THRESHOLD, --merge_ins_threshold MERGE_INS_THRESHOLD

                        Maximum distance of insertion signals to be

                        merged.[100]

 

Generation of SV clusters:

  -s MIN_SUPPORT, --min_support MIN_SUPPORT

                        Minimum number of reads that support a SV to be

                        reported.[10]

  -l MIN_SIZE, --min_size MIN_SIZE

                        Minimum size of SV to be reported.[30]

  -L MAX_SIZE, --max_size MAX_SIZE

                        Maximum size of SV to be reported.[100000]

  -sl MIN_SIGLENGTH, --min_siglength MIN_SIGLENGTH

                        Minimum length of SV signal to be extracted.[10]

 

Computing genotypes:

  --genotype            Enable to generate genotypes.

  --gt_round GT_ROUND   Maximum round of iteration for alignments searching if

                        perform genotyping.[500]

 

Advanced:

  --max_cluster_bias_INS MAX_CLUSTER_BIAS_INS

                        Maximum distance to cluster read together for

                        insertion.[100]

  --diff_ratio_merging_INS DIFF_RATIO_MERGING_INS

                        Do not merge breakpoints with basepair identity more

                        than [0.2] for insertion.

  --diff_ratio_filtering_INS DIFF_RATIO_FILTERING_INS

                        Filter breakpoints with basepair identity less than

                        [0.6] for insertion.

  --max_cluster_bias_DEL MAX_CLUSTER_BIAS_DEL

                        Maximum distance to cluster read together for

                        deletion.[200]

  --diff_ratio_merging_DEL DIFF_RATIO_MERGING_DEL

                        Do not merge breakpoints with basepair identity more

                        than [0.3] for deletion.

  --diff_ratio_filtering_DEL DIFF_RATIO_FILTERING_DEL

                        Filter breakpoints with basepair identity less than

                        [0.7] for deletion.

  --max_cluster_bias_INV MAX_CLUSTER_BIAS_INV

                        Maximum distance to cluster read together for

                        inversion.[500]

  --max_cluster_bias_DUP MAX_CLUSTER_BIAS_DUP

                        Maximum distance to cluster read together for

                        duplication.[500]

  --max_cluster_bias_TRA MAX_CLUSTER_BIAS_TRA

                        Maximum distance to cluster read together for

                        translocation.[50]

  --diff_ratio_filtering_TRA DIFF_RATIO_FILTERING_TRA

                        Filter breakpoints with basepair identity less than

                        [0.6] for translocation.

 

実行方法

bamファイルを指定する。

mkdir work_dir
cuteSV -t 20 input.bam output.vcf work_dir

vcf4.2 formatのSVコールが出力される。

引用

Long-read-based human genomic structural variation detection with cuteSV

Tao Jiang, Yongzhuang Liu, Yue Jiang, Junyi Li, Yan Gao, Zhe Cui, Yadong Liu, Bo Liu & Yadong Wang
Genome Biology volume 21, Article number: 189 (2020)

 

(microbial genomes)低分子量タンパク質のアノテーションを付ける SmORFinder

 

 Sberroら(2019)が行った最近の研究により、ヒトマイクロバイオーム内に存在するスモールタンパク質の広大な未踏空間が明らかになった。現在のところ、これらの小さなオープンリーディングフレーム(smORF)は既存のリファレンスゲノムではアノテーションされておらず、標準的なゲノムアノテーションツールでは正確に予測することができない。本研究では、Sberroらによって同定されたものをもとにスモールタンパク質を予測するアノテーションツールSmORFinderを導入した。 このツールは、各スモールタンパク質ファミリーのプロファイル隠れマルコフモデル(pHMM)と、トレーニングセットでは見られないsmORFファミリーに対してより一般化する可能性のあるディープラーニングモデルを組み合わせたものである。pHMMとディープラーニングモデルの両方の予測を組み合わせることで、より正確なsmORF予測が可能となり、これらの予測されたsmORFはRibo-SeqまたはMetaRibo-Seq翻訳シグナルに富むことがわかった。特徴重要度解析により、ディープラーニングモデルは、Shine-Dalgarno配列を識別し、各コドンのウォブル位置を優先順位付けし、コドンテーブルで見つかったコドンの同義語に強く対応するようにコドンをグループ化することを学習したことが明らかになった。また、その中から、機能未知のコアとなる多くのsmORFを同定している。また、何千ものRefSeq単離ゲノムとHMPメタゲノムのスモールタンパク質アノテーションを事前に計算し、これらのデータをスモールタンパク質アノテーションと解析のための他の有用なツールとともにウェブポータルを通じて利用できるようにした。これらの重要なスモールタンパク質を体系的に同定し、アノテーションを行うことで、研究者は生物学のこのエキサイティングな分野の理解を深めることができる。

 

webサーバー

計算済みの低分子量タンパク質のアノテーション(RefseqとHuman Microbiome Project - HMP)ダウンロードとユーザーゲノムの新規アノテーション、そして低分子量タンパク質の分析が可能になっている。

http://104.154.134.205:3838/DBsmORF/

f:id:kazumaxneo:20200807230525p:plain

 

ローカルでの使用

インストール

ubuntu18.04のpython3.7環境でテストした。

Github

pip install smorfinder

> smorf #紹介はhmmモデルがダウンロードされるため、help表示まで数分以上の時間がかかる。

smorf --help

$ smorf --help

Usage: smorf [OPTIONS] COMMAND [ARGS]...

 

  Command-line tool to predict and annotate small protein sequences in

  genomic sequencing data

 

Options:

  --help  Show this message and exit.

 

Commands:

  single  Run SmORFinder on a complete or draft genome assembly of a single

          species.

  meta    Run SmORFinder on a metagenomic assembly.

smorf single --help

$ smorf single --help

Usage: smorf single [OPTIONS] FASTA

 

  A click access point for the run module. This is used for creating the

  command line interface.

 

Options:

  -o, --outdir TEXT

  -pp, --prodigal-path PATH

  -shp, --dsn1-model-path PATH

  -shp, --dsn2-model-path PATH

  -shp, --smorf-hmm-path PATH

  -hp, --hmmsearch-path PATH

  --force / --no-force            Force overwriting of output directory.

  -idsn1, --dsn1-indiv-cutoff FLOAT

                                  Minimum cutoff necessary to keep prediction

                                  based on DSN1 significance cutoff alone.

                                  Between 0 and 1, default=0.9999

  -idsn2, --dsn2-indiv-cutoff FLOAT

                                  Minimum cutoff necessary to keep prediction

                                  based on DSN2 significance cutoff alone.

                                  Between 0 and 1, default=0.9999

  -iphmm, --phmm-indiv-cutoff FLOAT

                                  Minimum cutoff necessary to keep prediction

                                  based on pHMM significance cutoff alone.

                                  Between 0 and 1, default=1e-6

  -odsn1, --dsn1-overlap-cutoff FLOAT

                                  Minimum cutoff necessary to keep prediction

                                  based on DSN1 significance if both other

                                  models meet their respective cutoffs.

                                  Between 0 and 1, default=0.5

  -odsn2, --dsn2-overlap-cutoff FLOAT

                                  Minimum cutoff necessary to keep prediction

                                  based on DSN2 significance if both other

                                  models meet their respective cutoffs.

                                  Between 0 and 1, default=0.5

  -ophmm, --phmm-overlap-cutoff INTEGER

                                  Minimum cutoff necessary to keep prediction

                                  based on pHMM significance if both other

                                  models meet their respective cutoffs.

                                  Between 0 and 1, default=1

  --help                          Show this message and exit.

smorf meta --help

$ smorf meta --help

Usage: smorf meta [OPTIONS] FASTA

 

  A click access point for the run module. This is used for creating the

  command line interface.

 

Options:

  -o, --outdir TEXT

  -t, --threads INTEGER

  -pp, --prodigal-path PATH

  -shp, --dsn1-model-path PATH

  -shp, --dsn2-model-path PATH

  -shp, --smorf-hmm-path PATH

  -hp, --hmmsearch-path PATH

  --force / --no-force            Force overwriting of output directory.

  -idsn1, --dsn1-indiv-cutoff FLOAT

                                  Minimum cutoff necessary to keep prediction

                                  based on DSN1 significance cutoff alone.

                                  Between 0 and 1, default=0.9999

  -idsn2, --dsn2-indiv-cutoff FLOAT

                                  Minimum cutoff necessary to keep prediction

                                  based on DSN2 significance cutoff alone.

                                  Between 0 and 1, default=0.9999

  -iphmm, --phmm-indiv-cutoff FLOAT

                                  Minimum cutoff necessary to keep prediction

                                  based on pHMM significance cutoff alone.

                                  Between 0 and 1, default=1e-6

  -odsn1, --dsn1-overlap-cutoff FLOAT

                                  Minimum cutoff necessary to keep prediction

                                  based on DSN1 significance if both other

                                  models meet their respective cutoffs.

                                  Between 0 and 1, default=0.5

  -odsn2, --dsn2-overlap-cutoff FLOAT

                                  Minimum cutoff necessary to keep prediction

                                  based on DSN2 significance if both other

                                  models meet their respective cutoffs.

                                  Between 0 and 1, default=0.5

  -ophmm, --phmm-overlap-cutoff INTEGER

                                  Minimum cutoff necessary to keep prediction

                                  based on pHMM significance if both other

                                  models meet their respective cutoffs.

                                  Between 0 and 1, default=1

  --help                          Show this message and exit.

 

 

実行方法

単離ゲノムのfastaファイルを指定する。

smorf single -t 12 -o outdir input_genome.fasta

 出力

f:id:kazumaxneo:20200807225425p:plain

small output.fa

f:id:kazumaxneo:20200807230148p:plain

メタゲノムのドラフトアセンブリを指定する。

smorf meta -t 12 -o outdir input_metagenome.fasta

 

引用

Automated prediction and annotation of small proteins in microbial genomes

Matthew G. Durrant, Ami S. Bhatt

bioRxiv, Posted July 28, 2020

 

関連