forked from YongxinLiu/EasyMetagenome
-
Notifications
You must be signed in to change notification settings - Fork 0
/
Copy path2pipeline.sh
1443 lines (1107 loc) · 57.7 KB
/
2pipeline.sh
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
246
247
248
249
250
251
252
253
254
255
256
257
258
259
260
261
262
263
264
265
266
267
268
269
270
271
272
273
274
275
276
277
278
279
280
281
282
283
284
285
286
287
288
289
290
291
292
293
294
295
296
297
298
299
300
301
302
303
304
305
306
307
308
309
310
311
312
313
314
315
316
317
318
319
320
321
322
323
324
325
326
327
328
329
330
331
332
333
334
335
336
337
338
339
340
341
342
343
344
345
346
347
348
349
350
351
352
353
354
355
356
357
358
359
360
361
362
363
364
365
366
367
368
369
370
371
372
373
374
375
376
377
378
379
380
381
382
383
384
385
386
387
388
389
390
391
392
393
394
395
396
397
398
399
400
401
402
403
404
405
406
407
408
409
410
411
412
413
414
415
416
417
418
419
420
421
422
423
424
425
426
427
428
429
430
431
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
467
468
469
470
471
472
473
474
475
476
477
478
479
480
481
482
483
484
485
486
487
488
489
490
491
492
493
494
495
496
497
498
499
500
501
502
503
504
505
506
507
508
509
510
511
512
513
514
515
516
517
518
519
520
521
522
523
524
525
526
527
528
529
530
531
532
533
534
535
536
537
538
539
540
541
542
543
544
545
546
547
548
549
550
551
552
553
554
555
556
557
558
559
560
561
562
563
564
565
566
567
568
569
570
571
572
573
574
575
576
577
578
579
580
581
582
583
584
585
586
587
588
589
590
591
592
593
594
595
596
597
598
599
600
601
602
603
604
605
606
607
608
609
610
611
612
613
614
615
616
617
618
619
620
621
622
623
624
625
626
627
628
629
630
631
632
633
634
635
636
637
638
639
640
641
642
643
644
645
646
647
648
649
650
651
652
653
654
655
656
657
658
659
660
661
662
663
664
665
666
667
668
669
670
671
672
673
674
675
676
677
678
679
680
681
682
683
684
685
686
687
688
689
690
691
692
693
694
695
696
697
698
699
700
701
702
703
704
705
706
707
708
709
710
711
712
713
714
715
716
717
718
719
720
721
722
723
724
725
726
727
728
729
730
731
732
733
734
735
736
737
738
739
740
741
742
743
744
745
746
747
748
749
750
751
752
753
754
755
756
757
758
759
760
761
762
763
764
765
766
767
768
769
770
771
772
773
774
775
776
777
778
779
780
781
782
783
784
785
786
787
788
789
790
791
792
793
794
795
796
797
798
799
800
801
802
803
804
805
806
807
808
809
810
811
812
813
814
815
816
817
818
819
820
821
822
823
824
825
826
827
828
829
830
831
832
833
834
835
836
837
838
839
840
841
842
843
844
845
846
847
848
849
850
851
852
853
854
855
856
857
858
859
860
861
862
863
864
865
866
867
868
869
870
871
872
873
874
875
876
877
878
879
880
881
882
883
884
885
886
887
888
889
890
891
892
893
894
895
896
897
898
899
900
901
902
903
904
905
906
907
908
909
910
911
912
913
914
915
916
917
918
919
920
921
922
923
924
925
926
927
928
929
930
931
932
933
934
935
936
937
938
939
940
941
942
943
944
945
946
947
948
949
950
951
952
953
954
955
956
957
958
959
960
961
962
963
964
965
966
967
968
969
970
971
972
973
974
975
976
977
978
979
980
981
982
983
984
985
986
987
988
989
990
991
992
993
994
995
996
997
998
999
1000
[TOC]
# 宏基因组分析流程 Pipeline of metagenomic analysis
- 版本 Version：1.10, 2021/1/22
- 系统要求 System: Linux Ubuntu 16.04+ / CentOS 7+
- 依赖软件 Sofware: (Rstudio server 1.1+)、KneadData v0.7.4、MetaPhlAn v2.7.5、HUMAnN2 v2.8.1、......
- Chrome浏览器访问Rstudio服务器
# 一、数据预处理 Data preprocessing
## 1.1 准备工作 Prepare
- 首次使用请参照`1soft_db.sh`脚本，安装软件和数据库(大约3-5天)
- EasyMetagenome项目Github或服务器家目录(~)有项目文件夹meta,
1. 包含测序数据(seq/*.fq.gz)和和实验设计(result/metadata.txt); 若无可手动从备份U盘或之前项目上传；
2. 使用谷歌浏览器访问RStudio服务器，内部服务器具体IP地址课上通知；
3. 以下段落代码新建项目，右侧File中进入项目目录并上传流程文件pipeline.sh，再单击打开。
### 1.1.1 环境变量设置(每次开始分析前必须运行)
设置数据库、软件和工作目录
# 公共数据库database(db)位置，如管理员设置/db，个人下载至~/db
db=/db
# Conda软件software安装目录，`conda env list`命令查看，如~/miniconda2
soft=/conda2
# 设置工作目录work directory(wd)，如meta
wd=~/meta
# 创建并进入工作目录
mkdir -p $wd && cd $wd
# 方法1.加载自己安装环境
# conda activate metagenome_env
# 方法2.加载别人安装环境
export PATH=${soft}/envs/metagenome_env/bin/:$PATH
source ${soft}/bin/activate metagenome_env
# 指定某个R语言环境
alias Rscript="/anaconda2/bin/Rscript --vanilla"
### 1.1.2 起始文件——序列和元数据
创建3个常用子目录：序列，临时文件和结果
mkdir -p seq temp result
# 上传实验设计(元数据) metadata.txt 于结果result目录
wget http://210.75.224.110/temp/meta/meta2010/result/metadata.txt -O result/metadata.txt
# 检查文件格式，^I为制表符，$为Linux换行，^M$为Windows回车，^M为Mac换行符
cat -A result/metadata.txt
# 转换Windows回车为Linux换行
sed -i 's/\r//' result/metadata.txt
cat -A result/metadata.txt
用户使用filezilla上传测序文件至seq目录，本次从其它位置复制，或从网络下载测试数据
# 方法1. 从其它目录复制测序数据
cp -rf /db/meta2101/seq/*.gz seq/
# 方法2. 网络下载测试数据
cd seq/
awk '{system("wget -c http://210.75.224.110/temp/meta/meta2010/seq/"$1"_1.fq.gz")}' \
<(tail -n+2 ../result/metadata.txt)
awk '{system("wget -c http://210.75.224.110/temp/meta/meta2010/seq/"$1"_2.fq.gz")}' \
<(tail -n+2 ../result/metadata.txt)
cd ..
# 查看文件大小
ls -lsh seq
# -l 列出详细信息 (l: list)
# -sh 显示人类可读方式文件大小 (s: size; h: human readable)
### 1.1.3 了解工作目录和文件
显示文件结构
tree -L 2
# .
# ├── pipeline.sh
# ├── result
# │ └── metadata.txt
# ├── seq
# │ ├── C1_1.fq.gz
# │ ├── C1_2.fq.gz
# │ ├── N1_1.fq.gz
# │ └── N1_2.fq.gz
# └── temp
- pipeline.sh是分析流程代码；
- seq目录中有12个样本Illumina双端测序，24个序列文件；
- temp是临时文件夹，存储分析中间文件，结束可全部删除节约空间
- result是重要节点文件和整理化的分析结果图表，
- 实验设计metadata.txt也在此
## 1.2 (可选)FastQC质量评估
# 第一次使用软件要记录软件版本，文章方法中必须写清楚
fastqc --version # 0.11.8
# time统计运行时间，此处耗时1m，fastqc质量评估
# *.gz为原始数据，-t指定多线程
time fastqc seq/*.gz -t 2
质控报告见`seq`目录，详细解读请阅读[《数据的质量控制软件——FastQC》](https://mp.weixin.qq.com/s/tDMih7ISLJcL4F4sWBq3Vw)。
multiqc将fastqc的多个报告生成单个整合报告，方法批量查看和比较
# 记录软件版本
multiqc --version # 1.5
# 整理seq目录下fastqc报告，输出multiqc_report.html至result/qc目录
multiqc -d seq/ -o result/qc
查看右侧result/qc目录中multiqc_report.html，单击，选择`View in Web Browser`查看可交互式报告。
## 1.3 KneadData质控和去宿主
kneaddata是流程，它主要依赖trimmomatic质控和去接头，bowtie2比对宿主，然后筛选非宿主序列用于下游分析 。
# 记录核心软件版本
kneaddata --version # 0.6.1
trimmomatic -version # 0.39
bowtie2 --version # 2.3.5
# 可只选一行中部分代码点击Run，如选中下行中#号后面命令查看程序帮助
# kneaddata -h # 显示帮助
检查点：zless/zcat查看压缩测序数据，检查序列质量格式(质量值为大写字母为标准Phred33格式 ，小写字母可能可能有Phred64需要使用fastp转换)；检查双端序列ID是否重复，如果重名需要在质控前改名更正。参考**附录，质控kneaddata，去宿主后双端不匹配——序列改名**。
# 设置某个样本名为变量i，以后再无需修改
i=C1
# zless查看压缩文件，空格翻页，按q退出。
zless seq/${i}_1.fq.gz | head -n4
# zcat显示压缩文件，head指定显示行数
zcat seq/${i}_2.fq.gz | head -n4
- | 为管道符，上一个命令的输出，传递给下一个命令做输入
- gzip: stdout: Broken pipe：管道断开。这里是人为断开，不是错误
- 运行过程中需要仔细阅读屏幕输出的信息
### 1.3.1 (可选)单样品质控
若一条代码分割在多行显示时，最好全部选中运行，多行分割的代码行末有一个 “\” 。多行注释命令运行，可全选，按Ctrl+Shift+C进行注释的取消和添加。
- 以metadata中`C1`样品本质控为例
1. 输入文件：双端FASTQ测序数据，提供给参数-i，seq/${i}_1.fq.gz和 seq/${i}_2.fq.gz
2. 参考数据库：宿主基因组索引 -db ${db}/kneaddata/human_genome/Homo_sapiens
3. 输出文件：质控后的FASTQ测序数据，在目录temp/qc下面，${i}_1_kneaddata_paired_1.fastq和${i}_1_kneaddata_paired_1.fastq，用于后续分析
4. 软件位置：`conda env list`查看软件安装位置，请务必根据自己软件和数据库安装位置，修改软件trimmomatic和接头文件位置。
(可选)手动设置trimmomatic程序和接头位置
如：${soft}/envs/metagenome_env/share/trimmomatic/
# 查看multiqc结果中接头污染最严重的样本，再到fastqc报告中查看接头序列，复制前20个碱基检索确定接头文件
grep 'GATCGGAAGAGCACACGTCT' ${soft}/envs/metagenome_env/share/trimmomatic/adapters/*
# 根据实际情况选择单端SE或双端PE，一般为 TruSeq3-PE-2.fa，更准确的是问测序公司要接头文件
100万条序列8线程质控和去宿主，耗时~2m。
mkdir -p temp/qc
time kneaddata -i seq/${i}_1.fq.gz -i seq/${i}_2.fq.gz \
-o temp/qc -v -t 8 --remove-intermediate-output \
--trimmomatic /conda2/envs/metagenome_env/share/trimmomatic/ \
--trimmomatic-options 'ILLUMINACLIP:/conda2/envs/metagenome_env/share/trimmomatic/adapters/TruSeq3-PE.fa:2:40:15 SLIDINGWINDOW:4:20 MINLEN:50' \
--bowtie2-options '--very-sensitive --dovetail' -db ${db}/kneaddata/human_genome/Homo_sapiens
# 查看质控后的结果文件
ls -shtr temp/qc/${i}_1_kneaddata_paired_?.fastq
### 1.3.1 多样品并行质控
并行队列管理软件——“parallel”。
# 记录软件版本
parallel --version # 20190522
# 打will cite承诺引用并行软件parallel
parallel --citation
parallel软件说明和使用实例
# 根据样本列表`:::`并行处理，并行j=2个任务，每个任务t=3个线程，2~7m
# 运行下面这行，体会下parallel的工作原理
# ::: 表示传递参数；第一个::: 后面为第一组参数，对应于{1};
# 第二个::: 后面为第二组参数，对应于{2}，依次替换
parallel -j 3 --xapply "echo {1} {2}" ::: seq/*_1.fq.gz ::: seq/*_2.fq.gz
# --xapply保持文件成对，否则将为两两组合，显示如下：
parallel -j 2 "echo {1} {2}" ::: seq/*_1.fq.gz ::: seq/*_2.fq.gz
# 从文件列表使用
parallel -j 3 --xapply "echo seq/{1}_1.fq.gz seq/{1}_2.fq.gz" ::: `tail -n+2 result/metadata.txt|cut -f1`
并行质控和去宿主：单样本运行成功，且参数设置绝对路径。出现错误`Unrecognized option: -d64`参考**附录，质控Kneaddata，Java版本不匹配——重装Java运行环境**。
# 每步分析产生多个文件时新建子文件夹
mkdir -p temp/qc
# 每个线程处理百万序列约10分钟，多线程可加速 j x t 倍
time parallel -j 2 --xapply \
"kneaddata -i seq/{1}_1.fq.gz \
-i seq/{1}_2.fq.gz \
-o temp/qc -v -t 9 --remove-intermediate-output \
--trimmomatic ${soft}/envs/metagenome_env/share/trimmomatic/ \
--trimmomatic-options 'ILLUMINACLIP:TruSeq3-PE.fa:2:40:15 SLIDINGWINDOW:4:20 MINLEN:50' \
--bowtie2-options '--very-sensitive --dovetail' \
-db ${db}/kneaddata/human_genome/Homo_sapiens" \
::: `tail -n+3 result/metadata.txt|cut -f1`
# (可选)大文件清理，高宿主含量样本可节约>90%空间
rm -rf temp/qc/*contam* temp/qc/*unmatched*
质控结果汇总
# 采用kneaddata附属工具kneaddata_read_count_table
kneaddata_read_count_table \
--input temp/qc \
--output temp/kneaddata.txt
# 筛选重点结果列
cut -f 1,2,4,12,13 temp/kneaddata.txt | sed 's/_1_kneaddata//' > result/qc/sum.txt
cat result/qc/sum.txt
# 用R代码统计下质控结果
Rscript -e "data=read.table('result/qc/sum.txt', header=T, row.names=1, sep='\t'); summary(data)"
# R转换宽表格为长表格
Rscript -e "library(reshape2); data=read.table('result/qc/sum.txt', header=T,row.names=1, sep='\t'); write.table(melt(data), file='result/qc/sum_long.txt',sep='\t', quote=F, col.names=T, row.names=F)"
cat result/qc/sum_long.txt
# 可用 http://www.ehbio.com/ImageGP/ 绘图展示
## 1.4 (可选)质控后质量评估
整理bowtie2, trimmomatic, fastqc报告，接头和PCR污染率一般小于1%。结果见：result/qc/multiqc_report_1.html
# 1-2m
fastqc temp/qc/*_1_kneaddata_paired_*.fastq -t 4
multiqc -d temp/qc/ -o result/qc/
# v1.7以后开始使用Python3，v1.9缺少bowtie2比对结果的统计
# 二、基于读长分析 Read-based (HUMAnN2)
中文教程：https://mp.weixin.qq.com/s/XkfT5MAo96KgyyVaN_Fl7g
英文教程：https://github.com/LangilleLab/microbiome_helper/wiki/Metagenomics-Tutorial-(Humann2)
## 2.1 准备HUMAnN2输入文件
小技巧：循环批量处理样本列表
# 基于样本元数据提取样本列表命令解析
# 去掉表头
tail -n+2 result/metadata.txt
# 提取第一列样本名
tail -n+2 result/metadata.txt|cut -f 1
# 循环处理样本
for i in `tail -n+2 result/metadata.txt|cut -f1`;do echo "Processing "$i; done
# ` 反引号为键盘左上角Esc键下面的按键，一般在数字1的左边，代表运行命令返回结果
HUMAnN2要求双端序列合并的文件作为输入，for循环根据实验设计样本名批量双端序列合并。注意星号和问号，分别代表多个和单个字符，当然大家更不能溜号~~~行分割的代码行末有一个 \
mkdir -p temp/concat
# 双端合并为单个文件
for i in `tail -n+2 result/metadata.txt|cut -f1`;do
cat temp/qc/${i}*_1_kneaddata_paired_?.fastq \
> temp/concat/${i}_all.fq;
done
# 取子集，如统一取为6G，1.6亿行，此处取20万行演示流程
for i in `tail -n+2 result/metadata.txt|cut -f1`;do
head -n800000 temp/concat/${i}_all.fq > temp/concat/${i}.fq
done
# 查看样品数量和大小
ls -sh temp/concat/*.fq
## 2.2 HUMAnN2计算物种和功能组成
- 物种组成调用MetaPhlAn2, bowtie2比对至核酸序列，解决有哪些微生物存在的问题；
- 功能组成为humann2调用diamond比对至蛋白库11Gb，解决这些微生物参与哪些功能通路的问题；
- 输入文件：temp/concat/*.fq 每个样品质控后双端合并后的fastq序列
- 输出文件：temp/humann2/ 目录下
- C1_pathabundance.tsv
- C1_pathcoverage.tsv
- C1_genefamilies.tsv
- 整合后的输出：
- result/metaphlan2/taxonomy.tsv 物种丰度表
- result/metaphlan2/taxonomy.spf 物种丰度表（用于stamp分析）
- result/humann2/pathabundance_relab_unstratified.tsv 通路丰度表
- result/humann2/pathabundance_relab_stratified.tsv 通路物种组成丰度表
- stratified(每个菌对此功能通路组成的贡献)和unstratified(功能组成)
(可选)启动humann2环境：仅humann2布置于自定义环境下使用
# 方法1. conda加载环境
# conda activate humann2
# 方法2. source加载指定
# source ~/miniconda3/envs/humann2/bin/activate
(可选)检查数据库配置是否正确
humann2 --version # v2.8.1/0.11.2
humann2_config
(可选)单样本1.25M PE150运行测试，8p，2.5M，1~2h；0.2M, 34m；0.1M，30m；0.01M，25m
mkdir -p temp/humann2
i=C1
time humann2 --input temp/concat/${i}.fq \
--output temp/humann2 --threads 8
多样本并行计算
mkdir -p temp/humann2
time parallel -j 2 \
'humann2 --input temp/concat/{}.fq \
--output temp/humann2/ --threads 8 ' \
::: `tail -n+2 result/metadata.txt|cut -f1` > temp/log
cat temp/log
# (可选)大文件清理，humann2临时文件可达原始数据30~40倍
# 链接重要文件至humann2目录
for i in `tail -n+2 result/metadata.txt|cut -f1`;do
ln temp/humann2/${i}_humann2_temp/${i}_metaphlan_bugs_list.tsv temp/humann2/
done
# 删除临时文件
# rm -rf temp/concat/* temp/humann2/*_humann2_temp
## 2.3 物种组成表
mkdir -p result/metaphlan2
### 2.3.1 样品结果合并
# 合并、修正样本名、预览
merge_metaphlan_tables.py \
temp/humann2/*_metaphlan_bugs_list.tsv | \
sed 's/_metaphlan_bugs_list//g' \
> result/metaphlan2/taxonomy.tsv
head -n3 result/metaphlan2/taxonomy.tsv
# (可选)筛选>0.5%的分类，绘制维恩图
awk 'BEGIN{OFS=FS="\t"}{if(FNR==1) {for(i=2;i<=NF;i++) a[i]=$i;} \
else {for(i=2;i<=NF;i++) if($i>0.5) print $1, a[i];}}' \
result/metaphlan2/taxonomy.tsv \
> result/metaphlan2/taxonomy_high.tsv
wc -l result/metaphlan2/taxonomy_high.tsv
### 2.3.2 转换为stamp的spf格式
metaphlan_to_stamp.pl result/metaphlan2/taxonomy.tsv \
> result/metaphlan2/taxonomy.spf
head -n3 result/metaphlan2/taxonomy.spf
# 下载metadata.txt和taxonomy.spf使用stamp分析
# 网络分析见附录 metaphlan2 - 网络
### 2.3.3 (可选) Python绘制热图
# c设置颜色方案，top设置物种数量，minv最小相对丰度，s标准化方法，log为取10为底对数，xy为势图宽和高，图片可选pdf/png/svg格式
metaphlan_hclust_heatmap.py \
--in result/metaphlan2/taxonomy.tsv \
--out result/metaphlan2/heatmap.pdf \
-c jet --top 30 --minv 0.1 \
-s log -x 0.4 -y 0.2
# 帮助见 metaphlan_hclust_heatmap.py -h
### 2.3.4 (可选) R绘制热图
# 推荐在Windows中运行
# 若在服务器上，请保持工作目录不变，并修改R脚本目录位置
# 本地分析目录
# 若在服务器上，这一句不需运行，已注释掉
# cd /c/meta
# 调置脚本所有目录，服务器位于 /db/script/ 或 ~/db/script，windows为/c/meta/db/script
R=/db/script/
# 显示脚本帮助 help
Rscript ${R}/metaphlan_hclust_heatmap.R -h
# 按指定列合并、排序并取Top25种绘制热图
# -i输入MetaPhlAn2文件；-t 分类级别，-t 分类级别，可选Kingdom/Phylum/Class/Order/Family/Genus/Species/Strain，界门纲目科属种株，推荐门，目，属
# -n 输出物种数量，默认为25，最大为合并后的数量
# -o输出图表前缀，默认根据输入文件、物种级别和数量自动生成
# 科水平Top25热图
Rscript ${R}/metaphlan_hclust_heatmap.R \
-i result/metaphlan2/taxonomy.spf \
-t Family -n 25 \
-w 183 -e 118 \
-o result/metaphlan2/heatmap_Family
## 2.4 功能组成分析
mkdir -p result/humann2
### 2.4.1 功能组成合并、标准化和分层
合并通路丰度(pathabundance)，含功能和对应物种组成。可选基因家族(genefamilies 太多)，通路覆盖度(pathcoverage)。注意看屏幕输出`# Gene table created: result/humann2/pathabundance.tsv`
humann2_join_tables \
--input temp/humann2 \
--file_name pathabundance \
--output result/humann2/pathabundance.tsv
# 样本名调整：删除列名多余信息
head result/humann2/pathabundance.tsv
sed -i 's/_Abundance//g' result/humann2/pathabundance.tsv
head result/humann2/pathabundance.tsv
标准化为相对丰度relab(1)或百万比cpm
humann2_renorm_table \
--input result/humann2/pathabundance.tsv \
--units relab \
--output result/humann2/pathabundance_relab.tsv
head -n5 result/humann2/pathabundance_relab.tsv
分层结果：功能和对应物种表(stratified)和功能组成表(unstratified)
humann2_split_stratified_table \
--input result/humann2/pathabundance_relab.tsv \
--output result/humann2/
# 可以使用stamp进行统计分析
### 2.4.2 差异比较和柱状图
两样本无法组间比较，在pcl层面替换为HMP数据进行统计和可视化。
参考 https://bitbucket.org/biobakery/humann2/wiki/Home#markdown-header-standard-workflow
- 输入数据：通路丰度表格 result/humann2/pathabundance.tsv
- 输入数据：实验设计信息 result/metadata.txt
- 中间数据：包含分组信息的通路丰度表格文件 result/humann2/pathabundance.pcl
- 输出结果：result/humann2/associate.txt
在通路丰度中添加分组
## 提取样品列表
head -n1 result/humann2/pathabundance.tsv | sed 's/# Pathway/SampleID/' | tr '\t' '\n' > temp/header
## 对应分组
awk 'BEGIN{FS=OFS="\t"}NR==FNR{a[$1]=$2}NR>FNR{print a[$1]}' result/metadata.txt temp/header | tr '\n' '\t'|sed 's/\t$/\n/' > temp/group
# 合成样本、分组+数据
cat <(head -n1 result/humann2/pathabundance.tsv) temp/group <(tail -n+2 result/humann2/pathabundance.tsv) \
> result/humann2/pathabundance.pcl
组间比较，样本量少无差异，结果为4列的文件：通路名字，通路在各个分组的丰度，差异P-value，校正后的Q-value。演示数据2样本无法统计，此处替换为HMP的结果演示统计和绘图(上传hmp_pathabund.pcl，替换pathabundance.pcl为hmp_pathabund.pcl)。
wget -c http://210.75.224.110/temp/meta/meta2010/result/humann2/hmp_pathabund.pcl -O result/humann2/hmp_pathabund.pcl
# 设置输入文件名
pcl=result/humann2/hmp_pathabund.pcl
# 统计表格行、列数量
csvtk -t stat ${pcl}
# 按分组KW检验
humann2_associate --input ${pcl} \
--focal-metadatum Group --focal-type categorical \
--last-metadatum Group --fdr 0.05 \
--output result/humann2/associate.txt
wc -l result/humann2/associate.txt
head -n5 result/humann2/associate.txt
barplot展示腺苷核苷酸合成的物种组成
# --sort sum metadata 按丰度和分组排序
# 指定差异通路，如 P163-PWY / PWY-3781 / PWY66-409 / PWY1F-823
path=PWY-3781
humann2_barplot --sort sum metadata \
--input ${pcl} \
--focal-feature ${path} \
--focal-metadatum Group --last-metadatum Group \
--output result/humann2/barplot_${path}.pdf
### 2.4.3 转换为KEGG注释
需要下载utility_mapping数据库并配置成功才可以使用。详见软件和数据库安装1soft_db.sh。
支持GO、PFAM、eggNOG、level4ec、KEGG的D级KO等注释，详见`humann2_regroup_table -h`。
# 转换基因家族为KO(uniref90_ko)，可选eggNOG(uniref90_eggnog)或酶(uniref90_level4ec)
for i in `tail -n+2 result/metadata.txt|cut -f1`;do
humann2_regroup_table \
-i temp/humann2/${i}_genefamilies.tsv \
-g uniref90_ko \
-o temp/humann2/${i}_ko.tsv
done
# 合并，并修正样本名
humann2_join_tables \
--input temp/humann2/ \
--file_name ko \
--output result/humann2/ko.tsv
sed -i '1s/_Abundance-RPKs//g' result/humann2/ko.tsv
## 2.5 GraPhlAn图
# metaphlan2 to graphlan
export2graphlan.py --skip_rows 1,2 -i result/metaphlan2/taxonomy.tsv \
--tree temp/merged_abundance.tree.txt \
--annotation temp/merged_abundance.annot.txt \
--most_abundant 1000 --abundance_threshold 20 --least_biomarkers 10 \
--annotations 3,4 --external_annotations 7
# 参数说明见PPT，或运行 export2graphlan.py --help
# graphlan annotation
graphlan_annotate.py --annot temp/merged_abundance.annot.txt \
temp/merged_abundance.tree.txt temp/merged_abundance.xml
# output PDF figure, annoat and legend
graphlan.py temp/merged_abundance.xml result/metaphlan2/graphlan.pdf \
--external_legends
## 2.6 LEfSe差异分析和Cladogram
- 输入文件：物种丰度表result/metaphlan2/taxonomy.tsv
- 输入文件：样品分组信息 result/metadata.txt
- 中间文件：整合后用于LefSe分析的文件 result/metaphlan2/lefse.txt，这个文件可以提供给www.ehbio.com/ImageGP 用于在线LefSE分析
- LefSe结果输出：result/metaphlan2/目录下lefse开头和feature开头的文件
准备输入文件，修改样本品为组名，**非通用代码**，可手动修改
# 只保留组名C, N, 删除数字重复编号，去除注释行
head -n3 result/metaphlan2/taxonomy.tsv
sed '1 s/[0-9]*//g' result/metaphlan2/taxonomy.tsv \
| grep -v '#' > result/metaphlan2/lefse.txt
head -n3 result/metaphlan2/lefse.txt
LEfSe本地代码供参考，可选Xshell下运行或ImageGP在线分析
# 格式转换为lefse内部格式
lefse-format_input.py \
result/metaphlan2/lefse.txt \
temp/input.in -c 1 -o 1000000
# 运行lefse
run_lefse.py temp/input.in \
temp/input.res
# 绘制物种树注释差异
lefse-plot_cladogram.py temp/input.res \
result/metaphlan2/lefse_cladogram.pdf --format pdf
# 绘制所有差异features柱状图
lefse-plot_res.py temp/input.res \
result/metaphlan2/lefse_res.pdf --format pdf
# 绘制单个features柱状图
# 查看显著差异features，按丰度排序
grep -v '-' temp/input.res | sort -k3,3n
# 手动选择指定feature绘图，如Firmicutes
lefse-plot_features.py -f one \
--feature_name "k__Bacteria.p__Firmicutes" \
--format pdf \
temp/input.in temp/input.res \
result/metaphlan2/lefse_Firmicutes.pdf
# 批量绘制所有差异features柱状图
lefse-plot_features.py -f diff \
--archive none --format pdf \
temp/input.in temp/input.res \
result/metaphlan2/lefse_
若出现错误 Unknown property axis_bgcolor，则修改`lefse-plot_cladogram.py`里的`ax_bgcolor`替换成`facecolor`即可。
## 2.7 Kraken2物种注释
Kraken2可以快速完成读长(read)层面的物种注释，还可以进行重叠群(contig)、基因(gene)、宏基因组组装基因组(MAG/bin)层面的序列物种注释。
# 启动kraken2工作环境
source /conda2/envs/kraken2/bin/activate
# 记录软件版本
kraken2 --version # 2.1.1
### 2.7.1 Kraken2物种注释
- {1}代表样本名字
- 输入：temp/qc/{1}_1_kneaddata_paired*.fastq 质控后的FASTQ数据
- 参考数据库：-db ${db}/kraken2/mini/，默认标准数据库>50GB，这里使用8GB迷你数据库。
- 输出结果：每个样本单独输出，temp/kraken2/{1}_report和temp/kraken2/{1}_output
- 整合后的输出结果： result/kraken2/taxonomy_count.txt 物种丰度表
(可选) 单样本注释，5m
mkdir -p temp/kraken2
i=C1
# 1m，--use-mpa-style可直接输出metphlan格式，但bracken无法处理
time kraken2 --db ${db}/kraken2/mini/ --paired temp/qc/${i}_1_kneaddata_paired*.fastq \
--threads 8 --use-names --report-zero-counts \
--report temp/kraken2/${i}.report \
--output temp/kraken2/${i}.output
多样本并行生成report，2样本各8线程
time parallel -j 2 \
"kraken2 --db ${db}/kraken2/mini --paired temp/qc/{1}_1_kneaddata_paired*.fastq \
--threads 8 --use-names --report-zero-counts \
--report temp/kraken2/{1}.report \
--output temp/kraken2/{1}.output" \
::: `tail -n+2 result/metadata.txt|cut -f1`
使用krakentools转换report为mpa格式
for i in `tail -n+2 result/metadata.txt|cut -f1`;do
kreport2mpa.py -r temp/kraken2/${i}.report \
--display-header \
-o temp/kraken2/${i}.mpa
done
合并样本为表格
mkdir -p result/kraken2
# 输出结果行数相同，但不一定顺序一致，要重新排序
parallel -j 1 \
'tail -n+2 temp/kraken2/{1}.mpa | sort | cut -f 2 | sed "1 s/^/{1}\n/" > temp/kraken2/{1}_count ' \
::: `tail -n+2 result/metadata.txt|cut -f1`
# 提取第一样本品行名为表行名
header=`tail -n 1 result/metadata.txt | cut -f 1`
echo $header
tail -n+2 temp/kraken2/${header}.mpa | sort | cut -f 1 | \
sed "1 s/^/Taxonomy\n/" > temp/kraken2/0header_count
head -n3 temp/kraken2/0header_count
# paste合并样本为表格
ls temp/kraken2/*count
paste temp/kraken2/*count > result/kraken2/tax_count.mpa
### 2.7.2 Bracken估计丰度
参数简介：
- -d为数据库，与kraken2一致
- -i为输入kraken2报告文件
- r是读长，此处默认为100，通常为150
- l为分类级，本次种级别(S)丰度估计，可选域、门、纲、目、科、属、种：D,P,C,O,F,G,S
- t是阈值，默认为0，越大越可靠，但可用数据越少
- -o 输出重新估计的值
循环重新估计每个样品的丰度
# 设置估算的分类级别D,P,C,O,F,G,S
tax=P
mkdir -p temp/bracken
for i in `tail -n+2 result/metadata.txt|cut -f1`;do
# i=C1
bracken -d ${db}/kraken2/mini \
-i temp/kraken2/${i}.report \
-r 100 -l ${tax} -t 0 \
-o temp/bracken/${i}
done
结果描述：共7列，分别为物种名、ID、分类级、读长计数、补充读长计数、**总数、总数百分比**
name taxonomy_id taxonomy_lvl kraken_assigned_reads added_reads new_est_reads fraction_total_reads
Capnocytophaga sputigena 1019 S 4798 996 5794 0.23041
Capnocytophaga sp. oral taxon 878 1316596 S 239 21 260 0.01034
bracken结果合并成表
# 输出结果行数相同，但不一定顺序一致，要去年表头重新排序
# 仅提取第6列reads count，并添加样本名
parallel -j 1 \
'tail -n+2 temp/bracken/{1}|sort|cut -f 6| sed "1 s/^/{1}\n/" > temp/bracken/{1}.count ' \
::: `tail -n+2 result/metadata.txt|cut -f1`
# 提取第一样本品行名为表行名
h=`tail -n1 result/metadata.txt|cut -f1`
tail -n+2 temp/bracken/${h}|sort|cut -f1 | \
sed "1 s/^/Taxonomy\n/" > temp/bracken/0header.count
# 检查文件数，为n+1
ls temp/bracken/*count | wc
# paste合并样本为表格，并删除非零行
paste temp/bracken/*count > result/kraken2/bracken.${tax}.txt
# 统计行列，默认去除表头
csvtk -t stat result/kraken2/bracken.${tax}.txt
结果筛选
# microbiome_helper按频率过滤，-r可标准化，-e过滤
Rscript /db/script/filter_feature_table.R \
-i result/kraken2/bracken.${tax}.txt \
-p 0.01 \
-o result/kraken2/bracken.${tax}.0.01
# > 0.01(1%)的样本在出现，剩1391行
csvtk -t stat result/kraken2/bracken.${tax}.0.01
# 按物种名手动去除宿主污染，以人为例
# 种水平去除人类P:Chordata,S:Homo sapiens
grep 'Homo sapiens' result/kraken2/bracken.S.0.01
grep -v 'Homo sapiens' result/kraken2/bracken.S.0.01 \
> result/kraken2/bracken.S.0.01-H
# 门水平去除脊索动物
grep 'Chordata' result/kraken2/bracken.P.0.01
grep -v 'Chordata' result/kraken2/bracken.P.0.01 > result/kraken2/bracken.P.0.01-H
分析后清理每条序列的注释大文件
rm -rf temp/kraken2/*.output
### 2.7.3 多样性分析
# 设置脚本位置
sd=/db/script
# 提取种级别注释并抽平至最小测序量，计算6种alpha多样性指数
# 查看帮助
Rscript $sd/otutab_rare.R -h
# -d指定最小样本量，默认0为最小值，抽平文件bracken.S.norm，alpha多样性bracken.S.alpha
tax=S
Rscript $sd/otutab_rare.R \
--input result/kraken2/bracken.${tax}.txt \
--depth 0 --seed 1 \
--normalize result/kraken2/bracken.${tax}.norm \
--output result/kraken2/bracken.${tax}.alpha
# 绘制Alpha多样性指数箱线图，详见script/MetaStatPlot.sh
# Beta多样性距离矩阵计算
mkdir -p result/kraken2/beta/
usearch -beta_div result/kraken2/bracken.${tax}.norm \
-filename_prefix result/kraken2/beta/
# PCoA分析，详见script/MetaStatPlot.sh
### 2.7.4 物种组成
- script/MetaStatPlot.sh ### 2.7.4 物种组成 (bracken2结果的物种组成堆叠可视化)
- script/MetaStatPlot.sh 附录 ## kraken2的物种组成热图和箱线图
# 三、组装分析流程 Assemble-based
## 3.1 拼接 Assembly
### 3.1.1 MEGAHIT拼接
# 删除旧文件夹，否则megahit无法运行
rm -rf temp/megahit
# 组装，10~30m，TB级数据需几天至几周
time megahit -t 3 \
-1 `tail -n+2 result/metadata.txt|cut -f 1|sed 's/^/temp\/qc\//;s/$/_1_kneaddata_paired_1.fastq/'| tr '\n' ','|sed 's/,$//'` \
-2 `tail -n+2 result/metadata.txt|cut -f 1|sed 's/^/temp\/qc\//;s/$/_1_kneaddata_paired_2.fastq/'| tr '\n' ','|sed 's/,$//'` \
-o temp/megahit
# 检查点：查看拼接结果，8.2M，通常300M~5G
ls -sh temp/megahit/final.contigs.fa
# 预览重叠群最前6行，前60列字符
head -n6 temp/megahit/final.contigs.fa | cut -c1-60
# 备份重要结果
mkdir -p result/megahit/
ln -f temp/megahit/final.contigs.fa result/megahit/
# 删除临时文件
rm -rf temp/megahit/intermediate_contigs/
### 3.1.2 (可选) metaSPAdes精细拼接
# 精细但使用内存和时间更多，15~65m
time metaspades.py -t 6 -m 100 \
`tail -n+2 result/metadata.txt|cut -f 1|sed 's/^/temp\/qc\//;s/$/_1_kneaddata_paired_1.fastq/'|sed 's/^/-1 /'| tr '\n' ' '` \
`tail -n+2 result/metadata.txt|cut -f 1|sed 's/^/temp\/qc\//;s/$/_1_kneaddata_paired_2.fastq/'|sed 's/^/-2 /'| tr '\n' ' '` \
-o temp/metaspades
# 23M，contigs体积更大
ls -sh temp/metaspades/contigs.fasta
# 备份重要结果
mkdir -p result/metaspades/
ln -f temp/metaspades/contigs.fasta result/metaspades/
# 删除临时文件
rm -rf temp/metaspades
### 3.1.3 QUAST评估
quast.py temp/megahit/final.contigs.fa -o result/megahit/quast -t 2
# 生成report文本tsv/txt、网页html、PDF等格式报告
# (可选) megahit和metaspades比较
quast.py --label "megahit,metapasdes" \
temp/megahit/final.contigs.fa \
temp/metaspades/contigs.fasta \
-o temp/quast
# (可选)metaquast评估，更全面，但需下载相关数据库，受网速影响可能时间很长
# metaquast based on silva, and top 50 species genome to access
time metaquast.py result/megahit/final.contigs.fa -o result/megahit/metaquast
## 3.2 基因预测、去冗余和定量
# Gene prediction, cluster & quantitfy
### 3.2.1 metaProdigal基因预测
# 输入文件：拼装好的序列文件 result/megahit/final.contigs.fa
# 输出文件：prodigal预测的基因序列 temp/prodigal/gene.fa
mkdir -p temp/prodigal
# prodigal的meta模式预测基因，35s，>和2>&1记录分析过程至gene.log
time prodigal -i result/megahit/final.contigs.fa \
-d temp/prodigal/gene.fa \
-o temp/prodigal/gene.gff \
-p meta -f gff > temp/prodigal/gene.log 2>&1
# 查看日志是否运行完成，有无错误
tail temp/prodigal/gene.log
# 统计基因数量19221
grep -c '>' temp/prodigal/gene.fa
### 3.2.2 基因聚类/去冗余cd-hit
# 输入文件：prodigal预测的基因序列 temp/prodigal/gene.fa
# 输出文件：去冗余后的基因和蛋白序列：result/NR/nucleotide.fa
# result/NR/protein.fa
mkdir -p result/NR
# aS覆盖度，c相似度，G局部比对，g最优解，T多线程，M内存0不限制
# 2万基因2m，2千万需要2000h，多线程可加速
time cd-hit-est -i temp/prodigal/gene.fa \
-o result/NR/nucleotide.fa \
-aS 0.9 -c 0.95 -G 0 -g 0 -T 0 -M 0
# 统计非冗余基因数量19182，单次拼接结果数量下降不大，多批拼接冗余度高
grep -c '>' result/NR/nucleotide.fa
# 翻译核酸为对应蛋白序列，emboss
transeq -sequence result/NR/nucleotide.fa \
-outseq result/NR/protein.fa -trim Y
# 序列名自动添加了_1，为与核酸对应要去除
sed -i 's/_1 / /' result/NR/protein.fa
### 3.2.3 基因定量salmon
# 输入文件：去冗余后的基因和蛋白序列：result/NR/nucleotide.fa
# 输出文件：Salmon定量后的结果：result/salmon/gene.count
# result/salmon/gene.TPM
mkdir -p temp/salmon
salmon -v # 1.3.0
# 建索引, -t序列, -i 索引，10s
time salmon index \
-t result/NR/nucleotide.fa \
-p 9 \
-i temp/salmon/index
# 定量，l文库类型自动选择，p线程，--meta宏基因组模式, 2个任务并行12个样，共48s
# 注意parallel中待并行的命令必须是双引号
time parallel -j 2 \
"salmon quant \
-i temp/salmon/index -l A -p 8 --meta \
-1 temp/qc/{1}_1_kneaddata_paired_1.fastq \
-2 temp/qc/{1}_1_kneaddata_paired_2.fastq \
-o temp/salmon/{1}.quant" \
::: `tail -n+2 result/metadata.txt|cut -f1`
# 合并
mkdir -p result/salmon
salmon quantmerge \
--quants temp/salmon/*.quant \
-o result/salmon/gene.TPM
salmon quantmerge \
--quants temp/salmon/*.quant \
--column NumReads -o result/salmon/gene.count
sed -i '1 s/.quant//g' result/salmon/gene.*
# 预览结果表格
head -n3 result/salmon/gene.*
## 3.3 功能基因注释
# 输入数据：上一步预测的蛋白序列 result/NR/protein.fa
# 中间结果：temp/eggnog/protein.emapper.seed_orthologs
# temp/eggnog/output.emapper.annotations
# temp/eggnog/output
# COG定量表：result/eggnog/cogtab.count
# result/eggnog/cogtab.count.spf (用于STAMP)
# KO定量表：result/eggnog/kotab.count
# result/eggnog/kotab.count.spf (用于STAMP)
# CAZy碳水化合物注释和定量：result/dbcan2/cazytab.count
# result/dbcan2/cazytab.count.spf (用于STAMP)
# 抗生素抗性：result/resfam/resfam.count
# result/resfam/resfam.count.spf (用于STAMP)
# 这部分可以拓展到其它数据库
### 3.3.1 基因注释eggNOG(COG/KEGG/CAZy)
# https://github.com/eggnogdb/eggnog-mapper/wiki/eggNOG-mapper-v2
# 记录软件版本
emapper.py --version # 2.0.0
# diamond比对基因至eggNOG 5.0数据库, 10h
mkdir -p temp/eggnog
time emapper.py --no_annot --no_file_comments --override \
--data_dir ${db}/eggnog5 \
-i result/NR/protein.fa \
--cpu 9 -m diamond \
-o temp/eggnog/protein
# 比对结果功能注释, 1h
time emapper.py \
--annotate_hits_table temp/eggnog/protein.emapper.seed_orthologs \
--data_dir ${db}/eggnog5 \
--cpu 9 --no_file_comments --override \
-o temp/eggnog/output
# 添表头, 1列为ID，9列KO，16列CAZy，21列COG，22列描述
sed '1 i Name\tortholog\tevalue\tscore\ttaxonomic\tprotein\tGO\tEC\tKO\tPathway\tModule\tReaction\trclass\tBRITE\tTC\tCAZy\tBiGG\ttax_scope\tOG\tbestOG\tCOG\tdescription' \
temp/eggnog/output.emapper.annotations \
> temp/eggnog/output
summarizeAbundance生成COG/KO/CAZy丰度汇总表
mkdir -p result/eggnog
# 显示帮助，需要Python3环境，可修改软件第一行指定python位置
summarizeAbundance.py -h
# 汇总，9列KO，16列CAZy按逗号分隔，21列COG按字母分隔，原始值累加
summarizeAbundance.py \
-i result/salmon/gene.TPM \
-m temp/eggnog/output \
-c '9,16,21' -s ',+,+*' -n raw \
-o result/eggnog/eggnog
sed -i 's/^ko://' result/eggnog/eggnog.KO.raw.txt
# eggnog.CAZy.raw.txt eggnog.COG.raw.txt eggnog.KO.raw.txt
# 添加注释生成STAMP的spf格式，结合metadata.txt进行差异比较
# KO，v2021.1
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \
${db}/kegg/KO_description.txt \
result/eggnog/eggnog.KO.raw.txt | \
sed 's/^\t/Unannotated\t/' \
> result/eggnog/eggnog.KO.TPM.spf
# CAZy，v2020.7
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \
${db}/dbcan2/CAZy_description.txt result/eggnog/eggnog.CAZy.raw.txt | \
sed 's/^\t/Unannotated\t/' > result/eggnog/eggnog.CAZy.TPM.spf
# COG
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2"\t"$3} NR>FNR{print a[$1],$0}' \
${db}/eggnog/COG.anno \
result/eggnog/eggnog.COG.raw.txt | sed '1 s/^/Level1\tLevel2/'> \
result/eggnog/eggnog.COG.TPM.spf
### 3.3.2 (可选)碳水化合物dbCAN2
# 比对CAZy数据库, 用时2~18m; 加--sensitive更全但慢至1h
mkdir -p temp/dbcan2
time diamond blastp \
--db ${db}/dbCAN2/CAZyDB.07312020 \
--query result/NR/protein.fa \
--threads 9 -e 1e-5 --outfmt 6 --max-target-seqs 1 --quiet \
--out temp/dbcan2/gene_diamond.f6
# 整理比对数据为表格
mkdir -p result/dbcan2
# 提取基因与dbcan分类对应表
perl /db/script/format_dbcan2list.pl \
-i temp/dbcan2/gene_diamond.f6 \
-o temp/dbcan2/gene.list
# 按对应表累计丰度
summarizeAbundance.py \
-i result/salmon/gene.TPM \
-m temp/dbcan2/gene.list \
-c 2 -s ',' -n raw \
-o result/dbcan2/TPM
# 添加注释生成STAMP的spf格式，结合metadata.txt进行差异比较
awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \
${db}/dbcan2/CAZy_description.txt result/dbcan2/TPM.CAZy.raw.txt | \
sed 's/^\t/Unannotated\t/' > result/dbcan2/TPM.CAZy.raw.spf
# 检查未注释数量，有则需要检查原因
grep 'Unannotated' result/dbcan2/TPM.CAZy.raw.spf|wc -l
### 3.3.3 抗生素抗性CARD
数据库：https://card.mcmaster.ca/
参考文献：http://doi.org/10.1093/nar/gkz935
软件使用Github: https://github.com/arpcard/rgi
# 启动rgi环境
source ${soft}/bin/activate rgi
# 加载数据库
rgi load -i ${db}/card/card.json \
--card_annotation ${db}/card/card_database_v3.1.0.fasta --local
# 蛋白注释
mkdir -p result/card
cut -f 1 -d ' ' result/NR/protein.fa > temp/protein.fa
# --include_loose 会显著增加上百倍结果
time rgi main \
--input_sequence temp/protein.fa \
--output_file result/card/protein \
--input_type protein --clean \
--num_threads 9 --alignment_tool DIAMOND
结果说明：
- protein.json，在线可视化
- protein.txt，注释基因列表
# 四、挖掘单菌基因组/分箱(Binning)
## 4.1 MetaWRAP
# 主要使用MetaWRAP，演示基于官方测试数据
# 主页：https://github.com/bxlab/metaWRAP
# 挖掘单菌基因组，需要研究对象复杂度越低、测序深度越大，结果质量越好。要求单样本6GB+，复杂样本如土壤推荐数据量30GB+，至少3个样本
# 上面的演示数据12个样仅140MB，无法获得单菌基因组，这里使用官方测序数据演示讲解
# 软件和数据库布置需2-3天，演示数据分析过程超10h，标准30G样也需3-30天，由服务器性能决定。
### 4.1.1 准备数据和环境变量
# 流程: https://github.com/bxlab/metaWRAP/blob/master/Usage_tutorial.md
>
# 输入数据：质控后的FASTQ序列，文件名格式必须为*_1.fastq和*_2.fastq