博文
宏基因组注释和可视化神器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 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官网下载页面：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
当然还有需要许可证的收费版本：数据库就包含了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, 102mMEGAN使用指南（Linux）
原始序列处理：如 https://www.ebi.ac.uk/ena/browser/view/ERR793599
# 下载
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.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提取功能：
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.txtWin版安装和使用
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.0KEGG注释查看：
"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.0EGGNOG注释查看：
"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.0519228E7InterPro2GO注释结果查看：
"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.0TAX物种注释信息查看：默认输出只有灭有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可视化终端，上面我们已经讲的很明白了。大家请参照。
撰文：文涛 南京农业大学
责编：刘永鑫 中科院遗传发育所
猜你喜欢
- 10000+: 菌群分析
宝宝与猫狗 提DNA发Nature 实验分析谁对结果影响大 Cell微生物专刊 肠道指挥大脑 - 系列教程：微生物组入门 Biostar 微生物组 宏基因组
- 专业技能：生信宝典 学术图表 高分文章 不可或缺的人
- 一文读懂：宏基因组 寄生虫益处 进化树
- 必备技能：提问 搜索 Endnote
- 文献阅读 热心肠 SemanticScholar Geenmedical
- 扩增子分析：图表解读 分析流程 统计绘图
- 16S功能预测 PICRUSt FAPROTAX Bugbase Tax4Fun
- 在线工具：16S预测培养基 生信绘图
- 科研经验：云笔记 云协作 公众号
- 编程模板: Shell R Perl
- 生物科普: 肠道细菌 人体上的生命 生命大跃进 细胞暗战 人体奥秘
写在后面
为鼓励读者交流、快速解决科研困难，我们建立了“宏基因组”专业讨论群，目前己有国内外5000+ 一线科研人员加入。参与讨论，获得专业解答，欢迎分享此文至朋友圈，并扫码加主编好友带你入群，务必备注“姓名-单位-研究方向-职称/年级”。技术问题寻求帮助，首先阅读《如何优雅的提问》学习解决问题思路，仍末解决群内讨论，问题不私聊，帮助同行。
学习扩增子、宏基因组科研思路和分析实战，关注“宏基因组”
点击阅读原文，跳转最新文章目录阅读
https://mp.weixin.qq.com/s/5jQspEvH5_4Xmart22gjMA
https://wap.sciencenet.cn/blog-3334560-1219163.html
上一篇：SBB：土壤微生物群落的特征究竟由什么决定
下一篇：Cell：植物根系如何允许有益微生物定植的
