@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.htmlMEGAN提供了三种版本:Win MEGAN_Community_windows-x64_6_18_4.exeMAC MEGAN_Community_macos_6_18_4.dmgLinux MEGAN_Community_unix_6_18_4.sh如Linux版下载# 安装程序 102Mwget -c https://software-ab.informatik.uni-tuebingen.de/download/megan6/MEGAN_Community_unix_6_18_4.sh# NCBI-nr编号与物种和功能注释 970Mwget -c https://software-ab.informatik.uni-tuebingen.de/download/megan6/megan-map-Oct2019.db.zip# 核酸与物种信息 655MBwget -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(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.shbash MEGAN_Community_unix_6_18_4.sh# JVM must be at least 11. Please define INSTALL4J_JAVA_HOME to point to a suitable JVMjava -version # 1.8.0_201# 安装完上面conda会变为 openjdk version "11.0.1-internal" 2018-10-16
数据库下载# NCBI-nr编号与物种和功能注释 970Mwget -c https://software-ab.informatik.uni-tuebingen.de/download/megan6/megan-map-Oct2019.db.zip# 解压后为 5 Gb 的db文件unzip megan-map-Oct2019.db.zip# 核酸与物种信息 655MBwget -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,67Gwget -c ftp://ftp.ncbi.nlm.nih.gov/blast/db/FASTA/nr.gz#尝试使用pigz多线程解压缩,123Gtime unpigz -k -p 16 nr.gz # 8m, 26m#gunzip -c nr.gz > nrtime diamond makedb --in nr -d nr -p 32 # 8m, 102m
MEGAN使用指南(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.gzwget -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.gzpigz -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.daa10% 20% 30% 40% 50% 60% 70% 80% 90% Parsing file: ./diamond/C1.2.daa10% 20% 30% 40% 50% 60% 70% 80% 90% 100% (810.8s)Total reads: 443,738Alignments: 9,103,355100% (0.0s)100% (0.0s)100% (0.0s)Linking paired readsNumber of pairs: 178,830Binning reads: Initializing...Initializing binning...Using paired reads in taxonomic assignment...Using 'Naive LCA' algorithm for binning: TaxonomyUsing Best-Hit algorithm for binning: SEEDUsing Best-Hit algorithm for binning: EGGNOGUsing Best-Hit algorithm for binning: INTERPRO2GOBinning reads...Binning reads: Analyzing alignmentsTotal reads: 443,738With hits: 443,738Alignments: 9,103,355Assig. Taxonomy: 443,738Assig. SEED: 53,014Assig. EGGNOG: 85,333Assig. INTERPRO2GO: 204,180MinSupport set to: 221Binning reads: Applying min-support & disabled filter to Taxonomy...Min-supp. changes: 1,599Binning reads: Writing classification tablesNumb. Tax. classes: 92Numb. SEED classes: 3,055Numb. EGG. classes: 4,756Numb. INT. classes: 6,832Binning reads: SyncingClass. Taxonomy: 92Class. SEED: 3,055Class. EGGNOG: 4,756Class. INTERPRO2GO: 6,832100% (728.6s)Total time: 1547sPeak 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-S21EGGNOG:功能注释;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 GB8GB内存可以设置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 Bob06k__(null) p__(null) c__(null) o__(null) f__(null) g__(null) s__(null) ID1 462699.0k__(null) p__(null) c__(null) o__(null) f__(null) g__(null) s__(null) ID131567 80876.0k__Bacteria p__(Bacteria) c__(Bacteria) o__(Bacteria) f__(Bacteria) g__(Bacteria) s__(Bacteria) ID2 2493265.0k__Bacteria p__(Bacteria) c__(Bacteria) o__(Bacteria) f__(Bacteria) g__(Bacteria) s__(Bacteria) ID1783270 0.0k__Bacteria p__(Bacteria) c__(Bacteria) o__(Bacteria) f__(Bacteria) g__(Bacteria) s__(Bacteria) ID68336 6045.0k__Bacteria p__Bacteroidetes c__(Bacteroidetes) o__(Bacteroidetes) f__(Bacteroidetes) g__(Bacteroidetes) s__(Bacteroidetes) ID976 721448.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__(Bacteroidia) f__(Bacteroidia) g__(Bacteroidia) s__(Bacteroidia) ID200643 0.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__(Bacteroidales) g__(Bacteroidales) s__(Bacteroidales) ID171549 7600487.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__(Bacteroidaceae) s__(Bacteroidaceae) ID815 0.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__(Bacteroides) ID816 2.9878604E7k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides caccae ID47678 301910.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides caccae ID411901 100571.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides dorei ID357276 914270.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides dorei ID556260 101864.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides dorei ID997877 85211.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides dorei ID483217 225193.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides fragilis ID817 142742.0k__Bacteria p__Bacteroidetes c__Bacteroidia o__Bacteroidales f__Bacteroidaceae g__Bacteroides s__Bacteroides ovatus ID28116 97339.0k__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.gzpigz -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@ColorEditsTAX 1 1761.0 2274.0TAX 2 40023.0 80016.0TAX 1437201 46377.0 78736.0TAX 5125 5877.0 8092.0·······INTERPRO2GO 16380 17.0 16.0INTERPRO2GO 16381 1.0 0.0INTERPRO2GO 16382 0.0 2.0END_OF_DATA_TABLE#SampleID @SourceC1 ./C1.rmaS1 ./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.gzdo# 提取双端公共文件名,并输出检验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.logdone
使用nr数据库进行比对,使用34个线程跑,每个测序样本压缩文件只有10M,但是还是需要耗费将近50min一个一共12个样本,所以这个步骤耗时将近半天#建立比对文件夹,不见文件夹会--daa选项会提示错误cd ../mkdir ./diamond# 调用for循环批处理文件for filename in ./seq/_1.fq.gzdo# 提取双端公共文件名,并输出检验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.daadone
使用MEGAN提取结果mkdir ./diamond/MEGAN_RMA# 调用for循环批处理文件for filename in ./seq/_1.fq.gzdo# 提取双端公共文件名,并输出检验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}.rmadone
合并全部的物种和功能数据:#运行命令 ~/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@ColorEditsTAX -2 1.0TAX 1 342.0 736.0 685.0 278.0 260.0 449.0 477.0 568.0 161.0 409.0 381.0 574.0TAX 2049 0.0 134.0 188.0 580.0 41.0 188.0 0.0 0.0 0.0 228.0 14.0 588.0TAX 2 45382.0 33905.0 13366.0 16258.0 17194.0 40821.0 35587.0 38908.0 41649.0 31388.0 28311.0 37730.0TAX 143361 743.0TAX 1437201 0.0 0.0 8.0 4.0 5.0INTERPRO2GO 16381 8.0 1.0 0.0 0.0 0.0 0.0 0.0 0.0 1.0 0.0 28.0 5.0END_OF_DATA_TABLE#SampleID @Sourcep136C ./diamond/MEGAN_RMA/p136C.rmap136N ./diamond/MEGAN_RMA/p136N.rmap143C ./diamond/MEGAN_RMA/p143C.rmap143N ./diamond/MEGAN_RMA/p143N.rmap144C ./diamond/MEGAN_RMA/p144C.rmap144N ./diamond/MEGAN_RMA/p144N.rmap146C ./diamond/MEGAN_RMA/p146C.rmap146N ./diamond/MEGAN_RMA/p146N.rmap153C ./diamond/MEGAN_RMA/p153C.rmap153N ./diamond/MEGAN_RMA/p153N.rmap156C ./diamond/MEGAN_RMA/p156C.rmap156N ./diamond/MEGAN_RMA/p156N.rma
多样本批量物种和功能导出for filename in ./diamond/MEGAN_RMA//.rmado# 提取双端公共文件名,并输出检验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.txtdone
可视化和后续分析参考MEGAN可视化终端,上面我们已经讲的很明白了大家请参照撰文:文涛 南京农业大学责编:刘永鑫 中科院遗传发育所(图片来源网络,侵删)
0 评论