宏基因组注释和可视化神器MEGAN入门

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

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

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

官网 https://www.wsi.uni-tuebingen.de/lehrstuehle/algorithms-in-bioinformatics/software/megan6/

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 66 87 @Sizes 51246 2000 @TotalReads  @Collapse SEED  @Algorithm Taxonomy tree-from-summary @Parameters normalizedTo= @NodeStyle KEGG piechart 

MEGAN下载

MEGAN官网下载页面：https://software-ab.informatik.uni-tuebingen.de/download/megan6/welcome.html

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

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）

原始序列处理：如 https://www.ebi.ac.uk/ena/browser/view/ERR

# 下载 wget -c ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR793/ERR/ERR_1.fastq.gz wget -c ftp://ftp.sra.ebi.ac.uk/vol1/fastq/ERR793/ERR/ERR_2.fastq.gz # 以任意的宏基因组数据为例，ERR #重压缩测序文件 pigz -p 6 -dc ./ERR_2.fastq.gz | pigz -p 6 > ./C1_2.fastq.gz pigz -p 6 -dc ./ERR_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.txt 

• n 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 

提取功能：

提取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.0E7 "Amino Acids and Derivatives" .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" .0 "CitAB" 1.0 

KEGG注释查看：

"KEGG2011" 7.034416E7 "Carbohydrate Metabolism" .0 "K Glycolysis / Gluconeogenesis" .0 "K Citrate cycle (TCA cycle)" .0 "K Pentose phosphate pathway" .0 "K Pentose and glucuronate interconversions" .0 "K Fructose and mannose metabolism" .0 "K Galactose metabolism" .0 "K Ascorbate and aldarate metabolism" 53282.0 "K Starch and sucrose metabolism" .0 "K Amino sugar and nucleotide sugar metabolism" .0 "K Pyruvate metabolism" .0 "K Glyoxylate and dicarboxylate metabolism" .0 "K Propanoate metabolism" .0 "K Butanoate metabolism" .0 "K C5-Branched dibasic acid metabolism" 77196.0 "K Inositol phosphate metabolism" 33097.0 "Energy Metabolism" .0 "K Oxidative phosphorylation" .0 "K Photosynthesis" .0 "K Carbon fixation in photosynthetic organisms" .0 "K Carbon fixation pathways in prokaryotes" .0 "K Methane metabolism" .0 "K Nitrogen metabolism" .0 "K Sulfur metabolism" 95305.0 "Lipid Metabolism" .0 "K Fatty acid biosynthesis" .0 

EGGNOG注释查看：

"eggNOG" 7.0E7 "informationStorageAndProcessing" .0 "cellularProcessesAndSignaling" .0 "metabolism" 1.0E7 "[C] Energy production and conversion" .0 "[E] Amino acid transport and metabolism" .0 "[F] Nucleotide transport and metabolism" .0 "[G] Carbohydrate transport and metabolism" .0 "[H] Coenzyme transport and metabolism" .0 "[I] Lipid transport and metabolism" .0 "[P] Inorganic ion transport and metabolism" .0 "[Q] Secondary metabolites biosynthesis, transport and catabolism" .0 "Not assigned" 5.0E7 

InterPro2GO注释结果查看：

"InterPro2GO" 6.E7 "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" .0 "GO:0006259 DNA metabolic process" .0 "GO:0016226 iron-sulfur cluster assembly" 47578.0 "GO:0008152 metabolic process" 1.0E7 "GO:0008218 bioluminescence" 29.0 "GO:0009058 biosynthetic process" .0 "GO:0005975 carbohydrate metabolic process" .0 "GO:0006091 generation of precursor metabolites and energy" .0 "GO:0006629 lipid metabolic process" .0 "GO:0015948 methanogenesis" 190.0 "GO:0006807 nitrogen compound metabolic process" .0 "GO:0016310 phosphorylation" .0 "GO:0015979 photosynthesis" 276.0 "GO:0044281 small molecule metabolic process" .0 "GO:0006351 transcription, DNA-templated" .0 "GO:0006412 translation" .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 .0 k__(null) p__(null) c__(null) o__(null) f__(null) g__(null) s__(null) ID 80876.0 k__Bacteria p__(Bacteria) c__(Bacteria) o__(Bacteria) f__(Bacteria) g__(Bacteria) s__(Bacteria) ID2 .0 k__Bacteria p__(Bacteria) c__(Bacteria) o__(Bacteria) f__(Bacteria) g__(Bacteria) s__(Bacteria) ID 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 .0 k__Bacteria p__Bacteroidetes c__Bacteroidia o__(Bacteroidia) f__(Bacteroidia) g__(Bacteroidia) s__(Bacteroidia) ID 0.0 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__(Bacteroidales) g__(Bacteroidales) s__(Bacteroidales) ID .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.E7 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides caccae ID47678 .0 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides caccae ID .0 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides dorei ID .0 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides dorei ID .0 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides dorei ID 85211.0 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides dorei ID .0 k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides fragilis ID817 .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 .0 

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

MEGAN分析进阶

双样本比对

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

数据质控

#重压缩测序文件 pigz -p 6 -dc ./ERR_1.fastq.gz | pigz -p 12 > ./S1_1.fastq.gz pigz -p 6 -dc ./ERR_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 

@CreationDate Tue Feb 04 17:49:13 CST 2020 @ContentType Summary4 @Names C1 S1 C11 @BlastMode BlastX BlastX BlastX @Uids 10 40 10 @Sizes .0 .0 .0 @TotalReads  @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  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文件的比对

物种信息可视化

这里我向大家介绍几种图形：；

• 堆叠华夫饼图，多个华夫饼图替换了气泡图的气泡，使用数量来代表丰富，但是必须在尺度范围内进行限制，目前R语言没有成熟的工具实现这个图形：
• 物种分布词云，这种方式让我们迅速就可以找到丰度较高的物种，或者差异物种。

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

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

功能信息可视化

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

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

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-muBTjnl4-58)(https://note.youdao.com/yws/public/resource/d39e91e4f7df7cf/xmlnote/BC4D049479C14A73B8CAB2D17A/36805)]

数据导出

无论是物种注释，还是功能注释，最终的结果我们当然需要进行一个很好的保存。MEGAN的保存类型很多，这里咱们来学习：

• 点击：Save As保存为.megan文件，这个文件保存后我们就可以随时打开查看数据，再也不需要使用rma文件。
• 点击Export Image将图片保存
• 点击Export Legend保存图例
• 红色选框种可以保存物种或者功能的表格文文件
• Metadate保存文件名
• Tree保存传统的数文件，这个文件我们可以使用其他树可视化软件出图，比如：Graphlan，ggtree，iTOL等。
  
  

[外链图片转存失败,源站可能有防盗链机制,建议将图片保存下来直接上传(img-zRmlpyDR-60)(https://note.youdao.com/yws/public/resource/d39e91e4f7df7cf/xmlnote/3B51DA0D8B8510A0A03BCE1626/36807)]

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 88 37 52 64 60 27 32 14 09 50 42 07 @Sizes .0 92730.0 63883.0 64549.0 57390.0 .0 .0 .0 .0 .0 .0 .0 @TotalReads  @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  743.0 TAX  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可视化终端，上面我们已经讲的很明白了。大家请参照。

撰文：文涛 南京农业大学

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

猜你喜欢

• 10000+: 菌群分析
  宝宝与猫狗 提DNA发Nature 实验分析谁对结果影响大 Cell微生物专刊 肠道指挥大脑
  
• 系列教程：微生物组入门 Biostar 微生物组 宏基因组
• 专业技能：生信宝典 学术图表 高分文章 不可或缺的人
• 一文读懂：宏基因组 寄生虫益处 进化树
• 必备技能：提问 搜索 Endnote
• 文献阅读 热心肠 SemanticScholar Geenmedical
• 扩增子分析：图表解读 分析流程 统计绘图
• 16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun
• 在线工具：16S预测培养基 生信绘图
• 科研经验：云笔记 云协作 公众号
• 编程模板: Shell R Perl
• 生物科普: 肠道细菌 人体上的生命 生命大跃进 细胞暗战 人体奥秘

写在后面