有时间可以写一个megan的中文教程。详细介绍软件、数据库安装，示例数据分析和结果描述。这个也引用了几千次，可以学一下。该软件还支持跨平台，我们分析在Linux和Windows上安装和测试，供不同用户选择使用。

MEGAN-宏基因组功能和物种分类

MEGAN用于功能和物种分类官网。这里有付费版和免费版本。我们以免费版为例。

MEGAN功能简介

MEGAN6是用于交互式分析微生物组数据的综合工具箱。在一个交互式工具。使用NCBI进行物种分类或自定义数据库分类(或SILVA)进行分类分析；

使用InterPro2GO，SEED，eggNOG或KEGG进行功能注释；

简单可视化 条形图，词云，树形图和许多其他图表；

多元统计分析：PCoA，聚类和网络；

支持元数据(metadata)

MEGAN支持许多不同类型的数据输入

原理简要示意图

通过Diamond对序列进行比对(nr)，对输出的daa比对结果文件进行功能和物种使用MEGAN注释。当然不止可以使用Diamond比对结果，还可以使用blast的结果。

MEGAN特有文件格式：RMA

这是MEGAN自己的文件格式，用于存储序列和数据库比对结果，就是RMA格式文件，以.rma格式为后缀，比如BLAST结果，当然这里还有我们的 Diamond输出结果。这两个典型的序列比对输出类似，但是BLAST功能更加强大，Diamond在处理大的数据时速度更快。MEGAN的RMA文件也逐渐升级到RMA6，速度更快，体积更小(仅仅需要原来RMA文件的三分之一的体积)，原来的就RMA文件仍然可以在新版本的MEGAN中打开。可以在软件中通过File ——> Import From BLAST…导入。

RMA文件以许多标题行开头，每行以开头。 这些行可以以任何顺序出现。@Creator MEGAN (version 4.0alpha20, built 14 Oct 2010)

@CreationDate Wed Oct 27 17:10:52 CEST 2010

@ContentType Summary4

@Names 155_PE_1_fixed-paired ecoli-testrun-2000-nr

@Uids 1288068180866 1288190195887

@Sizes 51246 2000

@TotalReads 200000

@Collapse SEED 2000041

@Algorithm Taxonomy tree-from-summary

@Parameters normalizedTo=100000

@NodeStyle KEGG piechart

前两行是软件及其版本信息和作者信息。第三行标识

注明格式为Summary4，表示这是MEGAN 4的总结文件。

第四行列出了此文件代表的所有样本的名称，这里有两个样本。第5行为样本唯一标识符编号(如果有才展示)。第6行列出了原始

样本大小。第7行列出了序列的总数。第8行针对SEED分类指定展示中图形树中的节点分类。这里除了SEED数据库，可以用TAXONOMY或KEGG代替，例如

其他分类。第9行包含用于计算分类的算法的名称。第二个是参数。第10行列出了运行参数用于生成文件。第11行指定分类节点的样式。

MEGAN下载

MEGAN提供了三种版本：Win MEGAN_Community_windows-x64_6_18_4.exe

MAC MEGAN_Community_macos_6_18_4.dmg

Linux MEGAN_Community_unix_6_18_4.sh

如Linux版下载# 安装程序 102M

wget -c https://software-ab.informatik.uni-tuebingen.de/download/megan6/MEGAN_Community_unix_6_18_4.sh

# NCBI-nr编号与物种和功能注释 970M

wget -c https://software-ab.informatik.uni-tuebingen.de/download/megan6/megan-map-Oct2019.db.zip

# 核酸与物种信息 655MB

wget -c https://software-ab.informatik.uni-tuebingen.de/download/megan6/megan-nucl-Oct2019.db.zip

提供数据库将NCBI-nr数据库比对文件注释到分类和功能：(taxonomy，eggNOG，InterPro2GO和SEED)，但是免费版本就只能使用这只是到这四个，并在使用前解压缩：megan-map-Oct2019.db.zip

当然还有需要许可证的收费版本：数据库就包含了KEGG。点击此处申请密匙 https://computomics.com/megan.html

数据库也不同于社区版本(megan-map-Oct2019-ue.db.zip)。我尝试申请了使用密匙，但是三天了也还没消息。

MEGAN使用

MEGAN(linux版本安装)

直接在terminal中运行，会弹出图形界面，按照提示安装即可，如果不选择位置，则在home下生成一个megan的文件夹。

软件安装# 方法1. conda安装 http://bioconda.github.io/ 6.12.3-0 built 14 Aug 2018，构建一个环境，但不是最新版

conda create -n megan # 创建megan环境

conda activate megan # 进入megan环境

conda install megan # 安装megan，235Mb 版本太低

# 方法2. 直接安装

wget -c https://software-ab.informatik.uni-tuebingen.de/download/megan6/MEGAN_Community_unix_6_18_4.sh

bash MEGAN_Community_unix_6_18_4.sh

# JVM must be at least 11. Please define INSTALL4J_JAVA_HOME to point to a suitable JVM

java -version # 1.8.0_201

# 安装完上面conda会变为 openjdk version "11.0.1-internal" 2018-10-16

数据库下载# NCBI-nr编号与物种和功能注释 970M

wget -c https://software-ab.informatik.uni-tuebingen.de/download/megan6/megan-map-Oct2019.db.zip

# 解压后为 5 Gb 的db文件

unzip megan-map-Oct2019.db.zip

# 核酸与物种信息 655MB

wget -c https://software-ab.i   nformatik.uni-tuebingen.de/download/megan6/megan-nucl-Oct2019.db.zip

# 解压后为 4.2 Gb 的db文件

unzip megan-nucl-Oct2019.db.zip

# NCBI-nr完整蛋白序列库和建diamond索引，2020/2/4，67G

wget -c ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz

#尝试使用pigz多线程解压缩，123G

time unpigz -k -p 16 nr.gz # 8m, 26m

#gunzip -c nr.gz > nr

time diamond makedb --in nr -d nr -p 32 # 8m, 102m

MEGAN使用指南(Linux)# 下载

wget -c ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR793/ERR793599/ERR793599_1.fastq.gz

wget -c ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR793/ERR793599/ERR793599_2.fastq.gz

# 以任意的宏基因组数据为例，ERR793599

#重压缩测序文件

pigz -p 6 -dc ./ERR793599_2.fastq.gz | pigz -p 6 > ./C1_2.fastq.gz

pigz -p 6 -dc ./ERR793599_1.fastq.gz | pigz -p 6 > ./C1_1.fastq.gz

#去除barcode并进行指控

java -jar ~/sra/Trimmomatic-0.38/trimmomatic-0.38.jar PE -threads 6 \

-phred33 ./unpack/C1_1.fastq.gz ./unpack/C1_2.fastq.gz \

./trimmomatic/C1_forward_paired.fq.gz ./trimmomatic/C1_forward_unpaired.fq.gz ./trimmomatic/C1_reverse_paired.fq.gz ./trimmomatic/C1_reverse_unpaired.fq.gz \

ILLUMINACLIP:../../database/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:30 MINLEN:100 2> C1.log

#使用nr数据框对前端数据进行比对，每个样本可能需要至少3-5小时，由数据大小决定

diamond blastx -c 1 --db /home/wentao/Desktop/biostack/database/nr/db/nr.dmnd  -t /tmp -p 34 -q ./trimmomatic/C1_forward_paired.fq.gz --daa ./diamond/C1.1.daa

#使用nr数据框对后端数据进行比对

diamond blastx -c 1 --db /home/wentao/Desktop/biostack/database/nr/db/nr.dmnd  -t /tmp -p 34 -q ./trimmomatic/C1_reverse_paired.fq.gz --daa ./diamond/C1.2.daa

#转化双端daa文件为MEGAN特有额rma文件。

~/megan/tools/daa2rma -i ./diamond/C1.1.daa ./diamond/C1.2.daa --paired -ms 50 -me 0.01 -top 50  -mdb ~/db/megan/megan-map-Oct2019.db  -o ./diamond/C1.rma

运行过程文件展示Parsing file: ./diamond/C1.1.daa

10% 20% 30% 40% 50% 60% 70% 80% 90% Parsing file: ./diamond/C1.2.daa

10% 20% 30% 40% 50% 60% 70% 80% 90% 100% (810.8s)

Total reads:           443,738

Alignments:          9,103,355

100% (0.0s)

100% (0.0s)

100% (0.0s)

Linking paired reads

Number of pairs:       178,830

Binning reads: Initializing...

Initializing binning...

Using paired reads in taxonomic assignment...

Using 'Naive LCA' algorithm for binning: Taxonomy

Using Best-Hit algorithm for binning: SEED

Using Best-Hit algorithm for binning: EGGNOG

Using Best-Hit algorithm for binning: INTERPRO2GO

Binning reads...

Binning reads: Analyzing alignments

Total reads:          443,738

With hits:             443,738

Alignments:          9,103,355

Assig. Taxonomy:       443,738

Assig. SEED:            53,014

Assig. EGGNOG:          85,333

Assig. INTERPRO2GO:    204,180

MinSupport set to: 221

Binning reads: Applying min-support & disabled filter to Taxonomy...

Min-supp. changes:       1,599

Binning reads: Writing classification tables

Numb. Tax. classes:         92

Numb. SEED classes:      3,055

Numb. EGG. classes:      4,756

Numb. INT. classes:      6,832

Binning reads: Syncing

Class. Taxonomy:            92

Class. SEED:             3,055

Class. EGGNOG:           4,756

Class. INTERPRO2GO:      6,832

100% (728.6s)

Total time:  1547s

Peak memory: 24.6 of 125.8G

提取注释内容(物种和功能)：rma2info： 用于提取rma文件中的物种和功能注释信息。

提取物种注释数据：~/megan/tools/rma2info -i ./diamond/C1.rma -r2c Taxonomy -n true -v >  ./diamond/C1Taxonomy.txtn true 显示菌名称

paths true  显示层级菌名

—ranks: true 显示注释到那个等级，会在序列前面加上菌等级的标记，P，S等等

—list true 添加简单序列统计信息

—listMore true 添加运行参数等信息

我们把这些参数添加上再运行一次：~/megan/tools/rma2info -i ./diamond/C1.rma -c2c Taxonomy -r2c Taxonomy -n true --paths true --ranks true --list true --listMore true --bacteriaOnly true -v >  ./diamond/C1Taxonomy1.txt

提取功能：

SEED：功能注释数据库。07年有篇文章介绍了MEGAN分析SEED和KEGG：https://bmcbioinformatics.biomedcentral.com/articles/10.1186/1471-2105-12-S1-S21

EGGNOG：功能注释；

INTERPRO2GO：将序列分类为蛋白质家族，并预测重要结构域和位点的存在。InterPro是蛋白质家族，域和位点的集成资源。参考网址学习

提取EGGNOG注释：~/megan/tools/rma2info -i ./diamond/C1.rma -r2c EGGNOG -n true --paths true --ranks true --list true --listMore true -v >  ./diamond/C1eggnog.txt

提取SEED注释：~/megan/tools/rma2info -i ./diamond/C1.rma -r2c SEED -n true --paths true --ranks true --list true --listMore true -v >  ./diamond/C1SEED.txt

提取INTERPRO2GO注释：~/megan/tools/rma2info -i ./diamond/C1.rma -r2c INTERPRO2GO -n true --paths true --ranks true --list true --listMore true -v >  ./diamond/C1INTERPRO2GO.txt

Win版安装和使用

MEGAN安装

软件安装：MEGAN提供win下的32和64为版本可供下载，这里我们选择64位;32位Windows版本的MEGAN配置只能使用1.1 GB内存。对于所有其他版本在安装过程中，安装程序将允许您在安装过程中设置最大内存。在默认情况下，程序将建议使用2 GB。如果你的电脑有更多的可用内存，那么，把这个限制定得更高。例如，如果您有4 GB的主内存，则设置将MEGAN限制为3 GB。8GB内存可以设置5GB、16GB可设置12G。这是因为所给的内存就越多，程序运行得越快。MEGAN依赖java，这里我我们首先要下载安装java8以上的版本；

官网选择exe程序下载；

双击安装程序-同意协议，开始安装

设置最大内存使用量：我的电脑8G内存，可设置为5G，大家可以根据自己实际条件设置分配内存数量。

安装完成后，打开软件：自动加载NCBI的物种三级分类：

Win版使用指南

MEGAN主界面介绍

主界面用于展示门类树并且可以通过菜单栏去控制门类树。在初始情况下，也就是刚刚打开MEGAN，还没有导入RMA文件的时候，主界面展示NCBI的分类树，默认展示三层结构，全部展示可以使用菜单栏的工具：rank。

一旦数据读入，NCBI的分类树就会被我们导入文件的树所替代，树中枝节点大小代表了注释到该节点的物种所匹配的序列数量。双击节点我们可以看到一些信息，例如：匹配到该节点的序列数量，总序列数量。

下游节点可以被折叠和展开，这通过菜单栏可以解决。鼠标右键也可以访问菜单栏。

MEGAN输入介绍The File→New… item:打开空白项目；

The File→Open… item: 打开MEGAN文件 (e文件后缀名： .rma, .meg 或者 .megan).

The File→Open Recent submenu ：打开最近的项目；

The File→Open From Server… item: 从MeganServer服务器上打开文件，(这也是MEGAN的示例文件)；

The File→Compare… item: 打开比较对话框，对多个数据集进行比较；

The File→Import From BLAST… item: 打开序列比对或者序列文件(daa)；

文件输入：这里我们可以输入blast或者diamond序列比对后输出的daa文件。首先是示例文件，这里我们选择MEGAN服务器上的文件；，按照如下操作，我们可以看到有好几组示例数据，这里我们直接选择第一个做演示：(打开速度与网速有关，因为文件在网络服务器上)

选择之后就出现了如下界面：这是物种分类树。

可以通过这几个按钮查看物种和功能信息树(这里导入rma6文件)，现在MEGAN改版了，只要注释出来RMA文件，都会有物种和功能，不用选择数据库了：

无论是物种注释，还是功能都是以树的形式展示，也就是说一定有一个根，对于树，我们可以操作节点，选择收起子节点还是展示子节点：

当我们了解了该样本物种的功能组成后，需要导出对应数据。我们以SEED蛋白注释为例：摁住shift键，然后使用鼠标左键拉去需要保存的树，变为黄色的也就是目前可以保存的数据(其他功能和物种注释保存类似操作)：

摘取部分数据查看：第一列是功能单位，第二列是丰度信息(可以选择是否输出标准化信息)"SEED"     7.0343728E7

"Amino Acids and Derivatives"     1303462.0

"Arabinose Sensor and transport module"     157.0

"Autotrophy"     2491.0

"2-Ketogluconate Utilization"     5.0

"2-O-alpha-mannosyl-D-glycerate utilization"     2141.0

"Acetoin, butanediol metabolism"     25694.0

"Acetone Butanol Ethanol Synthesis"     36732.0

"Acetyl-CoA biosynthesis in plants"     21184.0

"Acetyl-CoA fermentation to Butyrate"     29880.0

"acinetobacter tca"     97660.0

"alpha carboxysome"     25116.0

"Alpha-acetolactate operon"     4.0

"Alpha-Amylase locus in Streptocococcus"     38.0

"beta carboxysome"     3994.0

"Beta-Glucoside Metabolism"     12943.0

"Butanol Biosynthesis"     28782.0

"Calvin-Benson cycle"     79542.0

"Calvin-Benson-Bassham cycle in plants"     50638.0

"Carboxysome"     21503.0

"Chitin and N-acetylglucosamine utilization"     149992.0

"CitAB"     1.0

KEGG注释查看："KEGG2011"     7.034416E7

"Carbohydrate Metabolism"     3721914.0

"K1000010 Glycolysis / Gluconeogenesis"     422395.0

"K1000020 Citrate cycle (TCA cycle)"     310392.0

"K1000030 Pentose phosphate pathway"     329485.0

"K1000040 Pentose and glucuronate interconversions"     340106.0

"K1000051 Fructose and mannose metabolism"     425046.0

"K1000052 Galactose metabolism"     816922.0

"K1000053 Ascorbate and aldarate metabolism"     53282.0

"K1000500 Starch and sucrose metabolism"     745091.0

"K1000520 Amino sugar and nucleotide sugar metabolism"     877177.0

"K1000620 Pyruvate metabolism"     492159.0

"K1000630 Glyoxylate and dicarboxylate metabolism"     212640.0

"K1000640 Propanoate metabolism"     242167.0

"K1000650 Butanoate metabolism"     198659.0

"K1000660 C5-Branched dibasic acid metabolism"     77196.0

"K1000562 Inositol phosphate metabolism"     33097.0

"Energy Metabolism"     1609134.0

"K1000190 Oxidative phosphorylation"     406172.0

"K1000195 Photosynthesis"     113692.0

"K1000710 Carbon fixation in photosynthetic organisms"     307788.0

"K1000720 Carbon fixation pathways in prokaryotes"     216803.0

"K1000680 Methane metabolism"     379958.0

"K1000910 Nitrogen metabolism"     382561.0

"K1000920 Sulfur metabolism"     95305.0

"Lipid Metabolism"     968704.0

"K1000061 Fatty acid biosynthesis"     164601.0

EGGNOG注释查看："eggNOG"     7.0344072E7

"informationStorageAndProcessing"     4313603.0

"cellularProcessesAndSignaling"     5433823.0

"metabolism"     1.0153629E7

"[C] Energy production and conversion"     1243786.0

"[E] Amino acid transport and metabolism"     2183752.0

"[F] Nucleotide transport and metabolism"     738606.0

"[G] Carbohydrate transport and metabolism"     3057650.0

"[H] Coenzyme transport and metabolism"     712124.0

"[I] Lipid transport and metabolism"     607947.0

"[P] Inorganic ion transport and metabolism"     1538295.0

"[Q] Secondary metabolites biosynthesis, transport and catabolism"     155732.0

"Not assigned"     5.0519228E7

InterPro2GO注释结果查看："InterPro2GO"     6.8526864E7

"GO:0071973 bacterial-type flagellum-dependent cell motility"     48865.0

"GO:0007155 cell adhesion"     23.0

"GO:0045454 cell redox homeostasis"     45236.0

"GO:0071840 cellular component organization or biogenesis"     740858.0

"GO:0006259 DNA metabolic process"     1580534.0

"GO:0016226 iron-sulfur cluster assembly"     47578.0

"GO:0008152 metabolic process"     1.0141539E7

"GO:0008218 bioluminescence"     29.0

"GO:0009058 biosynthetic process"     3593054.0

"GO:0005975 carbohydrate metabolic process"     3166996.0

"GO:0006091 generation of precursor metabolites and energy"     416010.0

"GO:0006629 lipid metabolic process"     503233.0

"GO:0015948 methanogenesis"     190.0

"GO:0006807 nitrogen compound metabolic process"     3833028.0

"GO:0016310 phosphorylation"     277450.0

"GO:0015979 photosynthesis"     276.0

"GO:0044281 small molecule metabolic process"     4018961.0

"GO:0006351 transcription, DNA-templated"     185610.0

"GO:0006412 translation"     328870.0

"IPR003821 1-deoxy-D-xylulose 5-phosphate reductoisomerase"     21365.0

"IPR006401 2,5-diamino-6-hydroxy-4-(5-phosphoribosylamino)pyrimidine 1-reductase, archaeal"     1.0

TAX物种注释信息查看：默认输出只有灭有level层次信息，可选输出stamp，可以输出带有层级信息的注释文件。Level_1    Level_2    Level_3    Level_4    Level_5    Level_6    Level_7    Observation Ids    Bob06

k__(null)    p__(null)    c__(null)    o__(null)    f__(null)    g__(null)    s__(null)    ID1    462699.0

k__(null)    p__(null)    c__(null)    o__(null)    f__(null)    g__(null)    s__(null)    ID131567    80876.0

k__Bacteria    p__(Bacteria)    c__(Bacteria)    o__(Bacteria)    f__(Bacteria)    g__(Bacteria)    s__(Bacteria)    ID2    2493265.0

k__Bacteria    p__(Bacteria)    c__(Bacteria)    o__(Bacteria)    f__(Bacteria)    g__(Bacteria)    s__(Bacteria)    ID1783270    0.0

k__Bacteria    p__(Bacteria)    c__(Bacteria)    o__(Bacteria)    f__(Bacteria)    g__(Bacteria)    s__(Bacteria)    ID68336    6045.0

k__Bacteria    p__Bacteroidetes    c__(Bacteroidetes)    o__(Bacteroidetes)    f__(Bacteroidetes)    g__(Bacteroidetes)    s__(Bacteroidetes)    ID976    721448.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__(Bacteroidia)    f__(Bacteroidia)    g__(Bacteroidia)    s__(Bacteroidia)    ID200643    0.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__(Bacteroidales)    g__(Bacteroidales)    s__(Bacteroidales)    ID171549    7600487.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__(Bacteroidaceae)    s__(Bacteroidaceae)    ID815    0.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__Bacteroides    s__(Bacteroides)    ID816    2.9878604E7

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__Bacteroides    s__Bacteroides caccae    ID47678    301910.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__Bacteroides    s__Bacteroides caccae    ID411901    100571.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__Bacteroides    s__Bacteroides dorei    ID357276    914270.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__Bacteroides    s__Bacteroides dorei    ID556260    101864.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__Bacteroides    s__Bacteroides dorei    ID997877    85211.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__Bacteroides    s__Bacteroides dorei    ID483217    225193.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__Bacteroides    s__Bacteroides fragilis    ID817    142742.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__Bacteroides    s__Bacteroides ovatus    ID28116    97339.0

k__Bacteria    p__Bacteroidetes    c__Bacteroidia    o__Bacteroidales    f__Bacteroidaceae    g__Bacteroides    s__Bacteroides uniformis    ID820    112694.0

如此我们得到了样本物种注释和功能注释的文件，可用于后续分析了。

MEGAN分析进阶

双样本比对

为了进行双样本比对，我们选择了另外一个样本按照同样的流程再来一遍；

对ERR793603号样本采取C1的流程，再来一遍：体验合并rma文件的功能和MEGAN图形界面可视化的功能。

这种模式也是我们通常基于测序数据分析的流程。重点

数据质控#重压缩测序文件

pigz -p 6 -dc ./ERR793603_1.fastq.gz | pigz -p 12 > ./S1_1.fastq.gz

pigz -p 6 -dc ./ERR793603_2.fastq.gz | pigz -p 12 > ./S1_2.fastq.gz

#去除测序接头和引物序列并进行质控(软件、输入文件和接头序列位置)

java -jar ~/sra/Trimmomatic-0.38/trimmomatic-0.38.jar  PE -threads 6 -phred33 ./unpack/C1_1.fastq.gz ./unpack/C1_2.fastq.gz  ./trimmomatic/C1_forward_paired.fq.gz ./trimmomatic/C1_forward_unpaired.fq.gz ./trimmomatic/C1_reverse_paired.fq.gz ./trimmomatic/C1_reverse_unpaired.fq.gz ILLUMINACLIP:../../database/NexteraPE-PE.fa:2:30:10 LEADING:3 TRAILING:3 SLIDINGWINDOW:4:30 MINLEN:100 2> C1.log

比对生成megan输入文件#使用nr数据框对前端数据进行比对

diamond blastx -c 1 --db /home/wentao/Desktop/biostack/database/nr/db/nr.dmnd  -t /tmp -p 34 -q ./trimmomatic/S1_forward_paired.fq.gz --daa ./diamond/S1.1.daa

#使用nr数据框对后端数据进行比对

diamond blastx -c 1 --db /home/wentao/Desktop/biostack/database/nr/db/nr.dmnd  -t /tmp -p 34 -q ./trimmomatic/S1_reverse_paired.fq.gz --daa ./diamond/S1.2.daa

#转化双端daa文件为MEGAN特有额rma文件。

~/megan/tools/daa2rma -i ./diamond/S1.1.daa ./diamond/S1.2.daa --paired -ms 50 -me 0.01 -top 50  -mdb ~/db/megan/megan-map-Oct2019.db  -o ./diamond/MEGAN_RMA/S1.rma

合并不同样本的.rma文件#运行命令，需要图型界面支持，如远程桌面，或本地安装配置Xing/Xmanager

~/megan/tools/compute-comparison -i ./diamond/MEGAN_RMA/S1.rma ./diamond/MEGAN_RMA/C1.rma  -o ./compared_all -v true

这个输出文件：这就是MEGAN的标准输出文件：以TAX为例：第一列是NCBI id

，第二列是count数量。@CreationDate    Tue Feb 04 17:49:13 CST 2020

@ContentType    Summary4

@Names    C1    S1    C11

@BlastMode    BlastX    BlastX    BlastX

@Uids    1580592065210    1580804730040    1580592065210

@Sizes    443738.0    751161.0    443738.0

@TotalReads    1638637

@Collapse    Taxonomy    28384    2    12884    2759    12908    2157    10239

@Parameters    minScore=0.0 maxExpected='10000.0' minPercentIdentity='0.0' topPercent=100.0 minSupportPercent=0.0 minSupport=1 lcaAlgorithm=naive minPercentReadToCover=0.0 minPercentReferenceToCover=0.0 minComplexity=0.0 longReads=false pairedReads=false identityFilter=false readAssignmentMode=readCount fNames= { Taxonomy SEED EGGNOG INTERPRO2GO }

@NodeStyle    Taxonomy    BarChart

@NodeStyle    SEED    PieChart

@NodeStyle    EGGNOG    PieChart

@NodeStyle    INTERPRO2GO    PieChart

@ColorTable    Fews8    White-Green

@ColorEdits

TAX    1    1761.0    2274.0

TAX    2    40023.0    80016.0

TAX    1437201    46377.0    78736.0

TAX    5125    5877.0    8092.0

·······

INTERPRO2GO    16380    17.0    16.0

INTERPRO2GO    16381    1.0    0.0

INTERPRO2GO    16382    0.0    2.0

END_OF_DATA_TABLE

#SampleID    @Source

C1    ./C1.rma

S1    ./S1.rma

导出详细数据参考C1即可。

通过上面的教程我们基本了解了MEGAN的可视化界面，下面我们就两个样本进行可视化操作和分析。

MEGAN界面可视化-两样本(.rma)比对

第一步，导入MEGAN的rma比对注释文件：打开megan，按照图示选择，我们进行两个rma文件的比对

打开之后之后出现这个默认图形：默认节点都是按照颜色填充的，仅仅显示两层物种结构：

物种信息可视化

第二步：首先我们进行简单缩放(鼠标滑轮)，修改配色(点击如图所示工具栏类似马赛克的图表，调整配色)，调节Rank按钮显示的门类水平，可以从界到种：

第三步：物种注释可视化方案：MEGAN提供的物种可视化的方案超出我们的想象，十余种可视化方案可供我们选择，甚至今天我还有几种没有在R中进行尝试，这是一个很好的契机，大家可以进行尝试：

点击图中标记按钮，选择可视化方案，这里我们选择第一个即可进入可视化窗口：

这里我向大家介绍几种图形：；堆叠华夫饼图，多个华夫饼图替换了气泡图的气泡，使用数量来代表丰富，但是必须在尺度范围内进行限制，目前R语言没有成熟的工具实现这个图形：

物种分布词云，这种方式让我们迅速就可以找到丰度较高的物种，或者差异物种。

进入可视化窗口，工具栏有颜色的都是可视化方案，可供大家选择，都可以点击看看(滚动鼠标，拉伸图形)：

可以旋转坐标轴根据样本或者物种注释进行展示图形：点击图中标记为红色 的按钮，进行转换坐标轴

有些图形是没有转换选择的，例如：网络图。

物种可视化的物种来源于我们的选择：意味着我们点击目标菌群，就可以生成只包含目标菌群的物种丰度信息，例如我们只选择bacteria以及去所有的子节点做展示：

功能信息可视化

通过免费版本的数据库注释我们可以得到三个功能，从这里可以看到是哪三个，缺少KEGG和VFDB，这里的三个功能的的操作类似于物种可视化，我简单做一个演示，注意的细节和上面都一样：

只要注意，选择了多少节点，那就可视化多少数据。

数据导出

无论是物种注释，还是功能注释，最终的结果我们当然需要进行一个很好的保存。MEGAN的保存类型很多，这里咱们来学习：点击：Save As保存为.megan文件，这个文件保存后我们就可以随时打开查看数据，再也不需要使用rma文件。

点击Export Image将图片保存

点击Export Legend保存图例

红色选框种可以保存物种或者功能的表格文文件

Metadate保存文件名

Tree保存传统的数文件，这个文件我们可以使用其他树可视化软件出图，比如：Graphlan，ggtree，iTOL等。

MEGAN多个样本分析方案

本次我使用刘老师的宏基因组示例数据，这批数据在每个在100M以内，但是有12个样本，通过对这12个样本的操作，我们重点来做批量MEGAN分析和文件导出。

seq中是全部的原始序列；我们首先进入seq路径下面：# 切换工作路径；

cd ~/Desktop/Shared_Folder/metagenome_study/ori_pipline/seq/

ls ./*.fq.gz

进行序列质控mkdir ../trimmomatic

# 调用for循环批处理文件

for filename in *_1.fq.gz

do

# 提取双端公共文件名，并输出检验

base=$(basename $filename _1.fq.gz)

echo $base

# 运行去接头程序

java -jar ~/sra/Trimmomatic-0.38/trimmomatic-0.38.jar  PE  -threads 32 \

${base}_1.fq.gz \

${base}_2.fq.gz \

../trimmomatic/${base}_forward_paired.fq.gz ../trimmomatic/${base}_forward_unpaired.fq.gz \

../trimmomatic/${base}_reverse_paired.fq.gz ../trimmomatic/${base}_reverse_unpaired.fq.gz \

ILLUMINACLIP:TruSeq2-PE.fa:2:40:15 \

LEADING:2 TRAILING:2 \

SLIDINGWINDOW:4:2 \

MINLEN:25 2> C1.log

done

使用nr数据库进行比对，使用34个线程跑，每个测序样本压缩文件只有10M，但是还是需要耗费将近50min一个。一共12个样本，所以这个步骤耗时将近半天。#建立比对文件夹，不见文件夹会--daa选项会提示错误

cd ../

mkdir ./diamond

# 调用for循环批处理文件

for filename in ./seq/*_1.fq.gz

do

# 提取双端公共文件名，并输出检验

base=$(basename $filename _1.fq.gz)

echo $base

# 前端序列开始比对

diamond blastx -c 1 --db /home/wentao/Desktop/biostack/database/nr/db/nr.dmnd  -t /tmp -p 34 -q ./trimmomatic/${base}_forward_paired.fq.gz --daa ./diamond/${base}.1.daa

# 后端序列

diamond blastx -c 1 --db /home/wentao/Desktop/biostack/database/nr/db/nr.dmnd  -t /tmp -p 34 -q ./trimmomatic/${base}_reverse_paired.fq.gz --daa ./diamond/${base}.2.daa

done

使用MEGAN提取结果mkdir ./diamond/MEGAN_RMA

# 调用for循环批处理文件

for filename in ./seq/*_1.fq.gz

do

# 提取双端公共文件名，并输出检验

base=$(basename $filename _1.fq.gz)

echo $base

#运行命令

~/megan/tools/daa2rma -i ./diamond/${base}.1.daa ./diamond/${base}.2.daa --paired -ms 50 -me 0.01 -top 50  -mdb ~/db/megan/megan-map-Oct2019.db -o ./diamond/MEGAN_RMA/${base}.rma

done

合并全部的物种和功能数据：#运行命令

~/megan/tools/compute-comparison -i ./diamond/MEGAN_RMA/*.rma  -o ./diamond/MEGAN_RMA/compared_all -v true

输出文件：@Creator    ComputeComparison

@CreationDate    Tue Feb 04 19:56:46 CST 2020

@ContentType    Summary4

@Names    p136C    p136N    p143C    p143N    p144C    p144N    p146C    p146N    p153C    p153N    p156C    p156N

@BlastMode    BlastX    BlastX    BlastX    BlastX    BlastX    BlastX    BlastX    BlastX    BlastX    BlastX    BlastX    BlastX

@Uids    1580805529488    1580806456437    1580806849152    1580807063364    1580807282560    1580807470927    1580807820232    1580808148014    1580808496509    1580808796150    1580809145942    1580809441407

@Sizes    111945.0    92730.0    63883.0    64549.0    57390.0    110147.0    102055.0    113277.0    107679.0    114086.0    114115.0    112538.0

@TotalReads    1164394

@Collapse    Taxonomy    28384    2    12884    2759    12908    2157    10239

@Parameters    minScore=0.0 maxExpected='10000.0' minPercentIdentity='0.0' topPercent=100.0 minSupportPercent=0.0 minSupport=1 lcaAlgorithm=naive minPercentReadToCover=0.0 minPercentReferenceToCover=0.0 minComplexity=0.0 longReads=false pairedReads=false identityFilter=false readAssignmentMode=readCount fNames= { Taxonomy SEED EGGNOG INTERPRO2GO }

@NodeStyle    Taxonomy    BarChart

@NodeStyle    SEED    PieChart

@NodeStyle    EGGNOG    PieChart

@NodeStyle    INTERPRO2GO    PieChart

@ColorTable    Fews8    White-Green

@ColorEdits

TAX    -2    1.0

TAX    1    342.0    736.0    685.0    278.0    260.0    449.0    477.0    568.0    161.0    409.0    381.0    574.0

TAX    2049    0.0    134.0    188.0    580.0    41.0    188.0    0.0    0.0    0.0    228.0    14.0    588.0

TAX    2    45382.0    33905.0    13366.0    16258.0    17194.0    40821.0    35587.0    38908.0    41649.0    31388.0    28311.0    37730.0

TAX    143361    743.0

TAX    1437201    0.0    0.0    8.0    4.0    5.0

INTERPRO2GO    16381    8.0    1.0    0.0    0.0    0.0    0.0    0.0    0.0    1.0    0.0    28.0    5.0

END_OF_DATA_TABLE

#SampleID    @Source

p136C    ./diamond/MEGAN_RMA/p136C.rma

p136N    ./diamond/MEGAN_RMA/p136N.rma

p143C    ./diamond/MEGAN_RMA/p143C.rma

p143N    ./diamond/MEGAN_RMA/p143N.rma

p144C    ./diamond/MEGAN_RMA/p144C.rma

p144N    ./diamond/MEGAN_RMA/p144N.rma

p146C    ./diamond/MEGAN_RMA/p146C.rma

p146N    ./diamond/MEGAN_RMA/p146N.rma

p153C    ./diamond/MEGAN_RMA/p153C.rma

p153N    ./diamond/MEGAN_RMA/p153N.rma

p156C    ./diamond/MEGAN_RMA/p156C.rma

p156N    ./diamond/MEGAN_RMA/p156N.rma

多样本批量物种和功能导出for filename in ./diamond/MEGAN_RMA//*.rma

do

# 提取双端公共文件名，并输出检验

base=$(basename $filename .rma)

echo $base

#运行命令

~/megan/tools/daa2rma -i ./diamond/${base}.1.daa

# 提取物种注释信息

~/megan/tools/rma2info -i ./diamond/MEGAN_RMA/${base}.rma -c2c Taxonomy -r2c Taxonomy -n true --paths true --ranks true --list true --listMore true --bacteriaOnly true -v >  ./diamond/MEGAN_RMA/${base}Taxonomy.txt

#提取EGGNOG注释：

~/megan/tools/rma2info -i ./diamond/MEGAN_RMA/${base}.rma -r2c EGGNOG -n true --paths true --ranks true --list true --listMore true -v >  ./diamond/MEGAN_RMA/${base}eggnog.txt

#提取SEED注释：

~/megan/tools/rma2info -i ./diamond/MEGAN_RMA/${base}.rma -r2c SEED -n true --paths true --ranks true --list true --listMore true -v > ./diamond/MEGAN_RMA/${base}SEED.txt

#提取INTERPRO2GO注释：

~/megan/tools/rma2info -i ./diamond/MEGAN_RMA/${base}.rma -r2c INTERPRO2GO -n true --paths true --ranks true --list true --listMore true -v >  ./diamond/MEGAN_RMA/${base}INTERPRO2GO.txt

done

可视化和后续分析参考MEGAN可视化终端，上面我们已经讲的很明白了。大家请参照。撰文：文涛 南京农业大学

责编：刘永鑫 中科院遗传发育所

猜你喜欢

写在后面

为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。

学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”

​