macでインフォマティクス

macでインフォマティクス

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

リボソームプロファイリングのクオリティメトリクスを提供する MappingQC

 

MappingQCは、リボソームプロファイリングデータのマッピングの品質の概要を示すいくつかの図を簡単に生成するツールである。 より具体的には、 P site offsetの計算、遺伝子分布、およびメタジェニック分類の概要を示す。 さらに、MappingQCは、データの標準的なトランスクリプトのトリプレット周期性とリンクされたトリプレットフェーズ(リボソームプロファイリングに典型的)の徹底的な分析を行う。 特に、phase distributionとRPFの長さ、リンクの相対的な位置、およびトリプレット同一性の間のリンクが考慮される。

 

Galaxy版(link)とローカル版がある。ここではローカル版使用の流れを簡単にまとめます。


インストール

依存(Githubより)

MappingQC relies on following Perl modules which have to be installed on your system:

  • DBI
  • Getopt::Long
  • Parallel::ForkManager
  • CWD
  • Data::Dumper (for debugging purposes)

Furthermore, mappingQC relies on following Python2 modules which have to be installed on your system:

  • getopt
  • defaultdict (collections)
  • sqlite3
  • pandas
  • numpy
  • matplotlib (including pyplot, colors, cm, gridspec, ticker and mplot3d)
  • seaborn

Github

#bioconda (link) ここでは仮想環境mappingqcに入れる
conda create -n mappingqc -c bioconda -y mqc python=2.7
conda activate mappingqc

map2bed.pl

$ map2bed.pl 

map2bed converts ART's map files to a BED file

 

USAGE: /usr/local/bin/map2bed.pl out_bed_file.bed in_map_file_1 [ in_map_file_2 ...]

 

(mqc) kazuma@kamisakumanoMBP:~/Downloads$ mQC.pl 

Working directory                                        : /Users/kazuma/Downloads

The following tmpfolder is used                          : /Users/kazuma/Downloads/tmp

 

 

MappingQC (Stand-alone version)

 

    MappingQC is a tool to easily generate some figures which give a nice overview of the quality of the mapping of ribosome profiling data. More specific, it gives an overview of the P site offset calculation, the gene distribution and the metagenic classification. Furthermore, MappingQC does a thorough analysis of the triplet periodicity and the linked triplet phase (typical for ribosome profiling) in the canonical transcript of your data. Especially, the link between the phase distribution and the RPF length, the relative sequence position and the triplet identity are taken into account.

        

    Input parameters:

    --help                  this helpful screen

    --work_dir              working directory to run the scripts in (default: current working directory)

    --experiment_name       customly chosen experiment name for the mappingQC run (mandatory)

    --samfile               path to the SAM/BAM file that comes out of the mapping script of PROTEOFORMER (mandatory)

    --cores                 the amount of cores to run the script on (integer, default: 5)

    --species               the studied species (mandatory)

    --ens_v                 the version of the Ensembl database you want to use

    --tmp                   temporary folder for storing temporary files of mappingQC (default: work_dir/tmp)

    --unique                whether to use only the unique alignments.

    Possible options: Y, N (default Y)

    --mapper                the mapper you used to generate the SAM file (STAR, TopHat2, HiSat2) (default: STAR)

    --maxmultimap           the maximum amount of multimapped positions used for filtering the reads (default: 16)

    --ens_db                path to the Ensembl SQLite database with annotation info. If you want mappingQC to download the right Ensembl database automatically for you, put in 'get' for this parameter (mandatory)

    --offset                the offset determination method.

                                Possible options:

                                - plastid: calculate the offsets with Plastid (Dunn et al. 2016)

                                - standard: use the standard offsets from the paper of Ingolia et al. (2012) (default option)

                                - from_file: use offsets from an input file

    --plastid_bam           the mapping bam file for Plastid offset generation (default: convert)

    --min_length_plastid    the minimum RPF length for Plastid offset generation (default 22)

    --max_length_plastid    the maximum RPF length for Plastid offset generation (default 34)

    --offset_file           the offsets input file

    --min_length_gd         minimum RPF length used for gene distributions and metagenic classification (default: 26).

    --max_length_gd         maximum RPF length used for gene distributions and metagenic classification (default: 34).

    --outfolder             the folder to store the output files (default: work_dir/mQC_output)

    --tool_dir              folder with necessary additional mappingQC tools. More information below in the dependencies section. (default: search for the default tool directory location in the active conda environment)

    --plotrpftool           the module that will be used for plotting the RPF-phase figure

                                Possible options:

                                - grouped2D: use Seaborn to plot a grouped 2D bar chart (default)

                                - pyplot3D: use mplot3d to plot a 3D bar chart. This tool can suffer sometimes from Escher effects, as it tries to plot a 3D plot with the 2D software of pyplot and matplotlib.

                                - mayavi: use the mayavi package to plot a 3D bar chart. This tool only works on local systems with graphical cards.

    --outhtml               custom name for the output HTML file (default: work_dir/mQC_experiment_name.html)

    --outzip                custom name for output ZIP file (default: work_dir/mQC_experiment_name.zip)

    

 

 

 

ERROR: do not forget the experiment name!

 

 

実行方法

samとデータベースを指定する。

mQC.pl --experiment_name yourexperimentname --samfile yoursamfile.sam --cores 20 --species human --ens_v 86 --ens_db ENS_hsa_86.db --unique N --offset plastid --plastid_bam yourbamfile.bam --tool_dir mqc_tools

 

引用

https://github.com/Biobix/mQC

 

RNA seqのクオリティメトリクスを提供する RNA-SeQC

2022/01/23 インストール追記

 

 RNA-seqは、細胞のトランスクリプトーム全体の特性評価を提供する。シーケンスのパフォーマンスとライブラリの品質の評価は、RNA-seqデータの解釈に不可欠だが、この問題に対処するツールはほとんない(論文執筆時点)。ここではデータ品質の重要な指標を提供するプログラムであるRNA-SeQCを紹介する。これらのメトリックには、歩留まり、アラインメント、duplicationの割合、GCバイアス、rRNAコンテンツ、アラインメント領域(エクソンイントロン、遺伝子内)、カバレッジの連続性、3 '/ 5'バイアス、検出可能な転写産物の数などが含まれる。このソフトウェアは、ライブラリ構築プロトコル、入力マテリアル、その他の実験パラメーターのマルチサンプル評価を提供する。ソフトウェアのモジュール性により、パイプラインの統合と、アラインできるリード数、duplication率、rRNA汚染などのデータ品質の主要な測定値の定期的な監視が可能になる。 RNA-SeQCにより、研究者は下流の分析にサンプルを含めるかどうかについて十分な情報に基づいた決定を下すことができる。要約すると、RNA-SeQCは、実験の設計、プロセスの最適化、および下流のコンピューター解析に不可欠な品質管理手段を提供する。

 

HP

https://software.broadinstitute.org/cancer/cga/rna-seqc

 

インストール

The latest stable build of RNA-SeQC is available on the GitHub Releases page, and contains static binaries for Linux and OSX.

依存

  • openjdk 7

本体 Github

#bioconda (link)
mamba create -n RNA-SeQC -y
conda activate RNA-SeQC
mamba install -c bioconda -y rna-seqc

rnaseqc

$ rnaseqc 

  rnaseqc [gtf] [bam] [output] {OPTIONS}

 

    RNASeQC 2.3.4

 

  OPTIONS:

 

      -h, --help                        Display this message and quit

      --version                         Display the version and quit

      gtf                               The input GTF file containing features

                                        to check the bam against

      bam                               The input SAM/BAM file containing reads

                                        to process

      output                            Output directory

      -s[sample], --sample=[sample]     The name of the current sample. Default:

                                        The bam's filename

      --bed=[BEDFILE]                   Optional input BED file containing

                                        non-overlapping exons used for fragment

                                        size calculations

      --fasta=[fasta]                   Optional input FASTA/FASTQ file

                                        containing the reference sequence used

                                        for parsing CRAM files

      --chimeric-distance=[DISTANCE]    Set the maximum accepted distance

                                        between read mates. Mates beyond this

                                        distance will be counted as chimeric

                                        pairs. Default: 2000000 [bp]

      --fragment-samples=[SAMPLES]      Set the number of samples to take when

                                        computing fragment sizes. Requires the

                                        --bed argument. Default: 1000000

      -q[QUALITY],

      --mapping-quality=[QUALITY]       Set the lower bound on read quality for

                                        exon coverage counting. Reads below this

                                        number are excluded from coverage

                                        metrics. Default: 255

      --base-mismatch=[MISMATCHES]      Set the maximum number of allowed

                                        mismatches between a read and the

                                        reference sequence. Reads with more than

                                        this number of mismatches are excluded

                                        from coverage metrics. Default: 6

      --offset=[OFFSET]                 Set the offset into the gene for the 3'

                                        and 5' windows in bias calculation. A

                                        positive value shifts the 3' and 5'

                                        windows towards eachother, while a

                                        negative value shifts them apart.

                                        Default: 150 [bp]

      --window-size=[SIZE]              Set the size of the 3' and 5' windows in

                                        bias calculation. Default: 100 [bp]

      --gene-length=[LENGTH]            Set the minimum size of a gene for bias

                                        calculation. Genes below this size are

                                        ignored in the calculation. Default: 600

                                        [bp]

      --legacy                          Use legacy counting rules. Gene and exon

                                        counts match output of RNA-SeQC 1.1.9

      --stranded=[stranded]             Use strand-specific metrics. Only

                                        features on the same strand of a read

                                        will be considered. Allowed values are

                                        'RF', 'rf', 'FR', and 'fr'

      -v, --verbose                     Give some feedback about what's going

                                        on. Supply this argument twice for

                                        progress updates while parsing the bam

      -t[TAG...], --tag=[TAG...]        Filter out reads with the specified tag.

      --chimeric-tag=[TAG]              Reads maked with the specified tag will

                                        be labeled as Chimeric. Defaults to 'mC'

                                        for STAR

      --exclude-chimeric                Exclude chimeric reads from the read

                                        counts

      -u, --unpaired                    Treat all reads as unpaired, ignoring

                                        filters which require properly paired

                                        reads

      --rpkm                            Output gene RPKM values instead of TPMs

      --coverage                        If this flag is provided, coverage

                                        statistics for each transcript will be

                                        written to a table. Otherwise, only

                                        summary coverage statistics are

                                        generated and added to the metrics table

      --coverage-mask=[SIZE]            Sets how many bases at both ends of a

                                        transcript are masked out when computing

                                        per-base exon coverage. Default: 500bp

      -d[threshold],

      --detection-threshold=[threshold] Number of counts on a gene to consider

                                        the gene 'detected'. Additionally, genes

                                        below this limit are excluded from 3'

                                        bias computation. Default: 5 reads

      "--" can be used to terminate flag options and force all following

      arguments to be treated as positional options

 

 

Argument validation error: No GTF file provided

 

テストラン

rnaseqc [gtf] [bam] [output] {OPTIONS}の順で指定する。

git clone --recursive https://github.com/broadinstitute/rnaseqc.git
cd rnaseqc/
rnaseqc test_data/downsampled.gtf test_data/downsampled.bam --bed test_data/downsampled.bed --coverage .
  • <gtf>          The input GTF file containing features.
  • <bam      The input SAM/BAM file containing reads to process.
  • <output   Output directory.

大きめのテストファイルはそのままではgit cloneされません。直接githubからダウンロードするかGit LFS をインストールしてから実行してください。 

出力

f:id:kazumaxneo:20191204010302p:plain

引用

RNA-SeQC: RNA-seq metrics for quality control and process optimization
David S. DeLuca,* Joshua Z. Levin, Andrey Sivachenko, Timothy Fennell, Marc-Danie Nazaire, Chris Williams, Michael Reich, Wendy Winckler, Gad Getz

Bioinformatics. 2012 Jun 1; 28(11): 1530–1532

 

関連


 

 

IRIS-EDA(立ち上げだけ紹介)

 

 生物学的システムの統合モデルを構築し、疾病を予防または治療するための実現可能な戦略を考案するには、適切な実験計画法とインタラクティブなインタフェースを備えた高度な計算ツールが必要である。 RNA-Seqは膨大な量の遺伝子発現データを作成しており、データ解析と解釈の需要は非常に大きい[ref.4]。遺伝子発現データの分析は、方法および実験を適切に設計し、そして多くの計算言語のうちの1つを使用して分析プロセスを実施する際のコンピューター経験によって促進される。しかしこれは、RNA-Seq研究を分析したいと思うユーザが限られた計算経験しか持っていないと障害となる。したがって、使いやすくインタラクティブな表現を持った分析と結果の視覚化方法の必要性が高まっている[ref.5]。

 サンプルレベルまたは細胞レベルでデータの特定の品質を決定するために、多種多様な計算方法を発現データに適用することができるが[ref.6-13]、ディファレンシャル遺伝子発現(DGE)分析が最も一般的に使用される。それは研究者が2つ以上の条件にわたって差次的に発現される遺伝子(DEG)を同定することを可能にし、そして遺伝子発現レベルの差を表現型の変動と相関させる意味のある方法を提供することができる。 DESeq [ref.14]、DESeq2 [ref.15]、edgeR [ref.16]、limma [ref.17]、Cuffdiff [ref.18]、Cuffdiff2 [ref.19]、sleuth [ref.20]など、多くのツールが開発および最適化されている。 DGE分析とDGE結果の視覚化にはかなりの努力が払われてきたがref.[21-28]、実験計画の課題、包括的な統合発見主導分析およびDGEツールの必要性、分析結果の視覚化に関連するインタラクティブな機能、機能性の欠如など、多くの落とし穴とボトルネックが残っている。

 これらのボトルネックに対処するために、本著者らはIRIS-EDAを作成した。これは発現データ解析のためのインタラクティブRNA-Seq解釈システムである。IRIS-EDAは遺伝子発現データを包括的に分析し、そしてインタラクティブな要約視覚化を容易に生成するための、ユーザーフレンドリーなインタラクティブプラットフォームを提供する。他の分析プラットフォームとは対照的に、IRIS-EDAはより包括的でマルチレベルの分析プラットフォームをユーザーに提供する。 IRIS-EDAは、1)シングルセルおよびバルクRNA-Seq分析機能、2)GEO submitの互換性、3)有用な発見主導およびDGE分析、7)4)7 DGE解析のための3つの統合されたツールと5)7つのインタラクティブな視覚化による実験計画法アプローチ(論文図1)。

 具体的には、IRIS-EDAは包括的なRNA-Seqデータ処理と解析をシームレスなワークフローで提供する。この調査アプローチでは、発現品質管理と、広く使用されている3つのRパッケージ、DESeq2、edgeR、およびlimmaの1つを通じてDGE分析と統合された発見主導型分析を使用する。それは、直感的な実験的設計オプション(例えば、対比較および要因比較、主効果およびグループ化主効果など)の選択、ならびにDGE分析においてカスタム設計マトリックスをアップロードするためのオプションをユーザーに提供する。 IRIS-EDAには、分析タイプごとにインタラクティブな視覚化機能が多数含まれているため、ユーザーはデータと結果を即座にグローバルに表示したり、出版物用の高解像度静止画像としてダウンロードすることができる。このツールは初めてFAIR Data Principles [ref.30]に基づいたフレームワークを実装し、ユーザーが自分のデータと結果をNCBIのGene Expression Omnibus(GEO)に送信できるようにする。

バルクおよびシングルセルRNA-Seq解析

 IRIS-EDAは遺伝子発現データ解析のための包括的なプラットフォームを提供するように設計された。シングルセルRNA-Seq(scRNA-Seq)データ解析はRNA-Seq解析の中で研究の成長分野であり、細胞の変異を考慮して遺伝子発現パターンに独自の洞察を提供することができる[ref,31、32]。従来のDGE分析に使用された方法は、適切なフィルタリングおよびDGE法と組み合わせた場合、scRNA-Seq DGE分析への適用性を実証した[ref.32]。したがって、IRIS-EDAは、わずかな修正を加えるだけで、scRNA-Seqデータの発見主導型およびDGE分析を容易にすることができる。すなわち、シングルセルデータの分析は、特にedgeRまたはlimmaのいずれかと組み合わせた場合に、100万回あたりの転写産物(TPM)> 1のデフォルト設定に基づく厳密なフィルターカットオフを使用することによって適切に実施できる。 10Xgenomics シングルセルデータのような全体的に低い発現レベルを期待する特定の種類のscRNA-Seqデータについては、違いを説明するために異なるアプローチが提供されている。特に、DESeq2正規化法は、最も信頼性の高い分析結果を提供するために、遺伝子のフィルタリングなしで使用される[ref.33]。 scRNA-Seqデータの分析に関する詳細は、S1 TextのSingle-cell RNA-Seqセクションにある。

必要な入力

 IRIS-EDAでは、使用するデータの種類に応じて、2つまたは3つのユーザー提供の入力ファイルが必要になる。(1)遺伝子発現推定マトリックス(EEM、リードカウントデータとも呼ばれる)、(2)因子レベルを含む条件マトリックスEEMの提供されたサンプルに対応する、および(3)scRNA-Seqデータのみのフィルタリングに使用される各遺伝子の塩基対の長さを示す遺伝子長マトリックス。データをアップロードするとき、ユーザーは遺伝子発現データタイプを選択する:バルクまたはシングルセルRNA-Seqデータ。 scRNA-Seqデータを使用している場合は、遺伝子長マトリックスの追加要件がWebサーバーに表示される。また、シングルセルデータの分析を最適化するためのデフォルトのパラメータ設定が、サーバー全体に表示される。 GFF / GTF / GFF3アノテーションファイルから遺伝子長を取得する方法はS1 TextのSingle-cell RNA-Seqセクションにある。

必要な入力を送信した後、3つの正規化アプローチのうちの1つを選択したり、ユーザーがデータを正規化しないことを選択することもできる。 IRIS-EDAで利用可能な3つの正規化方法は、正規対数変換、正規化対数変換、および分散安定化変換である。通常の対数変換では、各遺伝子の発現を正規化するために2を底とする対数関数を使用する。そうすることで、特に、多数のゼロが変換されていないプロットから収集された情報をほとんどもたらさない可能性があるスパース表現行列に対して、表現分布の視覚化を改善した。正則化された対数変換は、遺伝子数が少ないサンプル間の差異を最小限に抑え、ライブラリのサイズに基づいて正則化する方法を提供する[ref.15]。正規化された対数変換法は、ライブラリのサイズがそれほど変わらないデータセットに最も役立つ。分散安定化変換はライブラリーのサイズによっても正規化され、ほぼホモロジー的な表現行列を提供する[ref.15]。ライブラリサイズが大きく異なるデータセットの場合は、分散安定化変換方法が最も適している。

 

tutorial

http://bmbl.sdstate.edu/IRIS/

FAQ

http://bmbl.sdstate.edu/IRIS/

Github

R shinnyのアプリとして構築されており、localでも使用可能です。

 

使い方

http://bmbl.sdstate.edu/IRIS/にアクセスする。

f:id:kazumaxneo:20190422002037p:plain

 

Count dataは以下のようなフォーマットで準備する。

f:id:kazumaxneo:20190422002601p:plain

 

Metadataは以下のようなフォーマットで準備する。

f:id:kazumaxneo:20190422002654p:plain

 

localで使用することもできる。Rstudioでテストした。

# CRAN
packages <- c(
"crosstalk", "dplyr", "DT", "gtools", "plotly", "shiny", "plyr",
"shinyBS", "shinycssloaders", "shinythemes", "tibble", "tidyr",
"Rcpp", "Hmisc", "ggplot2", "locfit", "GGally", "pheatmap",
"reshape2", "backports", "digest", "fields", "psych", "stringr",
"tools", "openxlsx", "Rtsne", "WGCNA", "flashClust", "parallel",
"MCL", "kmed", "ape"
)
np <- packages[!(packages %in% installed.packages()[, "Package"])]
if(length(np)) install.packages(np)


# Bioconductor
bioc_packages <- c(
"DESeq2", "edgeR", "limma", "QUBIC", "geneplotter", "GO.db", "impute",
"preprocessCore", "AnnotationDbi"
)
np <- bioc_packages[!(bioc_packages %in% installed.packages()[,"Package"])]
if (!require("BiocManager")) install.packages("BiocManager")
BiocManager::install(np)

# GitHub
if (!require("devtools")) install.packages("devtools")
devtools::install_github("OSU-BMBL/BRIC", force = T)


shiny::runGitHub("iris", "OSU-BMBL")

立ち上げたところ。

f:id:kazumaxneo:20190423005923p:plain

 

引用

IRIS-EDA: An integrated RNA-Seq interpretation system for gene expression data analysis
Brandon Monier , Adam McDermaid , Cankun Wang, Jing Zhao, Allison Miller, Anne Fennell, Qin Ma

PLoS Comput Biol. 2019 Feb 14;15(2)

 

 

CheckMのplotコマンド

2021 4/2 追記

2021 8/30 checkm tetra 追記

 

checkmのゲノムアセンブリ評価コマンドについて以前紹介した。

ここでは、タイトルの通りCheckMのplotコマンドについて簡単に紹介する。このコマンドはbinningして得た一連のbinned.fastaディレクトリに対して実行でき、ラフに各binを評価することができる。

各binごとにplotしたファイルが出力される。

f:id:kazumaxneo:20191130193352p:plain

 500binくらいあっても1時間ほどで結果は出力される。そんなに待たされない。

 

wiki

CheckM wiki - plots

help

> checkm gc_plot -h

$ checkm gc_plot -h

usage: checkm gc_plot [-h] [--image_type {eps,pdf,png,ps,svg}] [--dpi DPI]

                      [--font_size FONT_SIZE] [-x EXTENSION] [--width WIDTH]

                      [--height HEIGHT] [-w GC_WINDOW_SIZE] [-b GC_BIN_WIDTH]

                      [-q]

                      bin_dir output_dir dist_value [dist_value ...]

 

Create GC histogram and delta-GC plot.

 

positional arguments:

  bin_dir               directory containing bins to plot (fasta format)

  output_dir            directory to hold plots

  dist_value            reference distribution(s) to plot; integer between 0 and 100

 

optional arguments:

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

  --image_type {eps,pdf,png,ps,svg}

                        desired image type (default: png)

  --dpi DPI             desired DPI of output image (default: 600)

  --font_size FONT_SIZE

                        Desired font size (default: 8)

  -x, --extension EXTENSION

                        extension of bins (other files in directory are ignored) (default: fna)

  --width WIDTH         width of output image (default: 6.5)

  --height HEIGHT       height of output image (default: 3.5)

  -w, --gc_window_size GC_WINDOW_SIZE

                        window size used to calculate GC histogram (default: 5000)

  -b, --gc_bin_width GC_BIN_WIDTH

                        width of GC bars in histogram (default: 0.01)

  -q, --quiet           suppress console output

 

checkm tetra_plot -h

$ checkm tetra_plot -h

usage: checkm tetra_plot [-h] [--image_type {eps,pdf,png,ps,svg}] [--dpi DPI]

                         [--font_size FONT_SIZE] [-x EXTENSION]

                         [--width WIDTH] [--height HEIGHT] [-w TD_WINDOW_SIZE]

                         [-b TD_BIN_WIDTH] [-q]

                         results_dir bin_dir output_dir tetra_profile

                         dist_value [dist_value ...]

 

Create tetranucleotide distance (TD) histogram and delta-TD plot.

 

positional arguments:

  results_dir           directory specified during qa command

  bin_dir               directory containing bins to plot (fasta format)

  output_dir            directory to hold plots

  tetra_profile         tetranucleotide profiles for each bin (see tetra command)

  dist_value            reference distribution(s) to plot; integer between 0 and 100

 

optional arguments:

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

  --image_type {eps,pdf,png,ps,svg}

                        desired image type (default: png)

  --dpi DPI             desired DPI of output image (default: 600)

  --font_size FONT_SIZE

                        Desired font size (default: 8)

  -x, --extension EXTENSION

                        extension of bins (other files in directory are ignored) (default: fna)

  --width WIDTH         width of output image (default: 6.5)

  --height HEIGHT       height of output image (default: 3.5)

  -w, --td_window_size TD_WINDOW_SIZE

                        window size used to calculate TD histogram (default: 5000)

  -b, --td_bin_width TD_BIN_WIDTH

                        width of TD bars in histogram (default: 0.01)

  -q, --quiet           suppress console output

 

Example: checkm tetra_plot ./output ./bins ./plots tetra.tsv 95

checkm len_plot -h

$ checkm len_plot -h

usage: checkm len_plot [-h] [--image_type {eps,pdf,png,ps,svg}] [--dpi DPI]

                       [--font_size FONT_SIZE] [-x EXTENSION] [--width WIDTH]

                       [--height HEIGHT] [-q]

                       bin_dir output_dir

 

Cumulative sequence length plot.

 

positional arguments:

  bin_dir               directory containing bins to plot (fasta format)

  output_dir            directory to hold plots

 

optional arguments:

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

  --image_type {eps,pdf,png,ps,svg}

                        desired image type (default: png)

  --dpi DPI             desired DPI of output image (default: 600)

  --font_size FONT_SIZE

                        Desired font size (default: 8)

  -x, --extension EXTENSION

                        extension of bins (other files in directory are ignored) (default: fna)

  --width WIDTH         width of output image (default: 6.5)

  --height HEIGHT       height of output image (default: 6.5)

  -q, --quiet           suppress console output

 

Example: checkm len_plot ./bins ./plots

checkm marker_plot -h

$ checkm marker_plot -h

usage: checkm marker_plot [-h] [--image_type {eps,pdf,png,ps,svg}] [--dpi DPI]

                          [--font_size FONT_SIZE] [-x EXTENSION]

                          [--width WIDTH] [--height HEIGHT]

                          [--fig_padding FIG_PADDING] [-q]

                          results_dir bin_dir output_dir

 

Plot position of marker genes on sequences.

 

positional arguments:

  results_dir           directory specified during qa command

  bin_dir               directory containing bins to plot (fasta format)

  output_dir            directory to hold plots

 

optional arguments:

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

  --image_type {eps,pdf,png,ps,svg}

                        desired image type (default: png)

  --dpi DPI             desired DPI of output image (default: 600)

  --font_size FONT_SIZE

                        Desired font size (default: 8)

  -x, --extension EXTENSION

                        extension of bins (other files in directory are ignored) (default: fna)

  --width WIDTH         width of output image (default: 6.5)

  --height HEIGHT       height of output image (default: 6.5)

  --fig_padding FIG_PADDING

                        white space to place around figure (in inches) (default: 0.2)

  -q, --quiet           suppress console output

 

Example: checkm marker_plot ./output ./bins ./plots

checkm par_plot -h

$ checkm par_plot -h

usage: checkm par_plot [-h] [--image_type {eps,pdf,png,ps,svg}] [--dpi DPI]

                       [--font_size FONT_SIZE] [-x EXTENSION] [--width WIDTH]

                       [--height HEIGHT] [-q]

                       results_dir bin_dir output_dir coverage_file

 

Parallel coordinate plot of GC and coverage.

 

positional arguments:

  results_dir           directory specified during qa command

  bin_dir               directory containing bins to plot (fasta format)

  output_dir            directory to hold plots

  coverage_file         file indicating coverage of each sequence (see coverage command)

 

optional arguments:

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

  --image_type {eps,pdf,png,ps,svg}

                        desired image type (default: png)

  --dpi DPI             desired DPI of output image (default: 600)

  --font_size FONT_SIZE

                        Desired font size (default: 8)

  -x, --extension EXTENSION

                        extension of bins (other files in directory are ignored) (default: fna)

  --width WIDTH         width of output image (default: 6.5)

  --height HEIGHT       height of output image (default: 6.5)

  -q, --quiet           suppress console output

 

Example: checkm par_plot ./output ./bins ./plots coverage.tsv

checkm cov_pca -h

$ checkm cov_pca -h

usage: checkm cov_pca [-h] [--image_type {eps,pdf,png,ps,svg}] [--dpi DPI]

                      [--font_size FONT_SIZE] [-x EXTENSION] [--width WIDTH]

                      [--height HEIGHT] [-q]

                      bin_dir output_dir coverage_file

 

PCA plot of coverage profiles.

 

positional arguments:

  bin_dir               directory containing bins to plot (fasta format)

  output_dir            directory to hold plots

  coverage_file         file indicating coverage of each sequence (see coverage command)

 

optional arguments:

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

  --image_type {eps,pdf,png,ps,svg}

                        desired image type (default: png)

  --dpi DPI             desired DPI of output image (default: 600)

  --font_size FONT_SIZE

                        Desired font size (default: 8)

  -x, --extension EXTENSION

                        extension of bins (other files in directory are ignored) (default: fna)

  --width WIDTH         width of output image (default: 6.5)

  --height HEIGHT       height of output image (default: 6.5)

  -q, --quiet           suppress console output

 

Example: checkm cov_pca ./bins ./plots coverate.tsv

checkm coverage -h

$ checkm coverage -h

usage: checkm coverage [-h] [-x EXTENSION] [-r] [-a MIN_ALIGN]

                       [-e MAX_EDIT_DIST] [-m MIN_QC] [-t THREADS] [-q]

                       bin_dir output_file bam_files [bam_files ...]

 

Calculate coverage of sequences.

 

positional arguments:

  bin_dir               directory containing bins (fasta format)

  output_file           print results to file

  bam_files             BAM files to parse

 

optional arguments:

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

  -x, --extension EXTENSION

                        extension of bins (other files in directory are ignored) (default: fna)

  -r, --all_reads       use all reads to estimate coverage instead of just those in proper pairs

  -a, --min_align MIN_ALIGN

                        minimum alignment length as percentage of read length (default: 0.98)

  -e, --max_edit_dist MAX_EDIT_DIST

                        maximum edit distance as percentage of read length (default: 0.02)

  -m, --min_qc MIN_QC   minimum quality score (in phred) (default: 15)

  -t, --threads THREADS

                        number of threads (default: 1)

  -q, --quiet           suppress console output

 

Example: checkm coverage ./bins coverage.tsv example_1.bam example_2.bam

 

 

準備

 一部のコマンドはcheckm出力を要求するので、はじめにbinned.fasta群に対してcheckmの標準コマンドを実行しておく。

checkm lineage_wf -t 20 -x fa metagenome/ output 1> log

 

実行方法

gc_plot

ゲノムビン内の配列のGC分布の評価に適したグラフの出力。binned_dir はメタゲノムのコンティグのディレクトリ(ここではmetagenome/)、output_dir は出力ディレクトリの指定。

checkm gc_plot -x fa binned_dir output_dir 95

f:id:kazumaxneo:20191130144902p:plain
左のグラフはGCヒストグラム。典型的な 純化されたゲノムはシングルピークの分布になる。右のグラフは、ゲノムビンの各シーケンスを、ゲノム全体の平均GC(x軸)およびシーケンス長(y軸)からの偏差の関数としてプロットしたもの。赤の破線は長さの関数として平均GCから予想される偏差を表す。
 

 

coding_plot

ゲノムビン内の配列のcoding density評価に適したグラフの出力。checkm_result_dir は

checkmのcheckm lineage_wfコマンドの出力ディレクトリ。binned_dir はメタゲノムのコンティグのディレクトリ(ここではmetagenome/)、output_dir は出力ディレクトリの指定。

checkm coding_plot -x fa checkm_result_dir binned_dir output_dir 95

f:id:kazumaxneo:20191130145135p:plain


 

 len_plot

ゲノムビンの累積シーケンス長プロット

checkm len_plot -x fa binned_dir output_dir

f:id:kazumaxneo:20191130145350p:plain

 

marker_plot

ゲノムビン配列上のマーカー遺伝子の位置をプロット。 これにより、マーカー遺伝子が連結されている範囲に関する情報が提供される。  マーカー遺伝子のない配列は表示されない。

checkm marker_plot -x fa checkm_result_dir binned_dir output_dir

f:id:kazumaxneo:20191130145442p:plain

 

par_plot

ゲノムビン内の各配列のGCカバレッジを示す平行座標プロットを生成する。 典型的なゲノムでは、すべてのシーケンスがプロット全体に同様のパスを生成する。 配列の分岐パスを持つシーケンスは汚染である可能性がある。このプロットには、ゲノムビン内のすべての配列のカバレッジファイルが必要になる。前もってcoverageコマンドを実行しておく。

#coverage計算
checkm coverage -x fa binned_dir output_coverage input.bam

#(optional)
checkm profile output_coverage > coverage_profile

#par_plot、coverageファイルを指定する。
checkm par_plot -x fa checkm_result binned_dir output_dir output_coverage

f:id:kazumaxneo:20191130145621p:plain

 

cov_pca

配列間のカバレッジプロファイル距離の主成分プロット(PCA)を生成する。 このプロットには、ゲノムビン内のすべての配列のカバレッジプロファイルを示すファイルが必要( 3サンプル以上必要)。

#coverage計算
checkm coverage -x fa binned_dir output_coverage \
sample1.bam sample2.bam sample3.bam

#par_plot、coverageファイルを指定する。
checkm par_plot -x fa checkm_result binned_dir output_dir output_coverage

 

2021 8/30

tetra_plot

checkm tetra -t 30 assembly.fasta tetranucleotide.tsv
checkm tetra_plot -x fasta checkM_outdir/ bin_dir/ output_dir tetranucleotide.tsv 95

f:id:kazumaxneo:20210830204217p:plain

 

 

他にもいくつかのplotコマンドがあります。上のリンク先から確認して下さい。 

引用

CheckM: assessing the quality of microbial genomes recovered from isolates, single cells, and metagenomes.

Parks DH, Imelfort M, Skennerton CT, Hugenholtz P, Tyson GW.

Genome Res. 2015 Jul;25(7):1043-55.

 

 

 

PHYLUCE

 

 保存された領域、または超保存 (ultraconserved) された領域(以下、保存された遺伝子座 (conserved loci) )のエンリッチメントは、非モデル生物(Faircloth et al、2012、2013、2015)の複数の時間スケールでの普遍的なphylogenomic analysesを可能にする(Faircloth et al、2012; Smith et al、 2014)。このアプローチの強みは何千もの種から何千もの遺伝子座の配列データを収集する能力から得られ、Classes(>200–300 Ma)のような系統発生の巨大な裂け目を超えた系統比較から、populations(<0.5–5 Ma)のようなより小さい関係の進化的多様性のような系統比較を可能にすることにある。種の進化の歴史を推測することがデータ収集の目的である場合、その後の分析作業は一般的に次のようなものになる。(i)シーケンシングリードをアセンブリする。これは数十から数百の個体に及ぶかもしれない (ii)アセンブリされたサンプルごとにコンティグ間の推定オルソログを同定する。一方で推定されたパラログを除去する (iii)異なる個体、他の実験から含まれる個体または個々のゲノム配列から採取された保存された遺伝子座データを含むデータセットを容易に生成する; (iv)全セットに渡り、オルソログ由来のシーケンシングリードを同定し、exportする; (v)データをアラインメントさせ、場合によっては余分な領域をトリミングし、phylogenetic inferenceに備える。 (vi)アラインメントされたデータに関する要約統計を計算し、(vii)配列またはアラインメントデータにユーティリティ機能を実行して、種々の系統発生推論プログラムを用いた下流分析へ向けた準備をする。 PHYLUCE(pronounced ‘phy-loo-chee’)は、これらのタスクを保存されエンリッチされた遺伝子座に対して実行する、計算効率が良くインストールしやすい、初めてのオープンソースのソフトウェアパッケージである。

 

マニュアル

https://phyluce.readthedocs.io/en/latest/installation.html

 

 PHYLUCEに関するツイート

 

インストール

ubuntu16.04のanacondaminiconda2-4.3.30に導入した。

build and test the binaries available through conda using 64-bit operating systems that include:

依存

3rd-party dependencies and packages installed

- python
- abyss 1.5.2
- bcftools
- bedtools
- biopython
- bwa
- bx-python
- dendropy 3.12.3
- gblocks
- lastz
- mafft
- muscle
- pandas
- picard
- pysam
- pyvcf
- raxml
- samtools
- seqtk
- trimal
- trinity # [not osx]
- velvet
- illumiprocessor
- spades
- itero

本体 Github

#Bioconda (link)
conda create -n phyluce -c bioconda -y phyluce

> phyluce_assembly_get_fastq_lengths -h

$ phyluce_assembly_get_fastq_lengths -h

usage: phyluce_assembly_get_fastq_lengths [-h] --input INPUT [--csv]

 

Get summary (length) data from fastq

 

optional arguments:

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

  --input INPUT  The directory of fastq files to summarize

  --csv          Give output in CSV

 

テストラン

Tutorial I: UCE Phylogenomics

1、サンプルの準備

mkdir uce-tutorial
cd uce-tutorial/
wget -O fastq.zip https://ndownloader.figshare.com/articles/1284521/versions/1

mkdir raw-fastq
mv fastq.zip raw-fastq
cd raw-fastq
unzip fastq.zip && rm fastq.zip
cd ../

 

2、リードQT

uce-tutorial/に以下の内容のconfigファイルを作成

cat >illumiprocessor.conf <<EOF 
# this is the section where you list the adapters you used. the asterisk
# will be replaced with the appropriate index for the sample.
[adapters]
i7:AGATCGGAAGAGCACACGTCTGAACTCCAGTCAC*ATCTCGTATGCCGTCTTCTGCTTG
i5:AGATCGGAAGAGCGTCGTGTAGGGAAAGAGTGTAGATCTCGGTGGTCGCCGTATCATT

# this is the list of indexes we used
[tag sequences]
BFIDT-166:GGAGCTATGG
BFIDT-016:GGCGAAGGTT
BFIDT-045:TTCTCCTTCA
BFIDT-011:CTACAACGGC

# this is how each index maps to each set of reads
[tag map]
Alligator_mississippiensis_GGAGCTATGG:BFIDT-166
Anolis_carolinensis_GGCGAAGGTT:BFIDT-016
Gallus_gallus_TTCTCCTTCA:BFIDT-045
Mus_musculus_CTACAACGGC:BFIDT-011

# we want to rename our read files something a bit more nice - so we will
# rename Alligator_mississippiensis_GGAGCTATGG to alligator_mississippiensis
[names]
Alligator_mississippiensis_GGAGCTATGG:alligator_mississippiensis
Anolis_carolinensis_GGCGAAGGTT:anolis_carolinensis
Gallus_gallus_TTCTCCTTCA:gallus_gallus
Mus_musculus_CTACAACGGC:mus_musculus

EOF

 

3、データのクリーニング。ここではillumiprocessor (github) が使われている(trimmomaticの並列処理ツール)。

illumiprocessor --input raw-fastq/ --output clean-fastq \
--config illumiprocessor.conf\ --cores 4

アダプターとlow qualityリードが除かれたシーケンシングデータはサブディレクトリのsplit-adapter-quality-trimmed/に出力される。

$ ls -alh clean-fastq/gallus_gallus/split-adapter-quality-trimmed/

total 100M

drwxr-xr-x 5 kazu kazu  170 Oct  8 03:33 .

drwxr-xr-x 6 kazu kazu  204 Oct  8 03:31 ..

-rw-r--r-- 1 kazu kazu 264K Oct  8 03:33 gallus_gallus-READ-singleton.fastq.gz

-rw-r--r-- 1 kazu kazu  49M Oct  8 03:33 gallus_gallus-READ1.fastq.gz

-rw-r--r-- 1 kazu kazu  52M Oct  8 03:33 gallus_gallus-READ2.fastq.gz

READ1とREAD2¥はペアエンドの順番が維持されているファイル。singletonはペアエンドの片方がQTで除かれてシングルになったファイル。

 

 4、アセンブリ。データのパスを示した以下のconfigファイルを作成。ここでは/data/uce-tutorial/で作業したので以下のようになっている。

cat >assembly.conf <<EOF 
[samples]
alligator_mississippiensis:/data/uce-tutorial/clean-fastq/alligator_mississippiensis/split-adapter-quality-trimmed/
anolis_carolinensis:/data/uce-tutorial/clean-fastq/anolis_carolinensis/split-adapter-quality-trimmed/
gallus_gallus:/data/uce-tutorial/clean-fastq/gallus_gallus/split-adapter-quality-trimmed/
mus_musculus:/data/uce-tutorial/clean-fastq/mus_musculus/split-adapter-quality-trimmed/
EOF

#de novo assembly
phyluce_assembly_assemblo_trinity \
--conf assembly.conf \
--output trinity-assemblies \
--clean \
--cores 12

エラーが出て終了する。 修正できたら追記します。

  

引用
PHYLUCE is a software package for the analysis of conserved genomic loci
Faircloth BC

Bioinformatics. 2016 Mar 1;32(5):786-8

 

ロングリードのマッピングから逆位を検出する npInv

 

 DNAのセグメントの向きが、染色体の残りの部分と比較してその先祖から反転している逆位多型(Inversion polymorphisms)は、ショウジョウバエの異なる系統のハイブリッドにおける染色体間の組換えの抑制因子として、スターテバントによって1917年に最初に発見された[ ref.1]。逆位は、非相同末端結合(NHEJ [ref.2])、非対立遺伝子相同組換え(NAHR)、またはフォークストールとテンプレートスイッチング(FoSTeS [逆位3])に大まかに分類できる。 NHEJは、DNAの二本鎖切断を修復するための経路である。逆位配列は、大きな相同配列なしでブレークポイントに直接ライゲーションする[ref.2]。 NAHRは、相同配列間で発生する異常な組換えメカニズムである。逆方向反復(IR)間の相同組換えは、介在配列を反転させ、逆位を作る[ref.4]。ほとんどすべて(12/14)の既知の大規模な逆位(> 1 Mb)多型はNAHRによって媒介される[ref.5]。( 一部略)

 NpInvというロングリードシーケンスデータからNAHRを介した逆位を検出およびジェノタイピングする新しいツールを提示する。 npInvへの入力は、BWA-MEM [ref.12]などのローカルアライナーから生成されたbam形式のアライメントファイルである。 npInvのパイプラインと擬似コードは、それぞれ論文図1と補足メソッドに示されている。簡単に言えば、npInvはアラインメントファイルをスキャンして、同じ染色体にマッピングされているが方向が異なるサブリードアラインメントのペアを含むリードを探す(論文図2)。 npInvは、このサブリードアラインメントペアを逆位のシグナルとして記録する。npInvは、クラスター内の逆位シグナルの位置と数に基づいて逆位イベントを検出するために、すべての逆位シグナルをクラスター化およびフィルター処理する。 npInvは、逆位をサポートするリード数と、非逆位対立遺伝子をサポートするリード(逆位ブレークポイントにまたがるリード)の両方を報告する。

 

インストール

ubuntu18.04 LTSのjava1.8環境でテストした。

依存

Java version 1.8+.

本体 Github

 リリースからビルドされたnpInv1.24.jarをダウンロードできる。

java -jar npInv1.24.jar -h

$ java -jar npInv1.24.jar -h

No Key Name: output

Program function: Read a SE bam file and get the inversion

Version: 1.24

--output[String] file to write

--input[String] file to read

optional:

--region[String] Specify the region for running.

                 Such as chr9:1-1000 OR chr9 OR all Default[all]

--minAln[int] minimum size for Alignment & Inv. Default[500]

--IRdatabase[String] An inverted repeat file for the reference in bed format. Default[none]

--min[int] minimum size of an inversion. Default[500]

--max[int] maximum size of an inversion. Default[10000]

--window[int] minimun window size (bp) to merge inversion breakpoints. Default[2000]

--threshold[int] minimum number of supporting reads for an inversion. Default[3]

--help Show usage

For example: java -jar npInv.jar --input sample.bam --output sample.VCF

 

 

実行方法

bamを指定する。

java -jar npInv.jar --input input.bam --output out.vcf

 

 

引用

npInv: accurate detection and genotyping of inversions using long read sub-alignment

Haojing Shao, Devika Ganesamoorthy, Tania Duarte, Minh Duc Cao, Clive J. Hoggart, Lachlan J. M. Coin
BMC Bioinformatics volume 19, Article number: 261 (2018)

 

関連


SolidBin

 

 微生物群集は多数の微生物のmixtureからなる。メタゲノミクス研究では、遺伝物質は微生物群集から直接抽出され、ハイスループットシークエンシング技術を使用してシークエンスされ、異なるリード長の大量のシークエンスリードが得られる。
 通常、メタゲノムデータ解析の最初のステップは、リードのオーバーラップ情報に基づいてリードをコンティグと呼ばれる比較的長いゲノムフラグメントにアセンブリする「リードアセンブリ」で、次に、同じビンにゲノムまたはcloselt relatedなゲノムを仕分ける。 3番目のステップでは、リードを既知のゲノム配列またはビン内のコンティグにマッピングし、ゲノムとゲノムビンの相対存在量プロファイルを各サンプルで推定する。最後に、統計学的方法を用いて、ゲノム/ビンの存在量プロファイルを個体の環境因子または表現型と関連付ける。メタゲノム研究は、ヒトの腸内の微生物群集がいくつかの疾患(Khor et al、2011; Jostins et al、2012)および免疫療法の効果(Wilck et al、2017; Chen et al)と有意に関連していることを示した。本稿では、メタゲノムデータ解析の通常の第2段階であるコンティグビニングに焦点を当てる。
 利用可能なコンティグビニング方法は、2つの大きなカテゴリーに分類することができる:MEGAN(Husons et al, 2007)のようなコンティグを既知のリファレンスデータベースにマッピングするtaxonomy依存の方法(教師あり方法とも呼ばれる)、MaxBin 、CONCOCT、MetaBAT、COCACOLAおよびMBC3Cなどの教師なし方法。これらの方法にはそれぞれ独自の限界がある。既知のゲノムは限られており、多数のゲノムが未だに未知であるため、多くの未知コンティグが存在する。さらに、アラインメントは計算上高価であり、したがってアラインメントに基づくコンティニングビニング方法は一般に遅い。教師なしビニング方法は上記の欠点のいくつかを部分的に克服することができる。既存の教師なしビニング方法は、いくつかのコンティグの既知のゲノムへのアラインメントなど、利用可能な追加の生物学的情報を十分に活用していないため、ビニング性能が向上する可能性がある。したがって、未知のゲノムからの大量のコンティグを含むデータセットに使用できるだけでなく、より良いビニング性能のためにアラインメントから利用可能な生物学的情報を利用することができる効果的なコンティグビニングアプローチを開発することが不可欠である。

 この問題に取り組むために、本著者らは、ペアワイズ組み入れによってクラスタリング性能を改善することができる様々な半教師クラスタリング研究(Wagstaff et al、2001; Gu et al、2013)に触発され、新規コンティグビニング方法SolidBinを開発した。

 

 

f:id:kazumaxneo:20190421123004j:plain

The general workflow of SolidBin.  論文より転載。

 

インストール

本体 Github

python3.6の仮想環境を作って導入することが推奨されている。

git clone https://github.com/sufforest/SolidBin
cd SolidBin/
conda create -n solidbin -y python=3.6 numpy pandas scikit-learn scipy
conda activate solidbin

>python SolidBin.py -h

$ python SolidBin.py -h

usage: SolidBin.py [-h] [--contig_file CONTIG_FILE]

                   [--coverage_profiles COVERAGE_PROFILES]

                   [--composition_profiles COMPOSITION_PROFILES]

                   [--priori_ml_list PRIORI_ML_LIST]

                   [--priori_cl_list PRIORI_CL_LIST] [--output OUTPUT]

                   [--log LOG] [--clusters CLUSTERS] [-a A] [-b B] [--use_sfs]

 

optional arguments:

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

  --contig_file CONTIG_FILE

                        The contigs file.

  --coverage_profiles COVERAGE_PROFILES

                        The coverage profiles, containing a table where each

                        row correspond to a contig, and each column correspond

                        to a sample. All values are separated with tabs.

  --composition_profiles COMPOSITION_PROFILES

                        The composition profiles, containing a table where

                        each row correspond to a contig, and each column

                        correspond to the kmer composition of particular kmer.

                        All values are separated with comma.

  --priori_ml_list PRIORI_ML_LIST

                        The edges encoding either the co-alignment or other

                        priori information, one row of an edge is in the

                        format: contig_name_A contig_name_B 1. The edge is

                        undirected.

  --priori_cl_list PRIORI_CL_LIST

                        The cannot-link information, one row of an edge is in

                        the format: contig_name_A contig_name_B 1. The edge is

                        undirected.

  --output OUTPUT       The output file, storing the binning result.

  --log LOG             Specify where to store log file

  --clusters CLUSTERS   Specify the number of clusters. If not specified, the

                        cluster number is estimated by single-copy genes.

  -a A                  weight of must-link constraints

  -b B                  weight of cannot-link constraints

  --use_sfs             use SFS as constraints

 

docker imageも 公開されている。

docker run -it -v /data/StrainMock/input:/input -v /data/StrainMock/output:/output solidbin python SolidBin.py --contig_file /input/StrainMock_Contigs_cutup_10K_nodup_filter_1K.fasta --composition_profiles /input/kmer_4.csv --coverage_profiles /input/cov_inputtableR.tsv --output /output/result.tsv --log /output/log.txt

 

 

実行方法 

 

 

引用

SolidBin: Improving Metagenome Binning with Semi-supervised Normalized Cut
Ziye Wang, Zhengyang Wang, Yang Young Lu, Fengzhu Sun, Shanfeng Zhu
Bioinformatics, Published: 12 April 2019