macでインフォマティクス

macでインフォマティクス

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

メタゲノムの組成を解析する CAMAMED

 

 

 メタゲノミクスは、分子ゲノミクス、微生物生態学、データ解析が交差する学際的な研究分野である。この分野の主な研究対象は、ある環境に存在する微生物のゲノム総量を指すメタゲノムである。メタゲノミクスは、ハイスループットゲノムシークエンシング技術をはじめとする微生物の培養に依存しない方法に基づいている。この方法では、腸などの特定の環境に生息する微生物(マイクロバイオーム)のDNAを抽出し、様々な計算機技術を用いて研究する。

 メタゲノムデータを解析するためのアプローチは大きく分けて2つある。(i) サンプル中の微生物の系統的多様性を記述する分類学的プロファイリング、(ii) サンプル中の微生物の機能的能力についての洞察を得るために、ゲノム配列を「機能的」グループにマッピングするための計算戦略を含む機能的プロファイリングである。

 マイクロバイオームの包括的な機能解析は、微生物群集の生化学的能力に関する我々の理解を大幅に向上させることができる(ref.3)。マイクロバイオームの機能性プロファイリングのための戦略として、リードをアセンブリしてより大きなコンティグにし、機能性プロファイリングのために遺伝子量を予測することが考えられる(ref.4)。これまで研究されていないマイクロバイオームに対しては、従来のアセンブリのほとんどが比較的精度が低いという事実にもかかわらず、アセンブリに基づく戦略は避けられない(ref.5)。対照的に、ヒトの腸のように既に広範囲に研究されている環境では、良好なマイクロバイオーム遺伝子カタログが以前に編集されている(ref.6)。このような場合には、マッピングに基づく解析をメタゲノムの機能解析に利用することができる。

 ショットガンメタゲノミクス研究では、微生物のコレクションは、培養や分離を行わずに直接DNAシークエンシングを行うことで研究される。遺伝子カタログにマッピングされた遺伝子の頻度を比較することにより、サンプル間の機能的差異を研究することができる。このように、遺伝子の豊富さデータは高いレベルの系統的変動の影響を受け、統計的な力を大きく低下させ、偽陽性を増加させる可能性がある(ref.7)。メタゲノムデータの系統的変動が遺伝子や微生物の観察された豊富さに影響を与える理由は多くある。重要な理由の一つは、各サンプルが異なる数のシークエンシングリードを持つことである。つまりシークエンシングの深さが異なることである(ref.8)。系統的変動の他の理由としては、サンプリング方法の不一致、DNA抽出、シーケンシングランの品質のばらつき、リードマッピングのエラー、リファレンス遺伝子カタログの不完全さなどが挙げられる(ref.9)。さらに、系統的変動は、微生物の平均ゲノムサイズの違い、種の豊富さ、リードに関連するGC含有量の違いに起因することがあり、これが観察された遺伝子の豊富さに影響を与える可能性がある((ref.7,10)。

 次世代シークエンシング(NGS)データもまた、本質的に組成的(compositional)である。組成(compositional)とは、各ヌクレオチド断片の相対的な豊富さが他の断片の豊富さに依存していることを意味する。この性質は、シーケンシングマシンおよび基礎となる方法に関係しており、結果として得られる配列は、増幅およびそれに続くヌクレオチドサンプリングに関与するバイアスの影響を受ける((ref.11)。したがって、組成は、全体の不明確な部分である測定(例えば、NGSシーケンシングによって生成されたメタゲノムカウントデータ)において、このような曖昧さの結果である。The compositional data analysis (CoDA) は、このバイアスを処理して解決することを目的としている((ref.12)。

 メタゲノムカウントデータもまた、他のNGSデータと比較して、より厳しい課題に直面している。これらの課題の1つは、異なるサンプルにおけるシークエンシングされたリードの数やシークエンシングの深さが非常に変動しやすいことである。第二の課題は、ゼロインフレーションと呼ばれるメタゲノムカウントデータにおけるゼロの割合が非常に高いことである(約50~90%)(ref.13)。また、メタゲノムデータは、他のNGSデータと比較して非常に高次元である(例えば、腸内マイクロバイオーム遺伝子カタログのサンプルでは、約10Mの遺伝子配列やfeaturesがある(ref.6))。一方で、DNAサンプリングの頻度が低いため、非常に希少な分類群は記録されておらず、テクニカルゼロと呼ばれている。また、構造的ゼロと呼ばれる欠落した個体群が記録されていない場合もある。もう一つの課題は、研究のサイズと、分類群の分布の大きなばらつき(過剰分散)である(ref.14)。

 正規化処理は系統的な変動や組成の偏りを識別して除去することができるため、データの前処理や解析には欠かせないステップである。高次元カウントデータに対しては多くの正規化手法が提案されてきたが、そのほとんどはメタゲノムデータでの性能が評価されていない(ref.7)。組成データに関わる課題に対処するために、様々なアプローチがまだ提案されている。例えば、シーケンシングの深さが不均一であるという問題を解決するために、2つのアプローチが紹介されている。第一に、異なるサンプルのリード数を再スケーリングしてライブラリサイズの固定値を達成すること、第二に、すべてのサンプルについて固定リード数を達成するためにリードを再サンプリングすることである(ref.2)。また、多くのCoDA法では、正規化の代わりに変換を使用している。これらの方法では、対数比変換を用いてデータを実空間にマッピングするため、従来の統計的手法を用いて更なる解析を行うことが可能となる。これらの方法は、データ変換にサブセット特徴量の幾何平均などの参照を使用しようとしている(ref.15)。

 説明したように、メタゲノムの組成データにおける主要な課題の一つは、sparsity(まばら)またはゼロインフレーションであり、これは遺伝子中心のカウントデータではより深刻になる。この問題を処理し、比較遺伝子量研究の精度を向上させるために、いくつかのパッケージやツールが開発されてきた。いくつかのパッケージは、メタゲノムデータのゼロを扱うために特別に開発されてきた。metagenomeSeqでは、遺伝子量データにゼロインフレーション対数正規モデルを使用している(ref.16,17)。この方法では、ゼロインフレーションはサンプル固有のものであり、シークエンシングの深さに依存すると仮定している(ref.18)。比率アプローチ(Ratio approach for identifying differential abundance (RAIDA))は、最初にカウントを対数正規分布で記述される相対頻度に変換する統計モデルを使用している。RAIDA は、ほとんどの特徴は変動していないと仮定しているため、分類群や遺伝子レベルでのメタゲノムデータの解析に適している。また、このツールは 2 つの異なる条件でのメタゲノムデータサンプルの比較分析のために開発されたもので、2 つ以上の条件に一般化することが可能である(ref.19)。また、メタゲノムデータの統計モデルとしては、ゼロインフレーション負二項式やゼロインフレーションβ回帰モデルなど、いくつかのゼロインフレーション統計モデルがある(ref.20,21)。

 しかし、他にも、ANCOM(ref.22)、ZIBSeq(ref.21)、CPL(ref.23)、DESeq2(ref.24)、edgeR(ref.25)など、CoDAのためのパッケージやツールが提供されている。これらの手法は、異なる統計的手法を駆使し、ハイスループットシーケンシングデータの組成バイアスやゼロインフレーションを処理しようとしている。DESeq2とedgeRを含む上記のパッケージのほとんどは、もともとRNA-seqデータ解析のために開発されたものであることに注意することが重要である。最も広く使用されているCoDAパッケージのいくつかとその特性を論文表1に整理した。

 また、組成データの正規化手法も数多く提案されている。これらの手法のいくつかは、RNA-seqデータ、アンプリコンシーケンシングによって生成されたオペレーショナル・タクソノミック・ユニット(OTU)データ、およびメタゲノムデータに対して提供されている。しかし、これらの方法の比較研究では、性能とデータの特性の間に大きな依存性があることが示されている。例えば、RNA-seqデータに対してより優れた性能を持つ方法が、必ずしもメタゲノムデータに適しているとは限らない(ref.7)。これらの正規化手法のいくつかを以下に説明する。Total sum scaling(TSS)は、正規化された値の合計が1となるように、個々のカウントをサンプル中の総カウントで除算することで得られるカウントデータの標準的な正規化方法である(ref.28)。固定のスケーリング係数を持つ TSS は、技術的なシーケンシングのバイアスにより OTU カウントに悪影響を与える可能性がある。Cumulative sum scaling (累積和スケーリング)(CSS) (ref.16)は、低頻度(比較的一定で独立した)四分位に基づいてサンプルを再スケーリングし、高頻度サンプルの影響を排除しない。さらに、CSSは対数正規分布モデルとして高く評価されており、多くの研究で分類群/遺伝子レベルでのCoDAに使用されている(ref.2,18,29,30)。 Trimmed mean of M-values (TMM) は、差のある発現を識別するために、サンプル間のスケールファクターを推定する。TMM正規化は、ほとんどの遺伝子の発現がサンプル間で変動していないことを前提としている(ref.31)。組成データを正規化するための別の方法のファミリーも存在し、これは対数比変換に基づいている。組成における特徴間の相互依存性は、個々の特徴の分析が、各サンプルを新しい空間に変換するリファレンスベースラインに対して行われ、統計分析がこの新しい空間で行われることを暗示している。基準の選択に基づいて、異なる対数比に基づく方法が開発された。中心対数比(CLR)、加法対数比(ALR)、相対対数式(RLE)変換は、基準を選択するために異なる戦略を使用している(ref.11,32)。論文表1は、いくつかのパッケージ、関連する正規化方法、特性、およびCoDAのためのそれらの利点の詳細を示している。

 メタゲノムデータの機能解析のためには、メタゲノムデータの組成を考慮したパイプラインを使用する必要がある。本論文では、メタゲノムデータの分類学的・機能的解析を行うためのマッピングベースのソフトウェアパイプラインであるCAMAMEDを紹介する。そのため、適切な正規化手法は、メタゲノムデータの解析において非常に重要である。このデータの特性(例えば、スパース、高次元性、レアファクション、アンダーサンプリング、シークエンスデプスさの違いなど)を考慮すると、metagenomeSeq とその正規化手法としてのCSSは、メタゲノム研究において最も広く利用されている手法の一つである(ref.2,17,18,33)。しかし、METAGENOMESeqの最も重要な欠点の一つは、サンプルサイズが小さい場合のdifferential abundance analysis(存在量変動解析)の偽陽性率が高いことであると報告されている(ref.34,35)。CAMAMEDでは、分類、遺伝子、KO、EC数、反応レベルの疎なメタゲノムデータの組成バイアスに対応するため、metagenomeSeqを使用している。CAMAMED は、Linux オペレーティングシステムをベースとした Python3 と Shell プログラミングを用いて実装されている。CAMAMEDは非専門家向けに設計されており、比較的簡単に実行できるようになっている。また、CAMAMEDをより簡単に利用できるように、2つのDockerイメージを構築し、インストールの詳細や依存関係を介さずにCAMAMEDパイプラインを実行できるようにした。これらのイメージは www.hub.docker.com にある。
 この研究では、80のメタゲノムショットガン配列の糞便サンプルを使用した。このデータセット(当初は「cohort1」と命名)には、健常者(コントロール)24名、大腸腺腫27名、大腸がん29名からのサンプルが含まれている。これらのサンプルは、Illuminaプラットフォームとペアエンドシークエンシング法を用いてシークエンシングされた(36)。また、サンプル中の遺伝子の豊富さを得るために、9.88×106の非重複遺伝子配列を持つ既報の遺伝子カタログを用いた。このカタログは統合遺伝子カタログ(IGC)と呼ばれ、ヒト腸内マイクロバイオームのヌクレオチドとタンパク質の遺伝子配列が掲載されている(ref.6)。

 

マニュアルPDFはレポジトリに含まれています。

https://github.com/mhnb/CAMAMED/blob/master/CAMAMED-manual.pdf

 

インストール

docker imagesを使ってテストした。

依存

Software Requirements

  • Linux operating system (Preferably Ubuntu)
  • Python 2 or 3 (Preferably Python ≥3.7 and it is also better to have Python3 by default)

Hardware Requirements

  • It is better to run on a computer with at least 32GB of RAM and ten cores of processor.

Github

Docker

#dockerhub (link)
docker pull camamed/camamed_pipeline_db:latest

 > ./camamed_pre_processing.py

#  ./camamed_pre_processing.py

 

This function is a preprocessing step in pipeline

 

Command: ./camamed_pre_processing.py method [options]

 

      method:

  -h Shows help related to this function

  -sra  Running SAR toolkit to convert SRA files to Fastq or Fasta files.  (**optional**)

Copy SRA samples in the /CAMAMED/sra_files folder.

Copy the sample names in the /CAMAMED/sra_files/sra_file_names.txt.

for example:

file1.sra

file2.sra

file3.sra

After executing SRA-Toolkit, Fastq or Fasta files are automatically copied to folder /Read_files.

For Example: ./camamed_pre_processing.py -sra

 

  -gec  Get information on sequences and gene catalogs.

options:

  -gc (gene_catalog): Gene catalog sequences must have Fasta format.

        First copy the gene catalog into the /CAMAMED folder.

  -rff (read_files_format)[fastq/fasta]:

  Copy the sample files in the /CAMAMED/Read_files folder.

  Enter the file name of the samples in the /CAMAMED/Read_files/sample_file_names.txt

For single end fastq files:

    file1.fastq

    file2.fastq

    file3.fastq

For paired end fastq files:

    file1_1.fastq

    file1_2.fastq

    file2_1.fastq

    file2_2.fastq

  -rt (read_type): paired end or single end [p/s]

  -is (insert_size): For paired end sequences (An integer number) or ignor enter -1 (default=-1).

For Example: ./camamed_pre_processing.py -gec -gc gene_catalog.fa -rff fastq -rt p

    ./camamed_pre_processing.py -gec -gc gene_catalog.fa -rff fastq -rt p -is 350

 

  -cdh  Running CD-HIT on the gene catalog to remove redundant genes. (**optional**)

You can use this option if you think your gene catalog has redundant sequences.

options:

  -sit (sequence_identity_threshold): This value can be in the range of 0.8 to 1 (default=0.9).

For Example: ./camamed_pre_processing.py -cdh

    ./camamed_pre_processing.py -cdh -sit 0.95

 

  -pgc  Peprocessing gene catalog.

At this point, the names of the genes are deleted from the gene catalog and stored in the /CAMAMED/files/gene_name.txt

and for the genes the gene1, gene2 and ... are respectively selected.

For Example: ./camamed_pre_processing.py -pgc

> ./camamed_quality_control.py -h

# ./camamed_quality_control.py -h

 

This function executes fastqc quality control and SeqKit tools on sample sequences.

At this state, if the sequences have Fastq format, they will be quality controlled and their statistical information extracted.

But if they have Fasta format, only their statistical information obtained.

Before running this step, the sample data should be in the /CAMAMED/Read_files folder

and the file names should be written in the text file /CAMAMED/Read_files/sample_file_names.txt, respectively.

 

Command: ./camamed_quality_control.py method

 

      method:

  -h Shows help related to this function.

  -fqc  Execut FastQC quality control only for Fastq files.

For Example: ./camamed_quality_control.py -fqc

 

  -sek  Execut SeqKit to extract information from sample files.

For Example: ./camamed_quality_control.py -sek

 

  -eti  Extract total information from SeqKit outputs.

For Example: ./camamed_quality_control.py -eti

> ./camamed_metaphlan_profiling.py -h

# ./camamed_metaphlan_profiling.py -h

 

In this function, using the Metaphlan2 software extract the abundance of all bacteria.

Metaphlan2 can produce taxonomic profiling at different levels.

Such as Kingdom, Phylum, Class, Order, Family, Genus, and Species.

 

Command: ./camamed_metaphlan_profiling.py method [options]

 

      method:

  -h Shows help related to this function.

  -mph  Execute metaphlan2.

options:

    -tl (taxa_level): For selecting Kingdom, Phylum, Class, Order, Family, Genus, or Species,

    Select one of {'k', 'p', 'c', 'o', 'f', 'g', 's'}

    -c (core_number): select number of cores for Metaphlan2 execution (default=1).

    The number of cores must be a positive integer, otherwise one is chosen.

For Example: ./camamed_metaphlan_profiling.py -mph -tl s -c 3

 

  -imp  Extract information from metaphlan2 output files.

For Example: ./camamed_metaphlan_profiling.py -imp

>  ./camamed_kegg_annotation.py -h

# ./camamed_kegg_annotation.py -h

 

This function extract annotated information from KEGG databases

 

Command: ./camamed_kegg_annotation.py method [options]

 

      method:

  -h Shows help related to this function

  -pgc  Preparing the Gene Catalog for extracting annotated information from the KEGG databases.

To obtain KOs associated with gene sequences using web services GhostKOALA and KAAS,

it would be better if the file size of the gene catalog is less than 300 MB.

The gene catalog is stored after the preprocessing and deletion of the duplicate genes (optional)

named main_gene_catalog.fa in the CAMAMED folder of the program.

If the gene catalog size is more than 300MB, you can use this option to convert it to smaller files.

For example, if the gene catalog has 300 gene sequences,

you can convert it to three files with 100 genes by entering the values 100 200.

For Example: ./camamed_kegg_annotation.py -pgc 100 200

 

  -egk  Extracting information from GhostKOALA or KAAS output files.

In this step, the data of the files generated by GhostKOALA or KAAS server are read and the required information is extracted.

To do this, first, copy the text files into CAMAMED/kegg_annotation and insert the filenames based on the example with space.

Be sure to arrange the file names based on the gene number.

For Example: ./camamed_kegg_annotation.py -egk gk_file1.txt gk_file2.txt

 

  -kla  KEGG local anotation  

We extracted all the annotated information related to the KOs and EC numbers and reactions to the date 2019/6/30.

from the KEGG database and saved in the 'kegg_annotation' folder in the text files that start with the def prefix.

At this step, we extract all the annotated information associated with the sequences step by step.

For Example: ./camamed_kegg_annotation.py -kla

 

  -kua  KEGG user annotation

We extracted all the annotated information related to the KOs and EC numbers and reactions to the date 2019/6/30

from the KEGG database and stored in the CAMAMED/kegg_annotation folder in the text files that start with the def prefix.

But if you want to extract this information yourself, you can use these options.

****These steps may take a long time****

options:

  -koec  Extract KO-related EC numbers from the KEGG database.

  -ecre  Extract EC-related reactions from the KEGG database.

  -reeq  Extract reaction-related equation from the KEGG database.

For Example: ./camamed_kegg_annotation.py -kua -koec

    ./camamed_kegg_annotation.py -kua -ecre

    ./camamed_kegg_annotation.py -kua -reeq

 

> ./camamed_data_normalization.py -h

# ./camamed_data_normalization.py -h

 

This function Extract samples information in normalized bacteria, gene, KO, EC number and reaction matrix.

 

Command: ./camamed_data_normalization.py method [options]

 

      method:

  -h Shows help related to this function

  -nlk  Normalizing data using local information previously extracted from the KEGG database.

option:

  -mtct  Minimum total counts for each taxon. This value is percentage of all taxa, for example 0.5 (default=0.001).

  -mtcg  Minimum total counts of mapped reads per gene in total samples (default=5).

For Example: ./camamed_data_normalization.py -nlk

    ./camamed_data_normalization.py -nlk -mtct 0.002 -mtcg 10

 

  -nuk  Normalizing data using information extracted by the user from the KEGG database.

option:

  -mtct  Minimum total counts for each taxon. This value is percentage of all taxa, for example 0.5 (default=0.001).

  -mtcg  Minimum total counts of mapped reads per gene in total samples (default=5).

For Example: ./camamed_data_normalization.py -nuk

    ./camamed_data_normalization.py -nuk -mtct 0.002 -mtcg 10

>  ./camamed_statistical_test.py -h

# ./camamed_statistical_test.py -h

 

This function perform the statistical test (Kruskal-Wallis H-test or ANOVA test) on normalized data

At this stage, we perform the statistical test on normalized bacteria, gene, KO, EC number, and Reaction data

that stored in the all_results folder.

To get started, you first select the label of class for each sample in the /Read_files/class_label.txt file.

The number of classes can be in the range [2:10]. For each sample, enter a separate row. for example:

class1

class2

class1

class3

 

Command: ./camamed_statistical_test.py method [options]

 

      method:

  -h Shows help related to this function

  -kwd  Running the statistical test on the default annotated data.

option:

  -tsty  Statistical type (Kruskal-Wallis H-test or ANOVA test)[krw/ano]

  -pval  The p-value to filter the output. This value can be in the interval [0:1] (default=0.05)

For Example: ./camamed_statistical_test.py -kwd -tsty krw

    ./camamed_statistical_test.py -kwd -tsty ano

    ./camamed_statistical_test.py -kwd -tsty krw -pval 0.01

    ./camamed_statistical_test.py -kwd -tsty ano -pval 0.01

  -kwu  Running the statistical test on the user extracted data.

option:

  -tsty  Statistical type (Kruskal-Wallis H-test or ANOVA test)[krw/ano]

  -pval  The p-value to filter the output. This value can be in the interval [0:1] (default=0.05)

For Example: ./camamed_statistical_test.py -kwu -tsty krw

    ./camamed_statistical_test.py -kwu -tsty ano

    ./camamed_statistical_test.py -kwu -tsty krw -pval 0.01

    ./camamed_statistical_test.py -kwu -tsty ano -pval 0.01

 

./camamed_mapping_mosaik.py -h

# ./camamed_mapping_mosaik.py -h

 

This function maps the reads to the genes catalog using the MOSAIK software.

To mapping the read sequences to the gene catalog, the following two steps should be implemented.

In the first step, the gene catalog should be converted into an acceptable form for the MOSAIK software.

In the second step, the reads are mapped to the gene catalog. 

 

Command: ./camamed_mapping_mosaik.py method [options]

 

      method:

  -h Shows help related to this function

  -cag  Creating an acceptable form of gene catalog.

options:

  -hw (hash word): Select the length of the hash word that can be in range 4 to 32 (default=15).

For Example: ./camamed_mapping_mosaik.py -cag

    ./camamed_mapping_mosaik.py -cag -hw 17

 

  -mrc  Mapping reads to the gene catalog using MOSAIK software.

options:

  -c (core_number): Select number of cores for MOSAIK execution (default=1).

    The number of cores must be a positive integer.

For Example: ./camamed_mapping_mosaik.py -mrc

    ./camamed_mapping_mosaik.py -mrc -c 5

 

 

 

 実行方法

1、CAMAMEDはマッピングベースの手法のため、最適なマッピング結果を得るためには、適切な遺伝子カタログを提供する必要がある。メタゲノムのリード配列は、MOSAIK v2.2.3ソフトウェアを使用して遺伝子カタログにマッピングされる。

 

2、遺伝子カタログにリードをマッピングした後、各サンプルの遺伝子頻度を計算する。機能プロファイリングには、遺伝子とゲノムの包括的なリストとそれらの生化学的アノテーションを含むKEGGデータベースを使用する。各遺伝子に関連するKEGGオルソロジー(KO)グループを抽出するには、KEGG Automatic Annotation Server(KAAS, Ver.2.1)(ref.42)を利用する。このツールを用いて、ヌクレオチド配列とアミノ酸配列をKEGGマッピングし、関連するKO群を検索する。また、アミノ酸配列をKEGGマッピングするGhostKOALA Ver.2.2(43)を使用してもよい。各遺伝子に関連するKOを抽出した後、各KOについて、酵素コミッション(EC)番号と反応IDを取得する。

 

3、遺伝子の頻度を抽出した後、これらのデータを組成バイアスのために補正する。このステップの前に、すべてのサンプル中のマッピングされたりーどが5未満の遺伝子は排除される。次に、各遺伝子の頻度を遺伝子の長さで割って、長さの異なる遺伝子の正規化されたアバンダンスを計算し、KO、EC数、反応頻度を計算する。そして、組成補正と正規化には、CSS法を用いたmetagenomeSeqパッケージを使用する。この方法を用いて、分類学データの種レベルでの組成の偏りを補正している。KO, EC数, 反応の正規化度数を計算するために、それらのサブセット遺伝子の正規化度数の和を使用する。

 

カレントと共有してイメージのラン

sudo docker run --rm -itv $PWD:/data camamed/camamed_pipeline_db:latest

 

1、準備

ペアエンドfastqファイルは、/camamed/sample_input_files/Read_files/paired_end_fastq/に用意されている。これを Read_files/にコピーする。

cp /camamed/sample_input_files/Read_files/paired_end_fastq/* /camamed/Read_files/
cp sample_input_files/gene_catalog/gene_catalog.fa

 

以下のコマンドを実行する。

./camamed_pre_processing.py -gec -gc gene_catalog.fa -rff fastq -rt p -is 350
  • -gec    Get information on sequences and gene catalogs. options:
      -gc (gene_catalog): Gene catalog sequences must have Fasta format.
            First copy the gene catalog into the /CAMAMED folder.
  • -rff   (read_files_format)[fastq/fasta]
  • -rt (read_type): paired end or single end [p/s]
  • -pgc  Peprocessing gene catalog.
    At this point, the names of the genes are deleted from the gene catalog and stored in the /CAMAMED/files/gene_name.txt
  • -is (insert_size): For paired end sequences (An integer number) or ignor enter -1 (default=-1).

Rad_files/sample_file_names.txtが出力される。

> cat cat Read_files/sample_file_names.txt

p_file1_1.fq

p_file1_2.fq

p_file2_1.fq

p_file2_2.fq

p_file3_1.fq

以後のコマンドではfastqのリストが自動で認識される。すなわち、CAMAMED/Read_files/sample_file_names.txtに記載されたfastqが使用される。

 

 

2、シークエンシングリードの前処理。Fastq/fastaに対して、Fastqcや SeqKit ツールを実行 。3つのコマンドを順番に実行する。

./camamed_quality_control.py -fqc
./camamed_quality_control.py -sek
./camamed_quality_control.py -eti
  •  -fqc   Execut FastQC quality control only for Fastq files. For Example: ./camamed_quality_control.py -fqc
  • -sek   Execut SeqKit to extract information from sample files. For Example: ./camamed_quality_control.py -sek

  • -eti   Extract total information from SeqKit outputs. For Example: ./camamed_quality_control.py -eti

 

 

3、Metaphlan2を用いた細菌の豊富さの抽出。Metaphlan2は様々なレベルでの分類学的プロファイリングを行うことができる。ここでは種レベルを選択。

./camamed_metaphlan_profiling.py -mph -tl s -c 12
./camamed_metaphlan_profiling.py -imp
  •  -mph   Execute metaphlan2. options:
        -tl (taxa_level): For selecting Kingdom, Phylum, Class, Order, Family, Genus, or Species, Select one of {'k', 'p', 'c', 'o', 'f', 'g', 's'}
  • -c   select number of cores for Metaphlan2 execution (default=1).
  • -imp   Extract information from metaphlan2 output files.

 

4、MOSAIKを用いてリード配列を遺伝子カタログにマッピングする。最初のステップでは、遺伝子カタログを MOSAIK ソフトウェアで使用可能な形式に変換する。第二のステップでは、リードを遺伝子カタログにマッピングする。

./camamed_mapping_mosaik.py -cag
./camamed_mapping_mosaik.py -mrc -c 12
  • -cag    Creating an acceptable form of gene catalog. options:
      -hw (hash word): Select the length of the hash word that can be in range 4 to 32 (default=15).

  • -mrc   Mapping reads to the gene catalog using MOSAIK software.

  • -c (core_number):  Select number of cores for MOSAIK execution (default=1).
        The number of cores must be a positive integer.

ここでトラブル発生。以下、未テスト。 

 

 KEGGデータベースからアノテーション情報を抽出。

./camamed_kegg_annotation.py -pgc 100 200
./camamed_kegg_annotation.py -egk gk_file1.txt gk_file2.txt

  

./camamed_kegg_annotation.py -kla
./camamed_data_normalization.py -nlk
./camamed_statistical_test.py -kwd -tsty krw

 

./camamed_kegg_annotation.py -kua -koec
./camamed_kegg_annotation.py -kua -ecre
./camamed_kegg_annotation.py -kua -reeq
./camamed_data_normalization.py -nuk
./camamed_statistical_test.py -kwu -tsty krw

 

引用

CAMAMED: a pipeline for composition-aware mapping-based analysis of metagenomic data
Mohammad H Norouzi-Beirami, Sayed-Amir Marashi, Ali M Banaei-Moghaddam, Kaveh Kavousi
NAR Genomics and Bioinformatics, Volume 3, Issue 1, March 2021

 

参考

分野は異なりますが、組成データについては以下の資料を読みました。


関連