macでインフォマティクス

macでインフォマティクス

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

anvi'oを使ってメタゲノム解析を行う

2020 4/22  追記

2020 5/20 コード修正

 

 ハイスループットシーケンシングとオミックス技術の進歩は、自然界に存在する微生物群集の研究に革命をもたらしている。微生物のライフスタイルを包括的に調査するためには、遺伝情報を対話的に整理して可視化し、複雑なデータの分解能を高めるための微妙な違いを取り入れる能力が必要となる。ここでは、複数のソースからのオミックスデータを単一の直感的な表現にリンクすることができるインタラクティブなインターフェースを備えた、メタゲノムアセンブリ内の微生物ゲノムの自動化とヒト主導の特性評価を提供する先進的な解析および可視化プラットフォームであるanvi'oを紹介する。その拡張可能な可視化アプローチは、各コンティグに関する多次元の情報を抽出し、データの探索、操作、報告のためのダイナミックで統一された作業環境を提供する。Anvi'oを使用して、公開されているデータセットを再解析し、1塩基変異のデノボの特徴付けを通じて、微生物集団内のゲノムの時間的変化を探り、培養種やシングルセルゲノムをメタゲノムやメタトランスクリプトームデータとリンクさせた。Anvi'oは、広範なバイオインフォマティクスのスキルを持たない研究者でも、大規模な「オミックスデータセット」の詳細な分析を実行し、伝えることができるようにするオープンソースのプラットフォームである。

 

HP

http://merenlab.org/software/anvio/

 


 

Anvi'oを使用すると、メタゲノムのビニング、一塩基変異の分析、バクテリアパンゲノムの研究、メタゲノムアセンブリ内のバクテリアゲノム数予測、また真核生物アセンブリプロジェクトからの汚染の除去までを行うことができる。

積極的なバージョンアップによって機能も徐々に変わってきています。注意して使って下さい。

インストール

公式dockerイメージを使って複数のubuntu18.04LTSマシンでテストした。

本体 Github

 

#bioconda (link)注意;依存が多いため時間がかかる
conda create -n anvio -y
conda activate anvio
conda install -c bioconda anvio -y

#homebrew (not tested)
brew tap merenlab/anvio
brew install merenlab/anvio/anvio

テスト

> anvi-self-test --suite mini

 

#依存が多いので、condaだと依存チェックに異常な時間がかかる。動かすだけならdockerが楽。

#docker (dockerhub) (link)

#latest
docker pull meren/anvio:latest

#lauch
docker run --rm -it -v `pwd`:`pwd` -w `pwd` -p 8080:8080 meren/anvio:latest

> anvi-setup-ncbi-cogs -h

> anvi-setup-ncbi-cogs -h

usage: anvi-setup-ncbi-cogs [-h] [--cog-data-dir COG_DATA_DIR] [--reset]

                            [--just-do-it] [-T NUM_THREADS]

 

Download and setup NCBI's Clusters of Orthologous Groups database.

 

optional arguments:

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

  --cog-data-dir COG_DATA_DIR

                        The directory for COG data to be stored. If you leave

                        it as is without specifying anything, the default

                        destination for the data directory will be used to set

                        things up. The advantage of it is that everyone will

                        be using a single data directory, but then you may

                        need superuser privileges to do it. Using this

                        parameter you can choose the location of the data

                        directory somewhere you like. However, when it is time

                        to run COGs, you will need to remember that path and

                        provide it to the program.

  --reset               Remove all the previously stored files and start over.

                        If something is feels wrong for some reason and if you

                        believe re-downloading files and setting them up could

                        address the issue, this is the flag that will tell

                        anvi'o to act like a real computer scientist

                        challenged with a computational problem.

  --just-do-it          Don't bother me with questions or warnings, just do

                        it.

  -T NUM_THREADS, --num-threads NUM_THREADS

                        Maximum number of threads to use for multithreading

                        whenever possible. Very conservatively, the default is

                        1. It is a good idea to not exceed the number of CPUs

                        / cores on your system. Plus, please be careful with

                        this option if you are running your commands on a SGE

                        --if you are clusterizing your runs, and asking for

                        multiple threads to use, you may deplete your

                        resources very fast.

anvi-gen-contigs-database -h

> anvi-gen-contigs-database -h

usage: anvi-gen-contigs-database [-h] -f FASTA [-n PROJECT_NAME]

                                 [-o DB_FILE_PATH] [--description TEXT_FILE]

                                 [-L INT] [-K INT] [--skip-gene-calling]

                                 [--prodigal-translation-table INT]

                                 [--external-gene-calls GENE-CALLS]

                                 [--ignore-internal-stop-codons]

                                 [--skip-mindful-splitting]

 

Generate a new anvi'o contigs database.

 

optional arguments:

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

 

MANDATORY INPUTS:

  Things you really need to provide to be in business.

 

  -f FASTA, --contigs-fasta FASTA

                        The FASTA file that contains reference sequences you

                        mapped your samples against. This could be a reference

                        genome, or contigs from your assembler. Contig names

                        in this file must match to those in other input files.

                        If there is a problem anvi'o will gracefully complain

                        about it.

  -n PROJECT_NAME, --project-name PROJECT_NAME

                        Name of the project. Please choose a short but

                        descriptive name (so anvi'o can use it whenever she

                        needs to name an output file, or add a new table in a

                        database, or name her first born).

 

OPTIONAL INPUTS:

  Things you may want to tweak.

 

  -o DB_FILE_PATH, --output-db-path DB_FILE_PATH

                        Output file path for the new database.

  --description TEXT_FILE

                        A plain text file that contains some description about

                        the project. You can use Markdwon syntax. The

                        description text will be rendered and shown in all

                        relevant interfaces, including the anvi'o interactive

                        interface, or anvi'o summary outputs.

  -L INT, --split-length INT

                        Anvi'o splits very long contigs into smaller pieces,

                        without actually splitting them for real. These

                        'virtual' splits improve the efficacy of the

                        visualization step, and changing the split size gives

                        freedom to the user to adjust the resolution of their

                        display when necessary. The default value is (20000).

                        If you are planning to use your contigs database for

                        metagenomic binning, we advise you to not go below

                        10,000 (since the lower the split size is, the more

                        items to show in the display, and decreasing the split

                        size does not really help much to binning). But if you

                        are thinking about using this parameter for ad hoc

                        investigations other than binning, you should ignore

                        our advice, and set the split size as low as you want.

                        If you do not want your contigs to be split, you can

                        set the split size to '0' or any other negative

                        integer (lots of unnecessary freedom here, enjoy!).

  -K INT, --kmer-size INT

                        K-mer size for k-mer frequency calculations. The

                        default k-mer size for composition-based analyses is

                        4, historically. Although tetra-nucleotide frequencies

                        seem to offer the the sweet spot of sensitivity,

                        information density, and manageable number of

                        dimensions for clustering approaches, you are welcome

                        to experiment (but maybe you should leave it as is for

                        your first set of analyses).

  --skip-mindful-splitting

                        By default, anvi'o attempts to prevent soft-splitting

                        large contigs by cutting proper gene calls to make

                        sure a single gene is not broken into multiple splits.

                        This requires a careful examination of where genes

                        start and end, and to find best locations to split

                        contigs with respect to this information. So, when the

                        user asks for a split size of, say, 1,000, it serves

                        as a mere suggestion. When this flag is used, anvi'o

                        does what the user wants and creates splits at desired

                        lengths (although some functionality may become

                        unavailable for the projects that rely on a contigs

                        database that is initiated this way).

 

GENES IN CONTIGS:

  Expert thingies.

 

  --skip-gene-calling   By default, generating an anvi'o contigs database

                        includes the identification of open reading frames in

                        contigs by running a bacterial gene caller. Declaring

                        this flag will by-pass that process. If you prefer,

                        you can later import your own gene calling results

                        into the database.

  --prodigal-translation-table INT

                        This is a parameter to pass to the Prodigal for a

                        specific translation table. This parameter corresponds

                        to the parameter `-g` in Prodigal, the default value

                        of which is 11 (so if you do not set anything, it will

                        be set to 11 in Prodigal runtime. Please refer to the

                        Prodigal documentation to determine what is the right

                        translation table for you if you think you need it.)

  --external-gene-calls GENE-CALLS

                        A TAB-delimited file to utilize external gene calls.

                        The file must have these columns: 'gene_callers_id' (a

                        unique integer number for each gene call, start from

                        1), 'contig' (the contig name the gene call is found),

                        'start' (start position, integer), 'stop' (stop

                        position, integer), 'direction' (the direction of the

                        gene open reading frame; can be 'f' or 'r'), 'partial'

                        (whether it is a complete gene call, or a partial one;

                        must be 1 for partial calls, and 0 for complete

                        calls), 'source' (the gene caller), and 'version' (the

                        version of the gene caller, i.e., v2.6.7 or v1.0). An

                        example file can be found via the URL

                        https://bit.ly/2qEEHuQ

  --ignore-internal-stop-codons

                        This is only relevant when you have an external gene

                        calls file. If anvi'o figures out that your custom

                        gene calls result in amino acid sequences with stop

                        codons in the middle, it will complain about it. You

                        can use this flag to tell anvi'o to don't check for

                        internal stop codons, EVEN THOUGH IT MEANS THERE IS

                        MOST LIKELY SOMETHING WRONG WITH YOUR EXTERNAL GENE

                        CALLS FILE. Anvi'o will understand that sometimes we

                        don't want to care, and will not judge you. Instead,

                        it will replace every stop codon residue in the amino

                        acid sequence with an 'X' character. Please let us

                        know if you used this and things failed, so we can

                        tell you that you shouldn't have really used it if you

                        didn't like failures at the first place (smiley).

anvi-run-hmms -h

> anvi-run-hmms -h

usage: anvi-run-hmms [-h] -c CONTIGS_DB [-H HMM PROFILE PATH]

                     [-I HMM PROFILE NAME] [--also-scan-trnas]

                     [-T NUM_THREADS] [--just-do-it]

 

This program deals with populating tables that store HMM hits in an anvi'o

contigs database.

 

optional arguments:

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

 

DB:

  An anvi'o contigs adtabase to populate with HMM hits

 

  -c CONTIGS_DB, --contigs-db CONTIGS_DB

                        Anvi'o contigs database generated by 'anvi-gen-

                        contigs'

 

HMM OPTIONS:

  If you have your own HMMs, or if you would like to run only a set of

  default anvi'o HMM profiles rather than running them all, this is your

  stop.

 

  -H HMM PROFILE PATH, --hmm-profile-dir HMM PROFILE PATH

                        You can use this parameter you can specify a directory

                        path that contain an HMM profile. This way you can run

                        HMM profiles that are not included in anvi'o. See the

                        online to find out about the specifics of this

                        directory structure .

  -I HMM PROFILE NAME, --installed-hmm-profile HMM PROFILE NAME

                        When you run this program without any parameter, it

                        runs all 4 HMM profiles installed on your system. If

                        you want only a specific one to run, you can select it

                        by using this parameter. These are the currently

                        available ones: "Protista_83" (type: singlecopy),

                        "Bacteria_71" (type: singlecopy), "Archaea_76" (type:

                        singlecopy), "Ribosomal_RNAs" (type: Ribosomal_RNAs).

 

tRNAs:

  Through this program you can also scan Transfer RNA sequences in your

  contigs database for free (instead of running `anvi-scan-trnas` later).

 

  --also-scan-trnas     Also scan tRNAs while you're at it.

 

PERFORMANCE:

  Stuff everyone forgets to set and then get upset with how slow science

  goes.

 

  -T NUM_THREADS, --num-threads NUM_THREADS

                        Maximum number of threads to use for multithreading

                        whenever possible. Very conservatively, the default is

                        1. It is a good idea to not exceed the number of CPUs

                        / cores on your system. Plus, please be careful with

                        this option if you are running your commands on a SGE

                        --if you are clusterizing your runs, and asking for

                        multiple threads to use, you may deplete your

                        resources very fast.

 

AUTHORITY:

  Because you are the boss.

 

  --just-do-it          Don't bother me with questions or warnings, just do

                        it.

anvi-display-contigs-stats -h

> anvi-display-contigs-stats -h

usage: anvi-display-contigs-stats [-h] [--report-as-text] [-o FILE_PATH]                                                                                                                                                                                

                                  [-I IP_ADDR] [-P INT] [--browser-path PATH]

                                  [--server-only] [--password-protected]

                                  CONTIG DATABASES) [CONTIG DATABASE(S ...]

 

Start the anvi'o interactive interactive for viewing or comparing contigs

statistics

 

positional arguments:

  CONTIG DATABASE(S)    Anvio'o Contig databases to display statistics, you

                        can give multiple databases by seperating them with

                        space.

 

optional arguments:

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

 

REPORT CONFIGURATION:

  Specify what kind of output you want.

 

  --report-as-text      If you give this flag, Anvi'o will not open new

                        browser to show Contigs database statistics and write

                        all stats to TAB separated file and you should also

                        give --output-file with this flag otherwise Anvi'o

                        will complain.

  -o FILE_PATH, --output-file FILE_PATH

                        File path to store results.

 

SERVER CONFIGURATION:

  For power users.

 

  -I IP_ADDR, --ip-address IP_ADDR

                        IP address for the HTTP server. The default ip address

                        (0.0.0.0) should work just fine for most.

  -P INT, --port-number INT

                        Port number to use for anvi'o services. If nothing is

                        declared, anvi'o will try to find a suitable port

                        number, starting from the default port number, 8080.

  --browser-path PATH   By default, anvi'o will use your default browser to

                        launch the interactive interface. If you would like to

                        use something else than your system default, you can

                        provide a full path for an alternative browser using

                        this parameter, and hope for the best. For instance we

                        are using this parameter to call Google's experimental

                        browser, Canary, which performs better with demanding

                        visualizations.

  --server-only         The default behavior is to start the local server, and

                        fire up a browser that connects to the server. If you

                        have other plans, and want to start the server without

                        calling the browser, this is the flag you need.

  --password-protected  If this flag is set, command line tool will ask you to

                        enter a password and interactive interface will be

                        only accessible after entering same password. This

                        option is recommended for shared machines like

                        clusters or shared networks where computers are not

                        isolated.

anvi-run-ncbi-cogs -h

> anvi-run-ncbi-cogs -h

usage: anvi-run-ncbi-cogs [-h] -c CONTIGS_DB [--cog-data-dir COG_DATA_DIR]

                          [-T NUM_THREADS] [--sensitive]

                          [--temporary-dir-path PATH]

                          [--search-with SEARCH_METHOD]

 

Run NCBI's COGs to associate genes in an anvi'o contigs database with

functions. COGs database was been designed as an attempt to classify proteins

from completely sequenced genomes on the basis of the orthology concept. It is

no longer actively developed, however, it is still very effective for daily

needs. You may want to consider Pfams or the eggNOG database for more

comprehensive functional insights.

 

optional arguments:

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

  -c CONTIGS_DB, --contigs-db CONTIGS_DB

                        Anvi'o contigs database generated by 'anvi-gen-

                        contigs'

  --cog-data-dir COG_DATA_DIR

                        The directory path for your COG setup. Anvi'o will try

                        to use the default path if you do not specify

                        anything.

  -T NUM_THREADS, --num-threads NUM_THREADS

                        Maximum number of threads to use for multithreading

                        whenever possible. Very conservatively, the default is

                        1. It is a good idea to not exceed the number of CPUs

                        / cores on your system. Plus, please be careful with

                        this option if you are running your commands on a SGE

                        --if you are clusterizing your runs, and asking for

                        multiple threads to use, you may deplete your

                        resources very fast.

  --sensitive           DIAMOND sensitivity. With this flag you can instruct

                        DIAMOND to be 'sensitive', rather than 'fast' during

                        the search. It is likely the search will take

                        remarkably longer. But, hey, if you are doing it for

                        your final analysis, maybe it should take longer and

                        be more accurate. This flag is only relevant if you

                        are running DIAMOND.

  --temporary-dir-path PATH

                        If you don't provide anything here, this program will

                        come up with a temporary directory path by itself to

                        store intermediate files, and clean it later. If you

                        want to have full control over this, you can use this

                        flag to define one..

  --search-with SEARCH_METHOD

                        What program to use for database searching. The

                        default search uses diamond. All available options

                        include: diamond, blastp.

> anvi-get-sequences-for-gene-calls -h

> anvi-get-sequences-for-gene-calls -h

usage: anvi-get-sequences-for-gene-calls [-h] [-c CONTIGS_DB]

                                         [--gene-caller-ids GENE_CALLER_IDS]

                                         [--delimiter CHAR]

                                         [--report-extended-deflines]

                                         [--wrap WRAP] [--export-gff3]

                                         [--get-aa-sequences]

                                         [-g GENOMES_STORAGE]

                                         [-G GENOME_NAMES] -o FILE_PATH

 

A script to get back sequences for gene calls

 

optional arguments:

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

 

OPTION #1: EXPORT FROM CONTIGS DB:

  -c CONTIGS_DB, --contigs-db CONTIGS_DB

                        Anvi'o contigs database generated by 'anvi-gen-

                        contigs'

  --gene-caller-ids GENE_CALLER_IDS

                        Gene caller ids. Multiple of them can be declared

                        separated by a delimiter (the default is a comma). In

                        anvi-gen-variability-profile, if you declare nothing

                        you will get all genes matching your other filtering

                        criteria. In other programs, you may get everything,

                        nothing, or an error. It really depends on the

                        situation. Fortunately, mistakes are cheap, so it's

                        worth a try.

  --delimiter CHAR      The delimiter to parse multiple input terms. The

                        default is ','.

  --report-extended-deflines

                        When declared, the deflines in the resulting FASTA

                        file will contain more information.

  --wrap WRAP           When to wrap sequences when storing them in a FASTA

                        file. The default is '120'. A value of '0' would be

                        equivalent to 'do not wrap'.

  --export-gff3         If this is true, the output file will be in GFF3

                        format.

  --get-aa-sequences    Store amino acid sequences instead.

 

OPTION #2: EXPORT FROM A GENOMES STORAGE:

  -g GENOMES_STORAGE, --genomes-storage GENOMES_STORAGE

                        Anvi'o genomes storage file

  -G GENOME_NAMES, --genome-names GENOME_NAMES

                        Genome names to 'focus'. You can use this parameter to

                        limit the genomes included in your analysis. You can

                        provide these names as a comma-separated list of

                        names, or you can put them in a file, where you have a

                        single genome name in each line, and provide the file

                        path.

 

OPTIONS COMMON TO ALL INPUTS:

  -o FILE_PATH, --output-file FILE_PATH

                        File path to store results.

anvi-import-taxonomy-for-genes -h

> anvi-import-taxonomy-for-genes -h

usage: anvi-import-taxonomy-for-genes [-h] -c CONTIGS_DB [-p PARSER] -i FILES)

                                      [FILE(S ...] [--just-do-it]

 

Import gene-level taxonomy into an anvi'o contigs database.

 

optional arguments:

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

  -c CONTIGS_DB, --contigs-db CONTIGS_DB

                        Anvi'o contigs database generated by 'anvi-gen-

                        contigs'

  -p PARSER, --parser PARSER

                        Parser to make sense of the input files. There are 3

                        parsers readily available: ['default_matrix',

                        'centrifuge', 'kaiju']. It is OK if you do not select

                        a parser, but in that case there will be no additional

                        contigs available except the identification of single-

                        copy genes in your contigs for later use. Using a

                        parser will not prevent the analysis of single-copy

                        genes, but make anvio more powerful to help you make

                        sense of your results. Please see the documentation,

                        or get in touch with the developers if you have any

                        questions regarding parsers.

  -i FILE(S) [FILE(S) ...], --input-files FILE(S) [FILE(S) ...]

                        Input file(s) for selected parser. Each parser (except

                        "blank") requires input files to process that you

                        generate before running anvio. Please see the

                        documentation for details.

  --just-do-it          Don't bother me with questions or warnings, just do

                        it.

anvi-merge -h

> anvi-merge -h

usage: anvi-merge [-h] -c CONTIGS_DB [-o DIR_PATH] [-S NAME]

                  [--description TEXT_FILE] [--skip-hierarchical-clustering]

                  [--enforce-hierarchical-clustering]

                  [--distance DISTANCE_METRIC] [--linkage LINKAGE_METHOD] [-W]

                  SINGLE_PROFILES) [SINGLE_PROFILE(S ...]

 

Merge multiple anvio profiles

 

positional arguments:

  SINGLE_PROFILE(S)     Anvo'o single profiles to merge

 

optional arguments:

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

  -c CONTIGS_DB, --contigs-db CONTIGS_DB

                        Anvi'o contigs database generated by 'anvi-gen-

                        contigs'

  -o DIR_PATH, --output-dir DIR_PATH

                        Directory path for output files

  -S NAME, --sample-name NAME

                        It is important to set a sample name (using only ASCII

                        letters and digits and without spaces) that is unique

                        (considering all others). If you do not provide one,

                        anvi'o will try to make up one for you based on other

                        information, although, you should never let the

                        software to decide these things).

  --description TEXT_FILE

                        A plain text file that contains some description about

                        the project. You can use Markdwon syntax. The

                        description text will be rendered and shown in all

                        relevant interfaces, including the anvi'o interactive

                        interface, or anvi'o summary outputs.

  --skip-hierarchical-clustering

                        If you are not planning to use the interactive

                        interface (or if you have other means to add a tree of

                        contigs in the database) you may skip the step where

                        hierarchical clustering of your items are preformed

                        based on default clustering recipes matching to your

                        database type.

  --enforce-hierarchical-clustering

                        If you have more than 25,000 splits in your merged

                        profile, anvi-merge will automatically skip the

                        hierarchical clustering of splits (by setting --skip-

                        hierarchical-clustering flag on). This is due to the

                        fact that computational time required for hierarchical

                        clustering increases exponentially with the number of

                        items being clustered. Based on our experience we

                        decided that 25,000 splits is about the maximum we

                        should try. However, this is not a theoretical limit,

                        and you can overwrite this heuristic by using this

                        flag, which would tell anvi'o to attempt to cluster

                        splits regardless.

  --distance DISTANCE_METRIC

                        The distance metric for the hierarchical clustering.

                        If you do not use this flag, the default distance

                        metric will be used for each clustering configuration

                        which is "euclidean".

  --linkage LINKAGE_METHOD

                        The same story with the `--distance`, except, the

                        system default for this one is ward.

  -W, --overwrite-output-destinations

                        Overwrite if the output files and/or directories

                        exist.

anvi-profile -h

> anvi-profile -h

usage: anvi-profile [-h] [-i INPUT_BAM] [-c CONTIGS_DB] [--blank-profile]

                    [-o DIR_PATH] [-W] [-S NAME] [--report-variability-full]

                    [--skip-SNV-profiling] [--profile-SCVs]

                    [--description TEXT_FILE] [--cluster-contigs]

                    [--skip-hierarchical-clustering]

                    [--distance DISTANCE_METRIC] [--linkage LINKAGE_METHOD]

                    [-M INT] [--max-contig-length INT] [-X INT] [-V INT]

                    [--list-contigs] [--contigs-of-interest FILE]

                    [-T NUM_THREADS] [--queue-size INT]

                    [--write-buffer-size-per-thread INT] [--force-multi]

 

Creates a single anvi'o profile database. The default input to this program is

a BAM file. When it is run on a BAM file, depending on the user parameters,

the program quantifies coverage per nucleotide position (and averages them out

per contig), calculates single-nucleotide, single-codon, and single-amino acid

variants, and stores these data into appropriate tables. Anvi'o single

profiles can be merged by the program `anvi-merge`.

 

optional arguments:

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

 

INPUTS:

  There are two possible inputs for anvio profiler. You must to declare

  either of these two.

 

  -i INPUT_BAM, --input-file INPUT_BAM

                        Sorted and indexed BAM file to analyze. Takes a long

                        time depending on the length of the file and

                        parameters used for profiling.

  -c CONTIGS_DB, --contigs-db CONTIGS_DB

                        Anvi'o contigs database generated by 'anvi-gen-

                        contigs'

  --blank-profile       If you only have contig sequences, but no mapping data

                        (i.e., you found a genome and would like to take a

                        look from it), this flag will become very hand. After

                        creating a contigs database for your contigs, you can

                        create a blank anvi'o profile database to use anvi'o

                        interactive interface with that contigs database

                        without any mapping data.

 

EXTRAS:

  Things that are not mandatory, but can be useful if/when declared.

 

  -o DIR_PATH, --output-dir DIR_PATH

                        Directory path for output files

  -W, --overwrite-output-destinations

                        Overwrite if the output files and/or directories

                        exist.

  -S NAME, --sample-name NAME

                        It is important to set a sample name (using only ASCII

                        letters and digits and without spaces) that is unique

                        (considering all others). If you do not provide one,

                        anvi'o will try to make up one for you based on other

                        information, although, you should never let the

                        software to decide these things).

  --report-variability-full

                        One of the things anvi-profile does is to store

                        information about variable nucleotide positions.

                        Usually it does not report every variable position,

                        since not every variable position is genuine

                        variation. Say, if you have 1,000 coverage, and all

                        nucleotides at that position are Ts and only one of

                        them is a C, the confidence of that C being a real

                        variation is quite low. anvi'o has a simple algorithm

                        in place to reduce the impact of noise. However, using

                        this flag you can disable it and ask profiler to

                        report every single variation (which may result in

                        very large output files and millions of reports, but

                        you are the boss). Do not forget to take a look at '--

                        min-coverage-for-variability' parameter

  --skip-SNV-profiling  By default, anvi'o characterizes single-nucleotide

                        variation in each sample. The use of this flag will

                        instruct profiler to skip that step. Please remember

                        that parameters and flags must be identical between

                        different profiles using the same contigs database for

                        them to merge properly.

  --profile-SCVs        Anvi'o can perform accurate characterization of codon

                        frequencies in genes during profiling. While having

                        codon frequencies opens doors to powerful evolutionary

                        insights in downstream analyses, due to its

                        computational complexity, this feature comes 'off' by

                        default. Using this flag you can rise against the

                        authority, as you always should, and make anvi'o

                        profile codons.

  --description TEXT_FILE

                        A plain text file that contains some description about

                        the project. You can use Markdwon syntax. The

                        description text will be rendered and shown in all

                        relevant interfaces, including the anvi'o interactive

                        interface, or anvi'o summary outputs.

 

HIERARCHICAL CLUSTERING:

  Do you want your splits to be clustered? Yes? No? Maybe? Remember: By

  default, anvi-profile will not perform hierarchical clustering on your

  splits; but if you use `--blank` flag, it will try. You can skip that by

  using the `--skip-hierarchical-clustering` flag.

 

  --cluster-contigs     Single profiles are rarely used for genome binning or

                        visualization, and since clustering step increases the

                        profiling runtime for no good reason, the default

                        behavior is to not cluster contigs for individual

                        runs. However, if you are planning to do binning on

                        one sample, you must use this flag to tell anvi'o to

                        run cluster configurations for single runs on your

                        sample.

  --skip-hierarchical-clustering

                        If you are not planning to use the interactive

                        interface (or if you have other means to add a tree of

                        contigs in the database) you may skip the step where

                        hierarchical clustering of your items are preformed

                        based on default clustering recipes matching to your

                        database type.

  --distance DISTANCE_METRIC

                        The distance metric for the hierarchical clustering.

                        Only relevant if you are using `--cluster-contigs`

                        flag. The default is "euclidean".

  --linkage LINKAGE_METHOD

                        The linkage method for the hierarchical clustering.

                        Just like the distance metric this is only relevant if

                        you are using it with `--cluster-contigs` flag. The

                        default is "ward".

 

NUMBERS:

  Defaults of these parameters will impact your analysis. You can always

  come back to them and update your profiles, but it is important to make

  sure defaults are reasonable for your sample.

 

  -M INT, --min-contig-length INT

                        Minimum length of contigs in a BAM file to analyze.

                        The minimum length should be long enough for tetra-

                        nucleotide frequency analysis to be meaningful. There

                        is no way to define a golden number of minimum length

                        that would be applicable to genomes found in all

                        environments, but we chose the default to be 1000, and

                        have been happy with it. You are welcome to

                        experiment, but we advise to never go below 1,000. You

                        also should remember that the lower you go, the more

                        time it will take to analyze all contigs. You can use

                        --list-contigs parameter to have an idea how many

                        contigs would be discarded for a given M.

  --max-contig-length INT

                        Just like the minimum contig length parameter, but to

                        set a maximum. Basically this will remove any contig

                        longer than a certain value. Why would anyone need

                        this? Who knows. But if you ever do, it is here.

  -X INT, --min-mean-coverage INT

                        Minimum mean coverage for contigs to be kept in the

                        analysis. The default value is 0, which is for your

                        best interest if you are going to profile multiple BAM

                        files which are then going to be merged for a cross-

                        sectional or time series analysis. Do not change it if

                        you are not sure this is what you want to do.

  -V INT, --min-coverage-for-variability INT

                        Minimum coverage of a nucleotide position to be

                        subjected to SNV profiling. By default, anvi'o will

                        not attempt to make sense of variation in a given

                        nucleotide position if it is covered less than 10X.

                        You can change that minimum using this parameter.

 

CONTIGS:

  Sweet parameters of convenience

 

  --list-contigs        When declared, the program will list contigs in the

                        BAM file and exit gracefully without any further

                        analysis.

  --contigs-of-interest FILE

                        It is possible to analyze only a group of contigs from

                        a given BAM file. If you provide a text file, in which

                        every contig of interest is listed line by line, the

                        profiler would engine only on those contigs in the BAM

                        file and ignore the rest. This can be used for

                        debugging purposes, or to engine on a particular group

                        of contigs that were identified as relevant during the

                        interactive analysis.

 

PERFORMANCE:

  Performance settings for profiler

 

  -T NUM_THREADS, --num-threads NUM_THREADS

                        Maximum number of threads to use for multithreading

                        whenever possible. Very conservatively, the default is

                        1. It is a good idea to not exceed the number of CPUs

                        / cores on your system. Plus, please be careful with

                        this option if you are running your commands on a SGE

                        --if you are clusterizing your runs, and asking for

                        multiple threads to use, you may deplete your

                        resources very fast.

  --queue-size INT      The queue size for worker threads to store data to

                        communicate to the main thread. The default is set by

                        the class based on the number of threads. If you have

                        *any* hesitation about whether you know what you are

                        doing, you should not change this value.

  --write-buffer-size-per-thread INT

                        How many items should be kept in memory before they

                        are written do the disk. The default is 500 per

                        thread. So a single-threaded job would have a write

                        buffer size of 500, whereas a job with 4 threads would

                        have a write buffer size of 4*500. The larger the

                        buffer size, the less frequent the program will access

                        to the disk, yet the more memory will be consumed

                        since the processed items will be cleared off the

                        memory only after they are written to the disk. The

                        default buffer size will likely work for most cases.

                        Please keep an eye on the memory usage output to make

                        sure the memory use never exceeds the size of the

                        physical memory.

  --force-multi         This is not useful to non-developers. It forces the

                        multi-process routine even when 1 thread is chosen.

anvi-interactive -h

> anvi-interactive -h

usage: anvi-interactive [-h] [-p PROFILE_DB] [-c CONTIGS_DB]                                                                                                                                                   

                        [-C COLLECTION_NAME] [--manual-mode] [-f FASTA]

                        [-d VIEW_DATA] [-t NEWICK] [--items-order FLAT_FILE]

                        [-V ADDITIONAL_VIEW] [-A ADDITIONAL_LAYERS]

                        [--gene-mode] [--inseq-stats] [-b BIN_NAME]

                        [--view NAME] [--title NAME]

                        [--taxonomic-level {t_domain,t_phylum,t_class,t_order,t_family,t_genus,t_species}]

                        [--split-hmm-layers] [--hide-outlier-SNVs]

                        [--state-autoload NAME] [--collection-autoload NAME]

                        [--export-svg FILE_PATH] [--show-views]

                        [--skip-check-names] [-o DIR_PATH] [--dry-run]

                        [--show-states] [--list-collections]

                        [--skip-init-functions] [--skip-auto-ordering]

                        [--distance DISTANCE_METRIC]

                        [--linkage LINKAGE_METHOD] [-I IP_ADDR] [-P INT]

                        [--browser-path PATH] [--read-only] [--server-only]

                        [--password-protected] [--user-server-shutdown]

 

Start an anvi'o server for the interactive interface

 

optional arguments:

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

 

DEFAULT INPUTS:

  The interavtive interface can be started with and without anvi'o

  databases. The default use assumes you have your profile and contigs

  database, however, it is also possible to start the interface using ad hoc

  input files. See 'MANUAL INPUT' section for required parameters.

 

  -p PROFILE_DB, --profile-db PROFILE_DB

                        Anvi'o profile database

  -c CONTIGS_DB, --contigs-db CONTIGS_DB

                        Anvi'o contigs database generated by 'anvi-gen-

                        contigs'

  -C COLLECTION_NAME, --collection-name COLLECTION_NAME

                        If you have a collection in your profile database, you

                        can use this flag to start the interactive interface

                        with a tree showing your bins in your collection,

                        instead of each split. This is very useful when you

                        have imported your external binning results into

                        anvi'o, and want to see the distribution of your bins

                        across samples. In these cases anvi'o will cluster

                        your bins and based on multiple metrics. Because this

                        particular clustering will be done on the fly within

                        anvi'o interactive class, you get to define a

                        disntance metric and a linkage method using --linkage

                        and --distance parameters if you want!

 

MANUAL INPUTS:

  Mandatory input parameters to start the interactive interface without

  anvi'o databases.

 

  --manual-mode         Using this flag, you can run the interactive interface

                        in an ad hoc manner using input files you curated

                        instead of standard output files generated by an

                        anvi'o run. In the manual mode you will be asked to

                        provide a profile database. In this mode a profile

                        database is only used to store 'state' of the

                        interactive interface so you can reload your visual

                        settings when you re-analyze the same files again. If

                        the profile database you provide does not exist,

                        anvi'o will create an empty one for you.

  -f FASTA, --fasta-file FASTA

                        A FASTA-formatted input file

  -d VIEW_DATA, --view-data VIEW_DATA

                        A TAB-delimited file for view data

  -t NEWICK, --tree NEWICK

                        NEWICK formatted tree structure

  --items-order FLAT_FILE

                        A flat file that contains the order of items you wish

                        the display using the interactive interface. You may

                        want to use this if you have a specific order of items

                        in your mind, and do not want to display a tree in the

                        middle (or simply you don't have one). The file format

                        is simple: each line should have an item name, and

                        there should be no header.

 

ADDITIONAL STUFF:

  Parameters to provide additional layers, views, or layer data.

 

  -V ADDITIONAL_VIEW, --additional-view ADDITIONAL_VIEW

                        A TAB-delimited file for an additional view to be used

                        in the interface. This file should contain all split

                        names, and values for each of them in all samples.

                        Each column in this file must correspond to a sample

                        name. Content of this file will be called 'user_view',

                        which will be available as a new item in the 'views'

                        combo box in the interface

  -A ADDITIONAL_LAYERS, --additional-layers ADDITIONAL_LAYERS

                        A TAB-delimited file for additional layers for splits.

                        The first column of this file must be split names, and

                        the remaining columns should be unique attributes. The

                        file does not need to contain all split names, or

                        values for each split in every column. Anvi'o will try

                        to deal with missing data nicely. Each column in this

                        file will be visualized as a new layer in the tree.

 

GENE MODE:

  Gene mode related parameters.

 

  --gene-mode           Initiate the interactive interface in 'gene mode'. In

                        this mode, the items are genes (instead of splits of

                        contigs). The following views are available: detection

                        (the detection value of each gene in each sample). The

                        mean_coverage (the mean coverage of genes). The

                        non_outlier_mean_coverage (the mean coverage of the

                        non-outlier nucleotide positions of each gene in each

                        sample (median absolute deviation is used to remove

                        outliers per gene per sample)). The

                        non_outlier_coverage_std view (standard deviation of

                        the coverage of non-outlier positions of genes in

                        samples). You can also choose to order items and

                        layers according to each one of the aforementioned

                        views. In addition, all layer ordering that are

                        available in the regular mode (i.e. the full mode

                        where you have contigs/splits) are also available in

                        'gene mode', so that, for example, you can choose to

                        order the layers according to 'detection', and that

                        would be the order according to the detection values

                        of splits, whereas if you choose 'genes_detections'

                        then the order of layers would be according to the

                        detection values of genes. Inspection and sequence

                        functionality are available (through the right-click

                        menu), except now sequences are of the specific gene.

                        Inspection has now two options available: 'Inspect

                        Context', which brings you to the inspection page of

                        the split to which the gene belongs where the

                        inspected gene will be highlighted in yellow in the

                        bottom, and 'Inspect Gene', which opens the inspection

                        page only for the gene and 100 nts around each side of

                        it (the purpose of this option is to make the

                        inspection page load faster if you only want to look

                        at the nucleotide coverage of a specific gene).

                        NOTICE: You can't store states or collections in 'gene

                        mode'. However, you still can make fake selections,

                        and create fake bins for your viewing convenience only

                        (smiley). Search options are available, and you can

                        even search for functions if you have them in your

                        contigs database. ANOTHER NOTICE: loading this mode

                        might take a while if your bin has many genes, and

                        your profile database has many samples, this is

                        because the gene coverages stats are computed in an

                        ad-hoc manner when you load this mode, we know this is

                        not ideal and we plan to improve that (along with

                        other things). If you have suggestions/complaints

                        regarding this mode please comment on this github

                        issue: https://goo.gl/yHhRei. Please refer to the

                        online tutorial for more information.

  --inseq-stats         Provide if working with INSeq/Tn-Seq genomic data.

                        With this, all gene level coverage stats will be

                        calculated using INSeq/Tn-Seq statistical methods.

  -b BIN_NAME, --bin-id BIN_NAME

                        Bin name you are interested in.

 

VISUALS RELATED:

  Parameters that give access to various adjustements regarding the

  interface.

 

  --view NAME           Start the interface with a pre-selected view. To see a

                        list of available views, use --show-views flag.

  --title NAME          Title for the interface. If you are working with a

                        RUNINFO dict, the title will be determined based on

                        information stored in that file. Regardless, you can

                        override that value using this parameter.

  --taxonomic-level {t_domain,t_phylum,t_class,t_order,t_family,t_genus,t_species}

                        The taxonomic level to use whenever relevant and/or

                        available. The default taxonomic level is t_genus, but

                        if you choose something specific, anvi'o will focus on

                        that whenever possible.

  --split-hmm-layers    When declared, this flag tells the interface to split

                        every gene found in HMM searches that were performed

                        against non-singlecopy gene HMM profiles into their

                        own layer. Please see the documentation for details.

  --hide-outlier-SNVs   During profiling, anvi'o marks positions of single-

                        nucleotide variations (SNVs) that originate from

                        places in contigs where coverage values are a bit

                        'sketchy'. If you would like to avoid SNVs in those

                        positions of splits in applicable projects you can use

                        this flag, and the interface would hide SNVs that are

                        marked as 'outlier' (although it is clearly the best

                        to see everything, no one will judge you if you end up

                        using this flag) (plus, there may or may not be some

                        historical data on this here:

                        https://github.com/meren/anvio/issues/309).

  --state-autoload NAME

                        Automatically load previous saved state and draw tree.

                        To see a list of available states, use --show-states

                        flag.

  --collection-autoload NAME

                        Automatically load a collection and draw tree. To see

                        a list of available collections, use --list-

                        collections flag.

  --export-svg FILE_PATH

                        The SVG output file path.

 

SWEET PARAMS OF CONVENIENCE:

  Parameters and flags that are not quite essential (but nice to have).

 

  --show-views          When declared, the program will show a list of

                        available views, and exit.

  --skip-check-names    For debugging purposes. You should never really need

                        it.

  -o DIR_PATH, --output-dir DIR_PATH

                        Directory path for output files

  --dry-run             Don't do anything real. Test everything, and stop

                        right before wherever the developer said 'well, this

                        is enough testing', and decided to print out results.

  --show-states         When declared the program will print all available

                        states and exit.

  --list-collections    Show available collections and exit.

  --skip-init-functions

                        When declared, function calls for genes will not be

                        initialized (therefore will be missing from all

                        relevant interfaces or output files). The use of this

                        flag may reduce the memory fingerprint and processing

                        time for large datasets.

  --skip-auto-ordering  When declared, the attempt to include automatically

                        generated orders of items based on additional data is

                        skipped. In case those buggers cause issues with your

                        data, and you still want to see your stuff and deal

                        with the other issue maybe later.

  --distance DISTANCE_METRIC

                        The distance metric for the hierarchical clustering.

                        Only relevant if you are running the interactive

                        interface in "collection" mode. The default is

                        "euclidean".

  --linkage LINKAGE_METHOD

                        The linkage method for the hierarchical clustering.

                        Only relevant if you are running the interactive

                        interface in "collection" mode. The default is "ward".

 

SERVER CONFIGURATION:

  For power users.

 

  -I IP_ADDR, --ip-address IP_ADDR

                        IP address for the HTTP server. The default ip address

                        (0.0.0.0) should work just fine for most.

  -P INT, --port-number INT

                        Port number to use for anvi'o services. If nothing is

                        declared, anvi'o will try to find a suitable port

                        number, starting from the default port number, 8080.

  --browser-path PATH   By default, anvi'o will use your default browser to

                        launch the interactive interface. If you would like to

                        use something else than your system default, you can

                        provide a full path for an alternative browser using

                        this parameter, and hope for the best. For instance we

                        are using this parameter to call Google's experimental

                        browser, Canary, which performs better with demanding

                        visualizations.

  --read-only           When the interactive interface is started with this

                        flag, all 'database write' operations will be

                        disabled.

  --server-only         The default behavior is to start the local server, and

                        fire up a browser that connects to the server. If you

                        have other plans, and want to start the server without

                        calling the browser, this is the flag you need.

  --password-protected  If this flag is set, command line tool will ask you to

                        enter a password and interactive interface will be

                        only accessible after entering same password. This

                        option is recommended for shared machines like

                        clusters or shared networks where computers are not

                        isolated.

  --user-server-shutdown

                        Allow users to shutdown an anvi'server via web

                        interface.

 > anvi-script-reformat-fasta -h

> anvi-script-reformat-fasta -h

usage: anvi-script-reformat-fasta [-h] [-l MIN_LENGTH]

                                  [--max-percentage-gaps PERCENTAGE]

                                  [-i TXT FILE] [-I TXT FILE] -o FASTA FILE

                                  [--simplify-names] [--prefix PREFIX]

                                  [-r REPORT FILE]

                                  FASTA FILE

 

Reformat FASTA file (remove contigs based on length, or based on a given list

of deflines, and/or generate an output with simpler names)

 

positional arguments:

  FASTA FILE

 

optional arguments:

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

  -l MIN_LENGTH, --min-len MIN_LENGTH

                        Minimum length of contigs to keep (contigs shorter

                        than this value will not be included in the output

                        file). The default is 0, so nothing will be removed if

                        you do not declare a minimum size.

  --max-percentage-gaps PERCENTAGE

                        Maximum fraction of gaps in a sequence (any sequence

                        with more gaps will be removed from the output FASTA

                        file). The default is 100.000000.

  -i TXT FILE, --exclude-ids TXT FILE

                        IDs to remove from the FASTA file. You cannot provide

                        both --keep-ids and --exclude-ids.

  -I TXT FILE, --keep-ids TXT FILE

                        If provided, all IDs not in this file will be excluded

                        from the reformatted FASTA file. Any additional

                        filters (such as --min-len) will still be applied to

                        the IDs in this file. You cannot provide both

                        --exclude-ids and --keep-ids.

  -o FASTA FILE, --output-file FASTA FILE

                        Output file path.

  --simplify-names      Edit deflines to make sure they contigs have simple

                        names.

  --prefix PREFIX       Use this parameter if you would like to add a prefix

                        to your contig names while simplifying them. The

                        prefix must be a single word (you can use underscor

                        character, but nothing more!).

  -r REPORT FILE, --report-file REPORT FILE

                        Report file path. When you run this program with

                        `--simplify-names` flag, all changes to deflines will

                        be reported in this file in case you need to go back

                        to this information later. It is not mandatory to

                        declare one, but it is a very good idea to have it.

anvi-export-splits-and-coverages -h

> anvi-export-splits-and-coverages -h

usage: anvi-export-splits-and-coverages [-h] -p PROFILE_DB -c CONTIGS_DB

                                        [-o DIR_PATH] [-O FILENAME_PREFIX]

                                        [--splits-mode] [--report-contigs]

                                        [--use-Q2Q3-coverages]

 

Export split or contig sequences and coverages across samples stored in an

anvi'o profile database. This program is especially useful if you would like

to 'bin' your splits or contigs outside of anvi'o and import the binning

results into anvi'o using `anvi-import-collection` program.

 

optional arguments:

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

  -p PROFILE_DB, --profile-db PROFILE_DB

                        Anvi'o profile database

  -c CONTIGS_DB, --contigs-db CONTIGS_DB

                        Anvi'o contigs database generated by 'anvi-gen-

                        contigs'

  -o DIR_PATH, --output-dir DIR_PATH

                        Directory path for output files

  -O FILENAME_PREFIX, --output-file-prefix FILENAME_PREFIX

                        A prefix to be used while naming the output files (no

                        file type extensions please; just a prefix).

  --splits-mode         Specify this flag if you would like to output

                        coverages of individual 'splits', rather than their

                        'parent' contig coverages.

  --report-contigs      By default this program reports sequences and their

                        coverages for 'splits'. By using this flag, you can

                        report contig sequences and coverages instead. For

                        obvious reasons, you can't use this flag with

                        `--splits-mode` flag.

  --use-Q2Q3-coverages  By default this program reports the mean coverage of a

                        split (or contig, see --report-contigs) for each

                        sample. By using this flag, you can report the mean

                        Q2Q3 coverage by excluding 25 percent of the

                        nucleotide positions with the smallest coverage

                        values, and 25 percent of the nucleotide positions

                        with the largest coverage values. The hope is that

                        this removes 'outlier' positions resulting from non-

                        specific mapping, etc. that skew the mean coverage

                        estimate.

> anvi-import-collection -h

> anvi-import-collection -h

usage: anvi-import-collection [-h] [-c CONTIGS_DB] [-p PAN_OR_PROFILE_DB] -C

                              COLLECTION_NAME [--bins-info BINS_INFO]

                              [--contigs-mode]

                              TAB DELIMITED FILE

 

Import an external binning result into anvi'o

 

positional arguments:

  TAB DELIMITED FILE    The input file that describes bin IDs for each split

                        or contig.

 

optional arguments:

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

  -c CONTIGS_DB, --contigs-db CONTIGS_DB

                        Anvi'o contigs database generated by 'anvi-gen-

                        contigs'

  -p PAN_OR_PROFILE_DB, --pan-or-profile-db PAN_OR_PROFILE_DB

                        Anvi'o pan or profile database (and even genes

                        database in appropriate contexts).

  -C COLLECTION_NAME, --collection-name COLLECTION_NAME

                        Collection name.

  --bins-info BINS_INFO

                        Additional information for bins. The file must contain

                        three TAB-delimited columns, where the first one must

                        be a unique bin name, the second should be a 'source',

                        and the last one should be a 7 character HTML color

                        code (i.e., '#424242'). Source column must contain

                        information about the origin of the bin. If these bins

                        are automatically identified by a program like

                        CONCOCT, this column could contain the program name

                        and version. The source information will be associated

                        with the bin in various interfaces so in a sense it is

                        not *that* critical what it says there, but on the

                        other hand it is, becuse we should also think about

                        people who may end up having to work with what we put

                        together later.

  --contigs-mode        Use this flag if your binning was done on contigs

                        instead of splits. Please refer to the documentation

                        for help.

 

 

 

実行方法

の手順に則り進める。

 

 1、NCBI COG(Clusters of Orthologus Groups) データベースの準備。こちらは初回のみ実行する。

anvi-setup-ncbi-cogs -T 40

dockerイメージを使っている場合、一度実行してcommitする。次回以降はそれを使えば手間が減る。

 

2、contig.fastaの準備とデータベース作成

2、 分析対象のメタゲノムのcontigg.fastaやbinned.fastaを準備する。配列が多すぎると階層的クラスタリングの時にエラーになるので注意する。回避するには短い配列を捨てて配列数を減らす。

 

補足ステップ。======================================

FASTAファイルのヘッダーのdeflinesを修正(option)する。また、サイズ選択も可能。。スペースなどあるとステップ2でエラーを起こす。

anvi-script-reformat-fasta -l 5000 -o contigs.fa input_contigs.fa

contigs.faが出力される。修正されなかった場合、ヘッダーをシンプルな名前に置換する。あとで使うアラインメントのbamファイルは、修正後のfastaを使って作っていないとエラーになる。

ヘッダやファイル名で割と一般的に使われるのが"-", "<space>", "-"などだが、これらはファイル内にあってもファイル名にあってもエラーを起こす。必ず置換しておく。アンダーバー”_”にしておけばエラーは起きない。

================================================ 

 

binned.fastaからデータベースを作成する。コンティグに関連する情報(ORFの位置、各コンティグのk-mer頻度、スプリットの開始位置と終了位置、Prodigalを使った遺伝子の機能的および分類学アノテーションなど)のデータベースとなる。

anvi-gen-contigs-database -f binned.fasta -o contigs.db -n 'An example contigs1 datbase'

複数あるなら順番に作成
anvi-gen-contigs-database -f binned2.fasta -o contigs2.db -n 'An example contigs2 datbase'
anvi-gen-contigs-database -f binned3.fasta -o contigs3.db -n 'An example contigs3 datbase'

 

3、コンティグデータベースを、プラットフォームに同梱されている HMM モデル(現時点では、複数のバクテリアのシングルコピー遺伝子コレクションが公開されている)からのヒットでデコレートする。できるだけ多くスレッドを当てる。dbが複数あるなら全て行う。dbを統合するまで以後も同じ。

anvi-run-hmms -c contigs.db -T 40
  •  -T    Maximum number of threads to use for multithreading whenever possible. Very conservatively, the default is 1. It is a good idea to not exceed the number of CPUs/ cores on your system. Plus, please be careful with this option if you are running your commands on a SGE --if you are clusterizing your runs, and asking for multiple threads to use, you may deplete your resources very fast.

4、コンティグの統計を表示(コンティグスデータベースとHMMモデルが既に作成されていること)

anvi-display-contigs-stats contigs.db

http://127.0.0.1:8080にアクセスしてstatisticsを確認する。

f:id:kazumaxneo:20200402094703p:plain

f:id:kazumaxneo:20200402094706p:plain

確認し終わったら"Ctrl + C"で停止。

 

5、NCBI COGを使ってコンティグスデータベースの遺伝子をアノテーションするためのプログラムanvi-run-ncbi-cogsを実行する。DIAMONDが動くので、できるだけ多くスレッドを当てる。

anvi-run-ncbi-cogs -c contigs.db --num-threads 40

 

6、NCBI COGを使ってコンティグスデータベースの遺伝子をアノテーションするためのプログラムanvi-run-ncbi-cogsを実行する。

anvi-get-sequences-for-gene-calls -c contigs.db -o gene-calls.fa

 

================================================

追加ステップ - centrifugiのランとtaxonomyのインポート(参考

各遺伝子のtaxonomyアノテーションを持っていて、それをデータベースに入れてキュレーションしたい時に実行する。kaijuやcentrifugeが使えるが、centrifugeだと以下のようにする。ステップ6のgene-calls.faを使い、各遺伝子へのtaxonomyアノテーションを行う。centrifugeはdockerイメージにも最初からインストールされている。データベースだけ用意すればよい。

# centrifugeプリビルドデータベース(初回のみ)
wget ftp://ftp.ccb.jhu.edu/pub/infphilo/centrifuge/data/p_compressed+h+v.tar.gz
tar -zxvf p_compressed+h+v.tar.gz
#p_compressed+h+v1.cf、p_compressed+h+v2.cf、p_compressed+h+v3.cfができる。ラン時は"-x p_compressed+h+v"と指定する。

#centrifugeのrun
centrifuge -f -x p_compressed+h+v gene-calls.fa -S centrifuge_hits.tsv -p 40
#=> centrifuge_report.tsvとcentrifuge_hits.tsvができる。

#anvi'oにcentrifugeの結果を取り込む。
anvi-import-taxonomy-for-genes -c contigs.db -i centrifuge_report.tsv centrifuge_hits.tsv -p centrifuge

エラーなくランできていれば、視覚化の際にtaxonomyのオプションが利用できるようになる。

================================================

 

7、sortしたbamファイルとbam.baiの準備。minimap2を使うなら、

#sample1
minimap2 -R "@RG\tID:X\tLB:Y\tSM:sample1\tPL:ILLUMINA" -t 40 -ax sr \
contigs.fa sample1_R1.fq.gz sample1_R2.fq.gz \
|samtools sort -@ 40 -m 2G -O BAM - > sample1.bam \
&& samtools index -@ 8 sample1.bam

#sample2
minimap2 -R "@RG\tID:X\tLB:Y\tSM:sample2\tPL:ILLUMINA" -t 40 -ax sr \
contigs.fa sample2_R1.fq.gz sample2_R2.fq.gz \
|samtools sort -@ 40 -m 2G -O BAM - > sample2.bam \
&& samtools index -@ 8 sample2.bam

#sample3
minimap2 -R "@RG\tID:X\tLB:Y\tSM:sample3\tPL:ILLUMINA" -t 40 -ax sr \
contigs.fa sample3_R1.fq.gz sample3_R2.fq.gz \
|samtools sort -@ 40 -m 2G -O BAM - > sample3.bam \
&& samtools index -@ 8 sample3.bam

 

8、サンプルに関する固有の情報のデータベースの作成。サンプルのbamごとに実施する。デフォルトでは 2,500-bp より長いコンティグを処理する( --min-contig-length)。

#sample1
anvi-profile -i sample1.bam -c contigs.db --output-dir PROFILES/SAMPLE1_Profile --sample-name sample1 -T 50

#sample2
anvi-profile -i sample2.bam -c contigs.db --output-dir PROFILES/SAMPLE2_Profile --sample-name sample2 -T 50

#sample3
anvi-profile -i sample3.bam -c contigs.db --output-dir PROFILES/SAMPLE3_Profile --sample-name sample3 -T 50

もし1サンプル(one bam file)しかないなら"--cluster-contigs"をつけてランする。

ステップ2ではコンティグのデータベースcontig.dbができたが、このステップではサンプルごとのデータベースであるPROFILE.dbができる。

f:id:kazumaxneo:20200403152637p:plain

うまくいけばHappyが出ます。

 

9、8の結果をマージし、クラスタリングを実行する。1サンプルのみの場合は不要。

anvi-merge \
PROFILES/SAMPLE1_Profile/PROFILE.db \
PROFILES/SAMPLE2_Profile/PROFILE.db \
PROFILES/SAMPLE3_Profile/PROFILE.db \
-o SAMPLES-MERGED -c contigs.db --enforce-hierarchical-clustering
  • --skip-hierarchical-clustering  If you are not planning to use the interactive interface (or if you have other means to add a tree of contigs in the database) you may skip the step where hierarchical clustering of your items are preformed based on default clustering recipes matching to your database type.
  • --enforce-hierarchical-clustering    If you have more than 25,000 splits in your merged  profile, anvi-merge will automatically skip the hierarchical clustering of splits (by setting --skip-hierarchical-clustering flag on). This is due to the  fact that computational time required for hierarchical clustering increases exponentially with the number of  items being clustered. Based on our experience we  decided that 25,000 splits is about the maximum we should try. However, this is not a theoretical limit,   and you can overwrite this heuristic by using this flag, which would tell anvi'o to attempt to cluster splits regardless.

versioon6+以降、ビニングは別コマンド anvi-cluster-contigsになり、このマージでは実行されなくなっている。この処理で階層的クラスタリングが実行されるが、--enforce-hierarchical-clustering をつけていてもクラスタリングに失敗することがある。失敗すると視覚化できないので、短い配列を減らすなどしてやり直す。

 

 

========================================================

追加ステップ1;コンティグを独自にビニングしている場合は、その結果をコレクションとしてマージされたプロファイルデータベースにインポートできる。

anvi-import-collection binning_results.txt -p SAMPLES-MERGED/PROFILE.db -c contigs.db --source "SOURCE_NAME"

binning_results.txtはTAB区切りのテキストファイルで、どのコンティグがどのビンに属しているかという情報を含んでいる必要がある。具体的には、コンティグ名<TAB>それが属するビン名、というタブ区切りファイルを用意する()。

 

追加ステップ2;カバレッジ情報と配列構成情報のエクスポート

anvi-export-splits-and-coverages -p SAMPLES-MERGED/PROFILE.db -c contigs.db -o output

========================================================

 

10、視覚化と分析。

anvi-interactive -p SAMPLES-MERGED/PROFILE.db -c contigs.db --server-only -P 8080

anvi-interactive を実行すると、各ビンの様々な特性、すなわち平均カバレッジ、分散、補完性、冗長性の推定値をその場で計算し、インタラクティブインターフェイスでサンプル全体の分布を表示する。

 

http://localhost:8080にアクセスして結果を確認する。 

左下のDrawボタンを押すと描画される。

f:id:kazumaxneo:20200403083928p:plain

単離サンプルのゲノム1つだけを分析した。複数のコンティグをクラスタリングし、環状のレイヤーにカバレッジGCなどを表現している。中央のデンドログラムがコンティグのクラスタリング結果を表現している。1つの環状ゲノムを表現している訳ではなく、全コンティグをカバレッジクラスタリングして環状に並べているということ (右上1/4が注釈になっていても何も問題はない)。

 

環状表現からlinear phylogramに変更した。

f:id:kazumaxneo:20200403183014p:plain

 

メタゲノムアセンブリ(37,297 sequecnes、total 488,010,807-bp) 

f:id:kazumaxneo:20200405102017p:plain

10 sample表示してます。

 

拡大した。カバレッジに基づいてクラスタリングされていることが分かる。

f:id:kazumaxneo:20200405102139p:plain

中心部のデンドログラムの枝部分にマウスホバーすると、特定のクレードだけハイライト表示される。

 

ハイライト表示してターゲットを定め、左クリックすることで名前をアサインできる。追加される名前はBinタブで制御する。

f:id:kazumaxneo:20200405103328p:plain

追加した後でも、Binタブの名前を変えれば変更可能。

 

BinタブでNew Binをクリックして名前を追加し、

f:id:kazumaxneo:20200405103551p:plain

 

Bin1とBin2を追加した。

f:id:kazumaxneo:20200405103529p:plain

 

興味あるクラスターが見つかったら、マウスホイールで拡大して右クリックする。

f:id:kazumaxneo:20200405102450p:plain

 

表示されるメニューの中のinspect splitからコンティグの全サンプルカバレッジを確認可能。長いコンティグはすべてsplitして扱われるため、splitt_xxxという表記になっている。

f:id:kazumaxneo:20200405103854p:plain

緑の線はGC%。下には予測されORFがオーバーレイ表示されている。ポジションは右上から変更可能。

 

右クリックからsplit配列のnrやrefseqへのblastingもできる。

f:id:kazumaxneo:20200405102454p:plain

 

ステップ6でtaxonomyをインポートしていると、taxonomyのレイヤーが最外周に追加されている。

f:id:kazumaxneo:20200405104844p:plain

 

Legendタブで色を修正可能。

f:id:kazumaxneo:20200405104635p:plain

 

LayerタブでトータルリードとSNVの色を変えた。

f:id:kazumaxneo:20200405105132p:plain

 

左のメニューから、カバレッジの表示方法について選択できる(詳細)。

f:id:kazumaxneo:20200405105746p:plain



 

abundance treeにした。カバレッジabundanceに基づいてサンプルの順番(円の中心から外周のリングの順番)が並べ替えられ、右上にそのクラスタリングのデンドログラムが追加された。

f:id:kazumaxneo:20200405105711p:plain

 

15サンプル

f:id:kazumaxneo:20200405114519p:plain

 

補足

重すぎる時は画面をリロードして下さい。ただし保存されてない設定は初期化されます。

 

MainのShow Advanced settingsを展開すると、図のサイズなど変更できる。

f:id:kazumaxneo:20200403183315p:plain

DENDROGRAMのアングル(0-270度)も変えることができる。

 

右下のボタンから出力すると、レジェンドつきでダウンロードディレクトリに.SVGが出力される。

f:id:kazumaxneo:20200422095023p:plain

引用

Anvi'o: an advanced analysis and visualization platform for 'omics data

Eren AM, Esen ÖC, Quince C, Vineis JH, Morrison HG, Sogin ML, Delmont TO

PeerJ. 2015 Oct 8;3:e1319

 

参考


関連


 

 

mobile element を検出する Mobster

 

 転移因子(ME)は自律的にコピーしたりゲノム上を移動したりすることができるDNA配列だが、その高度に反復的な配列構造のために検出が困難である。MEは、ゲノム構造を変化させる主要な進化ドライバーであるだけでなく、機能的に重要な領域に挿入され、遺伝子の機能を破壊することで、多くのヒト疾患の病原性変異を直接的にもたらしてきた。MEは、転移の様式によって2つの異なるクラスに分類される。クラスIのレトロトランスポゾンはコピー&ペーストによってRNA中間体を経由して移動し、クラスIIのDNAトランスポゾンはDNA中間体を持ち、一般的にはカット&ペーストによって移動する。これらのエレメントを合わせてヒトゲノムの大部分を構成しており、ヒトゲノム配列の45%から69%がこれらのトランスポゾンクラスのいずれかに属すると推定されている[ref.3, 4]。

 現在のところ、ヒトゲノム上で活性な、あるいは「ホット」なMEはわずかしか存在していないが、そのすべてがレトロトランスポゾンクラスに属し、自律型のL1ファミリー(6kb、50万コピー)、非自律型のAluファミリー(300bp、100万コピー)、SVAファミリー(2kb、3000コピー)などが含まれている。これらのMEファミリーは、DNAの新しい領域への挿入、DNAの形質転換、エクソンのシャッフル、加工された疑似遺伝子の生成などにより、ゲノム構造を変化させ続けている。転座の古代の遺物であっても、その配列の相同性が不均等なクロスオーバーを引き起こし、2つのMEコピー間でDNAの欠失や重複を引き起こす可能性があるため、ゲノムの変異に大きく貢献している。

 ME転移は、生殖細胞内または初期胚発生期にしばしば起こる。ヒトにおける最初のME挿入(MEI)は、血友病Aの2人の患者において、FVIII遺伝子のエクソン14に発見された[ref.10]。それ以来、90以上のMEIが発見されており、そのうち60個のAluエレメントの挿入、25個のL1sの挿入、7個のSVAの挿入が含まれている[ref.8]。さらに、MEはがんの発生に関与していることが知られており、いくつかの研究[ref.11-13]で腫瘍特異的なMEIイベントが発見されている。

 多型MEI(pMEI)を同定するために、ターゲットおよび次世代シークエンシング(NGS)解析の両方が開発されてきた。ヒトのNGSデータ中のpMEIを計算で検出しようとするこれまでの試みは、一般的に、pMEIを同定するために、不一致リードペアまたはクリップされたリードを使用する。HormozdiariらはVariationHunterを改良して多型のAlu挿入を特徴づけるようにした[ref.14]。一方、EwingとKazazianは多型のL1挿入を検出するためのパイプラインを開発した[ref.15]。Tea [ref.13]とRetroSeq [ref.16]は、MEIイベントのブレークポイントを微調整するために、不一致ペアに加えてクリップされたリードを使用することができる。最後に、Stewartらによる未発表のパイプラインでは、ペアエンドのIlluminaデータのペアエンドアプローチ[ref.17]に加えて、より長い454リードのpMEIを検出するためスプリットリード法を使用する。

 我々(本著者ら)は、WGSとWESの両方のデータにおいて、高精度でアクティブな非リファレンスMEIを検出することができるMobsterと名付けられた新しい方法を提示する。さらに、この手法は特定のMEIイベントのファミリーに限定されることなく、すべてのアクティブなMEIイベントのファミリーを検出することができる。この手法は、公開されているヒトのデータセットや、カバレッジの異なるシミュレーションデータにおいて、既存のツールを凌駕している。次に、ペアエンドWGSデータセット、ペアエンドWESデータセット、シングルエンドWESデータセットを含む様々なNGSデータタイプにMobsterを適用し、PCRバリデーションを行った。

 

インストール

ubuntu18.04LTSでテストした。

本体 Github

 

#bioconda (link)
conda create -n mobster-env -c bioconda -y mobster
conda activate mobster-env 

mobster -h

$ mobster -h

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

#MOBSTER                 #

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

Version: 0.2.4

Author: Djie Tjwan Thung

 

Predict non-reference Mobile Element Insertion (MEI) events using one properties file.

-properties [properties]

-in [input .bam file]. This value will override corresponding value in properties file. Multiple BAM files may be specified if seperated by a comma

-out [output prefix]. This value will override corresponding value in properties file.

-sn [sample name]. This value will override corresponding value in properties file. Multiple sample names may be specified if seperated by a comma

Default mapping tool: unspecified

0 [main] INFO Mobster  - Invalid arguments. Please try again.

 

 

実行方法

1サンプル。bamを指定する。

mobster \
-properties Mobster_latest.properties \
-in input.bam \
-sn test_sample \
-out mobster_test

 

複数サンプル

#family trio
mobster
\
-properties Mobster_latest.properties \
-in A1_child.bam,A1_father.bam,A1_mother.bam \
-sn A1_child,A1_father,A1_mother \
-out A1_trio_mobster

引用

Mobster: accurate detection of mobile element insertions in next generation sequencing data

Thung DT, de Ligt J, Vissers LE, Steehouwer M, Kroon M, de Vries P, Slagboom EP, Ye K, Veltman JA, Hehir-Kwa JY

Genome Biol. 2014;15(10):488.

 

SRAなどのシーケンシングデータを一括ダウンロードする grabseqs

2020 4/1  タイトル修正、誤字修正

2020 10/24 仮想環境を解くって導入するように修正

2021 5/23 conda => mambaに修正

 

 ハイスループットシーケンシングは、生物学的な疑問を解決するための強力な技術である。Grabseqsは、Sequence Read Archive(SRA)、Metagenomics Rapid Annotation through Subsystems Technology(MG-RAST)サーバー、iMicrobeを含む複数のリポジトリからデータとメタデータをダウンロードするための使いやすい単一のインターフェースを提供することで、一般に公開されているメタゲノムデータへのアクセスを効率化する。ユーザーは、1つの grabseqs コマンドで、任意のリポジトリから任意の数のサンプルやプロジェクトのデータやメタデータを標準化された形式でダウンロードすることができる。
 GrabseqsはPythonで実装され、MITライセンスの下でライセンスされたオープンソースのツールである。ソースコードhttps://github.com/louiejtaylor/grabseqsPython Package Index (PyPI)、Anaconda Cloud リポジトリから自由に入手できる。

 

インストール

ubuntu18.04LTSでpipを使ってテストした。

依存

  • Python 3 (external packages req'd: requests, requests-html, pandas, fake-useragent)

本体 Github

#grabseqsはpipにも対応(macosはpipのみ対応)
pip install grabseqs

#bioconda
mamba create -n grabseqs python=3.8 -y
conda activate grabseqs
mamba install grabseqs -c louiejtaylor -c bioconda -c conda-forge -y

#pigzやsra-toolsもないなら導入
mamba install pigz
mamba install -c bioconda -y sra-tools

> grabseqs --help

# grabseqs --help

usage: grabseqs [-h] [--version] {sra,imicrobe,mgrast} ...

 

Download metagenomic sequences from public datasets.

 

positional arguments:

  {sra,imicrobe,mgrast}

                        repositories available

    sra                 download from SRA

    imicrobe            download from iMicrobe

    mgrast              download from MG-RAST

 

optional arguments:

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

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

grabseqs sra -h

# grabseqs sra -h

usage: grabseqs sra [-h] [-m METADATA] [-o OUTDIR] [-r RETRIES] [-t THREADS]

                    [-f] [-l] [--parse_run_ids] [--use_fastq_dump]

                    [--custom_fqdump_args CUSTOM_FQD_ARGS] [--no_parsing]

                    id [id ...]

 

positional arguments:

  id                    One or more BioProject, ERR/SRR or ERP/SRP number(s)

 

optional arguments:

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

  -m METADATA           filename in which to save SRA metadata (.csv format,

                        relative to OUTDIR)

  -o OUTDIR             directory in which to save output. created if it

                        doesn't exist

  -r RETRIES            number of times to retry download

  -t THREADS            threads to use (for fasterq-dump/pigz)

  -f                    force re-download of files

  -l                    list (but do not download) samples to be grabbed

  --parse_run_ids       parse SRR/ERR identifers (do not pass straight to

                        fasterq-dump)

  --use_fastq_dump      use legacy fastq-dump instead of fasterq-dump (no

                        multithreaded downloading)

  --custom_fqdump_args CUSTOM_FQD_ARGS

                        'string' containing args to pass to fast(er)q-dump

  --no_parsing          Legacy option to not parse SRR IDs (now default)

grabseqs mgrast -h

# grabseqs mgrast -h

usage: grabseqs mgrast [-h] [-m METADATA] [-o OUTDIR] [-r RETRIES]

                       [-t THREADS] [-f] [-l]

                       rastid [rastid ...]

 

positional arguments:

  rastid       One or more MG-RAST project or sample identifiers

               (mgp####/mgm######)

 

optional arguments:

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

  -m METADATA  filename in which to save metadata (.csv format, relative to

               OUTDIR)

  -o OUTDIR    directory in which to save output. created if it doesn't exist

  -r RETRIES   number of times to retry download

  -t THREADS   threads to use (for pigz)

  -f           force re-download of files

  -l           list (but do not download) samples to be grabbed

grabseqs imicrobe -h

# grabseqs imicrobe -h

usage: grabseqs imicrobe [-h] [-m METADATA] [-o OUTDIR] [-r RETRIES]

                         [-t THREADS] [-f] [-l]

                         imicrobeid [imicrobeid ...]

 

positional arguments:

  imicrobeid   One or more iMicrobe project or sample identifiers (p##/s###)

 

optional arguments:

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

  -m METADATA  filename in which to save metadata (.csv format, relative to

               OUTDIR)

  -o OUTDIR    directory in which to save output. created if it doesn't exist

  -r RETRIES   number of times to retry download

  -t THREADS   threads to use (for pigz)

  -f           force re-download of files

  -l           list (but do not download) samples to be grabbed

 

 

実行方法

SRAの全サンプル(All runs)をダウンロードする。メタデータもダウンロードする。

grabseqs sra -t 8 -m metadata.csv -o outdir SRP#######
  • -o    directory in which to save output. created if it doesn't exist
  • -r     number of times to retry download
  • -t     threads to use (for fasterq-dump/pigz) 
  • -m   filename in which to save SRA metadata (.csv format, relative to OUTDIR)

 

複数指定にも対応している。さらにSRA/ERP ProjectとBioProjects、runsの混合も可能。

grabseqs sra -t 8 -m metadata.csv -o outdir -r 2 SRR######## ERP####### PRJNA######## ERR########

 

Run内の単一のシーケンシングデータだけダウンロードする。

grabseqs sra -t 8 -m metadata.csv -o outdir SRR6032562

#複数
grabseqs sra -t 8 -m metadata.csv -o outdir SRR6032562 SRR6032563 SRR6032564

 

ドライラン

grabseqs sra -l SRP########
  • -l    list (but do not download) samples to be grabbed 

>  grabseqs sra -l PRJNA590266

SRR10488335.fastq.gz

SRR10488336.fastq.gz

SRR10488337.fastq.gz

SRR10488338.fastq.gz

SRR10488339.fastq.gz

SRR10488340.fastq.gz

 

MG-RASTのプロジェクトをダウンロードする場合、サブコマンド"mgrast"を使う。

grabseqs mgrast -t 8 -m metadata.csv -o outdir -r 2 mgp##### mgm#######

 注;MG-RASTのプロジェクトの多くは、一般には公開されていません。もし、特定のアクセッション番号やプロジェクトに問題がある場合は、まずMG-RASTのウェブサイトに行き、手動でダウンロードできるかどうか確認してください。以前はダウンロード可能であった多くのサンプルが現在では入手不可能になっているため、まずこれを実行してください(FAQより)。

 

 

 

iMicrobeの例はGithubで確認して下さい。

引用

grabseqs: Simple downloading of reads and metadata from multiple next-generation sequencing data repositories
Louis J Taylor, Arwa Abbas, Frederic D Bushman
Bioinformatics, btaa167, Published: 10 March 2020

 

関連


 


 

 

 

 

 

pyGenomeTracks

 2020 8/4 論文引用追加

 

 出版可能な品質の複数のゲノムトラックを表示するプロットの生成は深刻な課題となる。望ましい正確な数値を作成するには、かなりの労力が必要である。これは通常、手作業で、またはベクターグラフィックソフトウェアを使用して行われる。
 pyGenomeTracks (PGT)は、複数のトラックを簡単に組み合わせることができるモジュール式のプロットツールである。再現性が高く、標準化された、高度にカスタマイズされた出版準備の整ったイメージの生成を可能にする。
 PGT は https://usegalaxy.eu のグラフィカルインターフェースとコマンドラインから利用できる。PGT は bioconda チャンネル経由で conda で提供され、pip で提供され、github: https://github.com/deeptools/pyGenomeTracks でオープンに開発されている。

 

pyGenomeTracksは、高度にカスタマイズ可能な高品質のゲノムブラウザトラックを作成することを目的としている。Hi-Cデータの他、以下のデータからplotを作成することが可能になっている。

  • bigwig
  • bed/gtf (many options)
  • bedgraph
  • epilogos
  • narrow peaks
  • links (represented as arcs)
  • Hi-C matrices

 

Documentation

example

Galaxy

 

インストール

依存

  • Python >=3.6
  • numpy >= 1.16
  • intervaltree >=2.1.0
  • pyBigWig >= 0.3.4
  • hicmatrix >= 0.14
  • pysam >= 0.8
  • matplotlib >= 3.1.1
  • gffutils >=0.9

本体 Github

https://github.com/deeptools/pyGenomeTracks

#bioconda (link)
conda create -n pygenometracks-env -y
conda activate pygenometracks-env
conda install -c bioconda -c conda-forge pygenometracks

#pip
pip install pyGenomeTracks

pygenometracks -h

$ pygenometracks -h

usage: pygenometracks --tracks tracks.ini --region chr1:1000000-4000000 -o image.png

 

Plots genomic tracks on specified region(s). Citation : Ramirez et al. High-

resolution TADs reveal DNA sequences underlying genome organization in flies.

Nature Communications (2018) doi:10.1038/s41467-017-02525-w

 

optional arguments:

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

  --tracks TRACKS       File containing the instructions to plot the tracks.

                        The tracks.ini file can be genarated using the

                        `make_tracks_file` program.

  --region REGION       Region to plot, the format is chr:start-end

  --BED BED             Instead of a region, a file containing the regions to

                        plot, in BED format, can be given. If this is the

                        case, multiple files will be created using a prefix

                        the value of --outFileName

  --width WIDTH         figure width in centimeters

  --height HEIGHT       Figure height in centimeters. If not given, the figure

                        height is computed based on the heights of the tracks.

                        If given, the track height are proportionally scaled

                        to match the desired figure height.

  --title TITLE, -t TITLE

                        Plot title

  --outFileName OUTFILENAME, -out OUTFILENAME

                        File name to save the image, file prefix in case

                        multiple images are stored

  --fontSize FONTSIZE   Font size for the labels of the plot

  --dpi DPI             Resolution for the image in case the ouput is a raster

                        graphics image (e.g png, jpg)

  --trackLabelFraction TRACKLABELFRACTION

                        By default the space dedicated to the track labels is

                        0.05 of the plot width. This fraction can be changed

                        with this parameter if needed.

  --trackLabelHAlign {left,right,center}

                        By default, the horizontal alignment of the track

                        labels is left. This alignemnt can be changed to right

                        or center.

  --version             show program's version number and exit

 

 

実行方法

1、trackファイルの作成

git clone https://github.com/deeptools/pyGenomeTracks.git
cd pyGenomeTracks/examples/
make_tracks_file --trackFiles input1.bed input2.bw -o tracks.ini

2、trackファイルと描画領域を指定してpyGenomeTracksを実行

pyGenomeTracks --tracks track.ini --region chr1:1000-20000 -o bigwig.png

 

 

テストラン

exampleランのコードをそのまま走らせてみる。

pyGenomeTracks --tracks bigwig_track.ini --region X:2,500,000-3,000,000 -o bigwig.png

f:id:kazumaxneo:20200330225924p:plain

 

pyGenomeTracks --tracks bigwig_with_genes.ini --region X:2,800,000-3,100,000 -o bigwig_with_genes.png 

f:id:kazumaxneo:20200330230018p:plain

 

pyGenomeTracks --tracks narrow_peak2.ini --region X:2760000-2802000 --trackLabelFraction 0.2 --dpi 130 -o master_narrowPeak23.png

f:id:kazumaxneo:20200330230223p:plain

 

pyGenomeTracks --tracks hic_track.ini -o hic_track.png --region chrX:2500000-3500000

f:id:kazumaxneo:20200330230355p:plain

 

引用

High-resolution TADs reveal DNA sequences underlying genome organization in flies

Fidel Ramírez, Vivek Bhardwaj, Laura Arrigoni, Kin Chung Lam, Björn A. Grüning, José Villaveces, Bianca Habermann, Asifa Akhtar & Thomas Manke
Nature Communications volume 9, Article number: 189 (2018)

 

2020 8/4

pyGenomeTracks: reproducible plots for multivariate genomic data sets
Lucille Lopez-Delisle, Leily Rabbani, Joachim Wolff, Vivek Bhardwaj, Rolf Backofen, Björn Grüning, Fidel Ramírez, Thomas Manke
Bioinformatics, Published: 03 August 2020

 

関連


 

 

パンゲノムにおいて有意な関連性の遺伝子を検出する Coinfinder

 

 原核生物および真核生物のパンゲノムのアクセサリー遺伝子は、遺伝子水平伝播、loss of gene、および選択の影響により蓄積する。 Coinfinderは、パンゲノム内の相同な遺伝子(遺伝子ファミリー)のセットが偶然に予想されるよりも頻繁に相互に関連または解離するかどうかを評価するソフトウェアプログラムである。 Coinfinderは、ユーザーが提供する系統樹を使用して、各アクセサリー遺伝子の系統依存性(系統発生分布)を評価する。これにより、Coinfinderは、偶然同じクレードに出現したために存在しない同時発生の遺伝子ペアに集中できる。むしろ、それらは系統発生全体で予想されるよりも頻繁に一緒に表示される傾向がある。 CoinfinderはC++、Python3およびRで実装されており、GNUライセンスの下でhttps://github.com/fwhelan/coinfinderから無料で入手できる。

 

インストール

依存

  • Cmake3.6 or greater 
  • Python3.6 or greater
  • Boost1.66 or greater 
  • OpenMP
  • gcc 5 or greater 
  • R libraries: caper, phytools, getopt, igraph, dplyr, cowplot, data.table, ggraph, flock, future
  • Bionconductor R library: ggtree 

本体 Github

#bioconda (link)
conda create -n coinfinder-env -y
conda activate coinfinder-env
conda install -c conda-forge -c bioconda -c defaults coinfinder -y

coinfinder -h

$ coinfinder -h

./confinder [OPTIONS]

File input- specify either: 

    -i or --input          The path to the gene_presence_absence.csv output from Roary

                           -or-

                           The path of the Alpha-to-Beta file with (alpha)(TAB)(beta)

    -I or --inputroary     Set if -i is in the gene_presence_absence.csv format from Roary

    -p or --phylogeny      Phylogeny of Betas in Newick format (required)

Max mode (mandatory for coincidence analysis):

    -a or --associate      Overlap; identify groups that tend to associate/co-occur (default).

    -d or --dissociate     Separation; identify groups that tend to dissociate/avoid.

Significance- specify: 

    -L or --level          Specify the significnace level cutoff (default: 0.05)

Significance correction- specify: 

    -m or --bonferroni     Bonferroni correction multiple correction (recommended & default)

    -n or --nocorrection   No correction, use value as-is

    -c or --fraction       (Connectivity analysis only) Use fraction rather than p-value

Alternative hypothesis- specify: 

    -g or --greater        Greater (recommended & default)

    -l or --less           Less

    -t or --twotailed      Two-tailed

Miscellaneous:

    -x or --num_cores      The number of cores to use (default: 2)

    -v or --verbose        Verbose output.

    -r or --filter         Permit filtering of saturated and low-abundance data.

    -U or --upfilthreshold Upper filter threshold for high-abundance data filtering (default: 1.0 i.e. any alpha in >=100/% of betas.

    -F or --filthreshold   Threshold for low-abundance data filtering (default: 0.05 i.e. any alpha in <=5% of betas.

    -q or --query          Query a specific gene.

    -T or --test           Runs the test cases and exits.

    -E or --all            Outputs all results, regardless of significance.

Output:

    -o or --output         The prefix of all output files (default: coincident).

 

 

If you use Coinfinder, please cite:

 

FJ Whelan, M Rusilowicz, & JO McInerney. "Coinfinder: Detecting Significant Associations and Dissociations in Pangenomes." doi: https://doi.org/10.1101/859371

 

 

テストラン

roaryで解析して得たgene_presence_absence.csvとnewick formatのツリーファイルを指定する。

git clone https://github.com/fwhelan/coinfinder-manuscript.git
cd coinfinder-manuscript/
coinfinder -i gene_presence_absence.csv -I -p core-gps_fasttree.newick -o out --associate
  • -i     The path to the gene_presence_absence.csv output from Roary
  • -I     Set if -i is in the gene_presence_absence.csv format from Roary
  • -x    The number of cores to use (default: 2)
  • -p     Phylogeny of Betas in Newick format (required)
  • --associate   Overlap; identify groups that tend to associate/co-occur (default).
  • --dissociate     Separation; identify groups that tend to dissociate/avoid.

出力

f:id:kazumaxneo:20200329170647p:plain

_pairs.tsvが有意な一致遺伝子ペアのタブ区切りリストになる。 

 

ネットワークファイル.GEXF v1.2(Graph Exchange XML Format) では、各遺伝子(ノード)が統計的に互いに共起している場合は、各遺伝子がエッジで別の遺伝子に接続される。ノードは系統における系統非依存性によって重み付けされる(すなわち、ノードが大きいほど系統的に独立した遺伝子である)。ノードは、接続されたコンポーネント、または互いに関連性を持つ遺伝子のセットによって色分けされる。GEXFはGEXFに対応したネットワークビューアで可視化できる。

このデータは、系統に関連した有無のヒートマップとして表示することもできる 。

out_heatmap.pdf

f:id:kazumaxneo:20200329170953p:plain

拡大

f:id:kazumaxneo:20200329170935p:plain

ヒートマップ内の遺伝子は、D値の高い順(最も系統に依存しないものから低いものまで)に並べられ、一致パターンに応じて色分けされる。ヒートマップは見やすくするために必要に応じて複数のファイルに分割される。

 

すべてのユニークな一致遺伝子とそのD値のリストはout_nodes.tsvにまとめられる。

f:id:kazumaxneo:20200329171516p:plain

 

引用

Coinfinder: detecting significant associations and dissociations in pangenomes Open Access
Fiona Jane Whelan, Martin Rusilowicz, James Oscar McInerney

Microbial Genomics, Published: 25 February 2020

 

関連

データをドラッグアンドドロップして素早く作図とデータ分析ができるGUIツール BioVinci

3/330 誤字修正、図を追記

2020 3/31 追記

2020 5/4 v2.0リリース追記

2020 5/12 統合TVのリンク追加

 

 

BioVinciのabout usより

BioTuringは、若くてモチベーションが高いバイオインフォマティシャンで構成されたチームで、次世代のバイオインフォマティクスツールを使ってライフサイエンスを育成したいという情熱を全員が共有している。2016年に設立されたこのチームは急速に動き出し、現在では製薬会社や学術機関にアナリティクスソリューションを提供している。

 人の病気を治す方法を見つけたい、自然を解読したいという私たちの究極の願いが、困難な試練を乗り越えるために私たちを支えている。チームとして、私たちは失敗からも成果からも学び、重要なバイオインフォマティクスの問題を解決するために最先端のアルゴリズムとソフトウェアを提供したいという熱い情熱を持っている。

 

 

 

2020 5/4

 

 

 User guide

https://vinci.bioturing.com/user-guide

Q&A

https://vinci.bioturing.com/FAQ

 

2020 5/12 

統合TVから解説した動画が公開されています。

 

利用方法

動画がたくさんあるので、ここでは簡単に流れだけ紹介します。macos10.14でテストしました。

 

注意;ライセンスについて色々調べたのですが理解しきれず、自分が試した時はversion 1.1.5(もしかすると本来は有料版?)がダウンロードできたのですが、現在はVersion: 0.1.1 betaのみダウンロードできるようです。ということで、ここではversion 1.1.5を紹介しますが、0.1.1 betaとはインターフェイスが異なっているので注意して下さい。 

2020 5/4

v2.0rリリースのお知らせメールが来ました。ダウンロードできます。試してみて下さい。

 

HP

https://vinci.bioturing.com にアクセスする。

f:id:kazumaxneo:20200328083251p:plain

ログインしてブラウザ上で使うか、ローカルアプリをダウンロードする。ここではローカル版を使う。DOWNLOADボタンをクリック。

 

メールアドレス等を記載してDOWNLOADボタンをクリック。

f:id:kazumaxneo:20200328083834p:plain

macosとwindows10が選べる。macosを選んだ。

 

ダウンロードしたら自己解凍形式の.dmgを実行する。 

f:id:kazumaxneo:20200328084326p:plain

アプリをapplicationsにコピーする。必要に応じてDockにエイリアスも追加しておく。

起動したところ(version 1.1.5 )

f:id:kazumaxneo:20200328084428p:plain

 

初期はExampleのタブになっている。

f:id:kazumaxneo:20200328085450p:plain

何か表示してみる。Box plot with annotationsを選択。

f:id:kazumaxneo:20200328123434p:plain

 

表示された。

f:id:kazumaxneo:20200328123528p:plain

 

拡大するには右下のスケールボタンを押す。

f:id:kazumaxneo:20200328123612p:plain

 

左上のBackから初期画面に戻った。自分のデータを使って作図するにはWorkSetsタブに移動する。 一度exampleのデータを表示しておくと、Worksetsタブでそのグラフがテンプレート表示されるようになる。

f:id:kazumaxneo:20200328100445p:plain

ここではBar plotを選択した。グラフの種類を決めず、新しいプリセットを作ることもできる(右上のプラスボタンから)。

 

白紙の画面に切り替わったら、左上のメニューにあるAdd dataボタンをクリックしてデータを追加する。

f:id:kazumaxneo:20200328101638p:plain

 

 

その場でデータを打ち込む。または、

f:id:kazumaxneo:20200328101839p:plain

 

Uploadに切り替えてデータファイルをウィンドウ内にドラッグアンドドロップする。認識可能なファイルフォーマットは、タブ区切りのTSVファイルかコンマ区切りのCSVファイル、そしてExcel形式(.xlsx)になる。

f:id:kazumaxneo:20200328101614p:plain

準備したExcelファイルをウィンドウ内にドラッグアンドドロップした。

 

=====================================================

参考

今回読み込んだデータ(RNA seqのリードカウントファイルの一部)

f:id:kazumaxneo:20200328102829p:plain

=====================================================

 

白紙の画面に戻る。

f:id:kazumaxneo:20200328102305p:plain

 

左端にシートとデータが表示される。data1をvalue(s)の項目にドラッグアンドドロップした。

f:id:kazumaxneo:20200328102954p:plain

 

valueがdata1の棒グラフとして視覚化された。

f:id:kazumaxneo:20200328103003p:plain

 

続いてdata2を同じvalue (s)にプロットする。

f:id:kazumaxneo:20200328103314p:plain

隣に並んで視覚化された。

f:id:kazumaxneo:20200328103316p:plain

 

data3-5のみ視覚化した。

f:id:kazumaxneo:20200328103528p:plain

 

図はインタラクティブに操作可能。データと重なってしまっている凡例をドラッグして右端に移動した。

f:id:kazumaxneo:20200328103717p:plain

 

凡例のdata3をクリックして一時的に表示をオフにした。

f:id:kazumaxneo:20200328103814p:plain

 

グラフの上で右クリックすることで表示形式を変更できる。

f:id:kazumaxneo:20200328104310p:plain

 

Position => Fill

f:id:kazumaxneo:20200328104329p:plain

 

Position => Stack

f:id:kazumaxneo:20200328104516p:plain

 

Position => Overlap

f:id:kazumaxneo:20200328104536p:plain

 

元に戻す。

Position => Side by side

f:id:kazumaxneo:20200328104623p:plain

 

 Components => Densitiy plot

f:id:kazumaxneo:20200328104715p:plain

 

Components => Both

f:id:kazumaxneo:20200328104752p:plain

 

右クリック => Scale から対数グラフに変更できる。

f:id:kazumaxneo:20200328132814p:plain

 

 

 

 

目盛りやフォントなどを修正していく。グラフの左上の赤ボタンを押し、編集モードにする。

f:id:kazumaxneo:20200328131425p:plain

 

目盛のフォントの種類、サイズ、文字の角度、ボールド・斜体など変更できる。

f:id:kazumaxneo:20200328131451p:plain

文字をドラッグして移動することもできる。 

 

編集モードでは、背景色や、
f:id:kazumaxneo:20200328131908p:plain

 

グリッドの線の有無も変更可能。

f:id:kazumaxneo:20200328131912p:plain

 

バーの背景やボーダーの色も変えられる。

f:id:kazumaxneo:20200328133126p:plain

 

新しくタイトルテキストなどを加えたい場合、一番上のAdd textをクリックする。

f:id:kazumaxneo:20200328133440p:plain

 

一番上のDoneを押すと編集モードを抜ける。

f:id:kazumaxneo:20200328131730p:plain

 plotサイズは変更できないようなので、右下からサイズを小さくして相対的にプロットを大きくする。

 

上のexportボタンから様々な形式で出力できる。

f:id:kazumaxneo:20200328130655p:plain

 

TOPに戻ろうとするとセーブするか聞いてくる。セーブを選択して戻り、Worksets に行くと、セーブしたグラフが表示されるようになる。

f:id:kazumaxneo:20200328131039p:plain

 

開く。初期は白紙だが、右上にスナップショットが表示されているので、

f:id:kazumaxneo:20200328131158p:plain

 

クリックすると最新の状態ですぐ再開できる。

f:id:kazumaxneo:20200328131247p:plain

このように、複数データある場合は右上のスナップショットを選んで再開する。

 

 

さらに機能を説明していく。グラフの種類は下のボタンから瞬時に変更できる。

f:id:kazumaxneo:20200328133827p:plain

 

グラフを変更する時、セーブするか聞いてくるので必ずセーブを選択する(元のグラフに戻してもカスタム設定がリセットされてしまうため)。

 

Pie chart

f:id:kazumaxneo:20200328134005p:plain

 

Violin plot

f:id:kazumaxneo:20200328134011p:plain

 

右下のテーマボタンからは、配色が統一されたテーマに変えられる。

f:id:kazumaxneo:20200328134350p:plain

 

テーマClsssic

f:id:kazumaxneo:20200328134446p:plain

 

テーマ;Journal(Color) (別のデータです)

f:id:kazumaxneo:20200328134606p:plain

 

テーマ;Journal(Black and white)

f:id:kazumaxneo:20200328134657p:plain

 

example

Violin plot 

f:id:kazumaxneo:20200330095834p:plain


Heatmap
f:id:kazumaxneo:20200330094554p:plain

 

Scatter plot

f:id:kazumaxneo:20200330094621p:plain

 

Welch's t test

f:id:kazumaxneo:20200330095406p:plain

(Student's t-testもあり)

 

One-way ANOVA

f:id:kazumaxneo:20200330094819p:plain

(Two-way ANOVAもあり)

 

そのほかの機能は動画を見て確認して下さい。統計関係の解説は 上にリンクを張ったUser guideに詳しく載っています。t検定のほか、回帰分析、各種クラスタリング、分散分析(ANOVA)などにも対応しています。

 

追記

作図や統計に関してはこの辺の記事も参考になるんじゃないかと思います。

引用

https://vinci.bioturing.com

 

 

 

 

 

 

単一のメタゲノムアセンブリゲノム(MAGs)とシーケンシングデータからバクテリアの増殖率を推定する iRep

 

 培養に依存しない微生物群集の研究により、微生物群集の複雑さと代謝の可能性に対する理解が深まった。ただし、コミュニティへの個々のマイクロバイオームメンバーの貢献を理解するには、どの細菌が活発に複製しているかを判断することが重要になる。ドラフト品質のゲノムシーケンスと単一のタイムポイントでのメタゲノムシーケンスを使用して微生物集団の複製率を推定するアルゴリズム、iRepを開発した。アルゴリズムは、単一の複製元からの双方向ゲノム複製から生じるシーケンスカバレッジトレンドに基づいて、複製のインデックス(iRep)を計算する。微生物の複製率がヒト幼児の抗生物質投与後に増加することを示すため、この方法を適用した。また、未培養の地下水に関連する候補フィラ放射線バクテリアが、地球化学の実質的な変化を経験している地下コミュニティでまれにしか急速に複製しないことも示す。この方法は、さまざまな条件に対する生物の反応を追跡し、活発に成長している個体群を特定し、研究のモデリングに使用する複製率を測定するために、すべてのゲノム関係マイクロバイオーム研究に適用できる。
 自然集団の分裂細胞は、平均して複数コピーのゲノムを含んでいる(論文図1)。成長している非同期の細菌集団では、細胞にはさまざまな程度に複製されるゲノムが含まれており、複製の起点から終点までの平均ゲノムコピー数が徐々に減少する(ref.1)。この減少は、完全なゲノム全体のDNAシーケンスカバレッジの変化を測定することで検出できる。バクテリアのゲノム複製は単一の複製起点から双方向に進行するため、複製の起点と終点はこのカバレッジパターンに基づいて推定できる(ref.2)。多種多様な細菌のGC skewおよびゲノムカバレッジの分析により、この複製メカニズムが広く適用可能であることが示されている。さらに、細菌培養のアーリィステージの研究では、細胞がゲノム複製を複数ラウンド同時に開始することにより、より速い分裂を達成できることが明らかになった(ref.9)。
 Koremらは細菌の複製率を測定するために、複製の終点と比較した起点でのシーケンスカバレッジの比率を使用した(ref.8)。複製起点と終点がそれぞれカバレッジピークとトラフに対応するため、著者はメソッドにPTR(peak-to-trough ratio)という名前を付けた。彼らはPTRを適用して、ヒトマイクロバイオーム内の特定の細菌の複製率を計算したが、シーケンシングリードを対象の細菌の完全な閉じた環状リファレンスゲノムにマッピングするための要件は、大きな制限である。細菌の大部分は未培養のままで、リファレンスゲノムが不足している。
 メタゲノミクスでは、リファレンスゲノムを欠く細菌や古細菌のドラフトゲノムを定期的に生成する(論文図1および補足図1)。多くの場合、これらの生物はほとんど知られていない微生物門からのものであり、データベースに完全なゲノムがある生物とは大きく異なる。単一のエコシステムから数百または数千のドラフトまたはほぼ完全なゲノムを回復できる場合がある。カバレッジベースのレプリケーションレート分析を拡張して、これらのドラフトゲノムのシーケンスカバレッジトレンドに基づく測定を可能にする方法を紹介する。フラグメントの順序が不明であるという事実にもかかわらず、この方法は機能する。 PTRとは異なり、このアプローチは、、大多数の細菌の完全なゲノムが利用できない土壌などの複雑なシステムを含む、事実上すべての自然または人工のエコシステムに適用できる。

 

https://twitter.com/search?q=iRep%20genome&src=typed_query

 

テストラン

macos101.4でテストした()。

本体 Github

https://github.com/christophertbrown/iRep

#bioconda (link)
conda create -n iRep -y
conda activate iRep
conda install -c bioconda -y iRep

#pip
pip install iRep

> iRep -h

$ iRep -h

usage: iRep [-h] -f [F [F ...]] -s [S [S ...]] -o O [--pickle] [-mm MM] [--sort] [-M M] [--no-plot] [--no-gc-correction] [-ff] [-t T]

 

# calculate the Index of Replication (iRep)

 

optional arguments:

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

  -f [F [F ...]]      fasta(s)

  -s [S [S ...]]      sorted sam file(s) for each sample (e.g.: bowtie2 --reorder)

  -o O                prefix for output files (table and plots)

  --pickle            save pickle file (optional)

  -mm MM              max. # of read mismatches allowed (default: 1)

  --sort              optional - sort the sam file

  -M M                max. memory (GB) for sorting sam (default: 100)

  --no-plot           do not plot output

  --no-gc-correction  do not correct coverage for GC bias before calculating iRep

  -ff                 overwrite files

  -t T                threads (default: 6)

bPTR -h

$ bPTR -h

usage: bPTR [-h] [-f [F [F ...]]] [-s [S [S ...]]] -m M [-c C] -o O [-pickle PICKLE] -plot PLOT [-mm MM] [-p P] [--sort] [-b B] [-ff] [-t T]

 

# est. growth rate from peak-to-trough coverage ratio

 

optional arguments:

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

  -f [F [F ...]]  fasta(s)

  -s [S [S ...]]  sorted sam file(s) for each sample (e.g.: bowtie2 --reorder)

  -m M            method for detecting Ori/Ter of replication: gc_skew or coverage

  -c C            pre-computed data from growth_ptr.py (optional: pickle file)

  -o O            filename for output table

  -pickle PICKLE  filename for output pickle file (optional)

  -plot PLOT      filename for coverage profile plots (default: no plots)

  -mm MM          maximum number of mapping mismatches allowed (default: no limit)

  -p P            number of permutations to perform (default: None)

  --sort          sort the sam file

  -b B            max memory (GB) for sorting sam (default: 100)

  -ff             overwrite files

  -t T            threads (default: 6)

 

 

 

テストラン

bPTRは完全長のゲノムに、iRepは高品質のドラフトゲノム(>=75% complete, <=175 fragments/Mbp sequence, and <=2% contamination)に使う。

 

iRepのラン。ドラフトゲノムのFASTAファイルとbowtie2でソート(--reorderをつける)してアラインメントしたSAMファイルを指定する。

git clone https://github.com/christophertbrown/iRep.git
cd iRep/
iRep -f sample_data/l_gasseri.fna -s sample_data/l_gasseri*sam -o test.iRep
  • -f     fasta(s)
  • -s    sorted sam file(s) for each sample (e.g.: bowtie2 --reorder)
  • -o    prefix for output files (table and plots)

出力

結果のtsvファイルと複製率計算に使用されたゲノムカバレッジプロットを示すPDFが出力される。

> cat test.iRep.tsv

$ cat test.iRep.tsv 

## index of replication (iRep) - thresholds: min cov. = 5, min wins. = 0.98, min r^2 = 0.9, max fragments/Mbp = 175, GC correction min r^2 = 0.0

# genome sample_data/l_gasseri.fna-vs-l_gasseri_sample1-shrunk.sam sample_data/l_gasseri.fna-vs-l_gasseri_sample2-shrunk.sam

sample_data/l_gasseri.fna 1.906559215638506 2.361137517138996

#

## un-filtered index of replication (iRep)

# genome sample_data/l_gasseri.fna-vs-l_gasseri_sample1-shrunk.sam sample_data/l_gasseri.fna-vs-l_gasseri_sample2-shrunk.sam

sample_data/l_gasseri.fna 1.906559215638506 2.361137517138996

#

## raw index of replication (no GC bias correction)

# genome sample_data/l_gasseri.fna-vs-l_gasseri_sample1-shrunk.sam sample_data/l_gasseri.fna-vs-l_gasseri_sample2-shrunk.sam

sample_data/l_gasseri.fna 1.9181206042659673 2.380084859373945

#

## r^2

# genome sample_data/l_gasseri.fna-vs-l_gasseri_sample1-shrunk.sam sample_data/l_gasseri.fna-vs-l_gasseri_sample2-shrunk.sam

sample_data/l_gasseri.fna 0.9988996508720132 0.9942789514575139

#

## coverage

# genome sample_data/l_gasseri.fna-vs-l_gasseri_sample1-shrunk.sam sample_data/l_gasseri.fna-vs-l_gasseri_sample2-shrunk.sam

sample_data/l_gasseri.fna 9.352049457279216 8.93781095577987

#

## % windows passing filter

# genome sample_data/l_gasseri.fna-vs-l_gasseri_sample1-shrunk.sam sample_data/l_gasseri.fna-vs-l_gasseri_sample2-shrunk.sam

sample_data/l_gasseri.fna 100.00 100.00

#

## fragments/Mbp

# genome sample_data/l_gasseri.fna-vs-l_gasseri_sample1-shrunk.sam sample_data/l_gasseri.fna-vs-l_gasseri_sample2-shrunk.sam

sample_data/l_gasseri.fna 1 1

#

## GC bias

# genome sample_data/l_gasseri.fna-vs-l_gasseri_sample1-shrunk.sam sample_data/l_gasseri.fna-vs-l_gasseri_sample2-shrunk.sam

sample_data/l_gasseri.fna 0.062006224268363996 0.000816622280244355

#

## GC r^2

# genome sample_data/l_gasseri.fna-vs-l_gasseri_sample1-shrunk.sam sample_data/l_gasseri.fna-vs-l_gasseri_sample2-shrunk.sam

sample_data/l_gasseri.fna 0.009408843274407497 0.0004338746957109896

#

test.iRep.pdf

f:id:kazumaxneo:20200327192837p:plain

f:id:kazumaxneo:20200327192841p:plain

f:id:kazumaxneo:20200327192846p:plain


closely relatedなゲノムが存在するメタゲノムのデータからアセンブルしていることを想定している。bPTRもiRepも単離した菌のシーケンシングデータ向けではないので注意してください。
引用

Measurement of bacterial replication rates in microbial communities

Brown CT, Olm MR, Thomas BC, Banfield JF

Nat Biotechnol. 2016 Dec;34(12):1256-1263