-
Notifications
You must be signed in to change notification settings - Fork 104
/
pipeline_mac.sh
1478 lines (1178 loc) · 62.1 KB
/
pipeline_mac.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]
# 易扩增子EasyAmplicon
# 作者 Authors: 刘永鑫(Yong-Xin Liu), 陈同(Tong Chen)等
# 版本 Version: v1.21
# 更新 Update: 2024-4-12
# 系统要求 System requirement: Windows 10+ / Mac OS 10.12+ / Ubuntu 20.04+
# 引文 Reference: Liu, et al. 2023. EasyAmplicon: An easy-to-use, open-source, reproducible, and community-based
# pipeline for amplicon data analysis in microbiome research. iMeta 2: e83. https://doi.org/10.1002/imt2.83
# 设置工作(work directory, wd)和软件数据库(database, db)目录
# 添加环境变量，并进入工作目录 Add environmental variables and enter work directory
# **每次打开Rstudio必须运行下面4行 Run it**，可选替换${db}为EasyMicrobiome安装位置
wd=~/github/EasyAmplicon
db=~/github/EasyMicrobiome
# Mac 与 win 的设置不同
bin=${db}/mac
chmod 755 ${bin}/*
sed -i 's/\r//g' ${db}/script/*.sh
sed -i 's/\r//g' ${db}/script/*.R
sed -i 's/\r//g' ${db}/script/*.py
sed -i 's/\r//g' ${bin}/*.sh
sed -i 's/\r//g' ${bin}/*.r
chmod 755 ${db}/script/*
export PATH=${bin}:`pwd`/${bin}:${db}/script:${PATH}
cd ${wd}
## 0. Mac安装一些基本命令
# Mac是 Unix 系统，与 Linux 用起来大同小异，但有一些工具还是有差别
# 如果您之前没有做过配置，需要先运行下面这些命令，重新安装一些
# 常用工具和与 Linux 服务器相同的工具版本
# 以免出现命令不兼容
# /bin/bash -c "$(curl -fsSL https://raw.githubusercontent.com/Homebrew/install/HEAD/install.sh)"
/bin/zsh -c "$(curl -fsSL https://gitee.com/cunkai/HomebrewCN/raw/master/Homebrew.sh)"
brew install coreutils
brew install gawk
brew install gsed
## 1. 起始文件 start files
# 1. 分析流程pipeline.sh
# 2. 样本元信息metadata.txt，保存于result目录
# 3. 测序数据fastq文件保存于seq目录，通常以`.fq.gz`结尾，每个样品一对文件
# 4. 创建临时文件存储目录，分析结束可删除
# -p: 表示若目录不存在则新建，若存在则不做任何操作
mkdir -p temp
### 1.1. 元数据/实验设计 metadata
# 准备样本元数据result/metadata.txt
# csvtk统计表行(样本数，不含表头)列数，-t设置列分隔为制表符，默认为;
csvtk -t stat result/metadata_raw.txt
# 元数据至少3列，首列为样本ID(SampleID)，结尾列为描述(Description)
# cat查看文件，-A显示符号，"|"为管道符实现命令连用，head显示文件头，-n3控制范围前3行
# 如果报错，说明用的是 mac 的 cat，需要参数未 -e 的这一行
# cat -e result/metadata_raw.txt | head -n3
cat -A result/metadata_raw.txt | head -n3
# windows用户结尾有^M，运行sed命令去除，再用cat -A检查
sed 's/\r//' result/metadata_raw.txt > result/metadata.txt
cat -A result/metadata.txt | head -n3
### 1.2. 测序数据 sequencing data
# # 本段代码可在RStudio中Ctrl + Shift + C 取消注释“#”后运行
# # (可选)下载测序数据，按GSA的CRA(批次)和CRR(样品)编号下载数据
# # 示例下载单个文件并改名
# mkdir -p seq
# wget -c ftp://download.big.ac.cn/gsa/CRA002352/CRR117575/CRR117575_f1.fq.gz -O seq/KO1_1.fq.gz
# # 按实验设计编号批量下载并改名
# awk '{system("wget -c ftp://download.big.ac.cn/gsa/"$5"/"$6"/"$6"_f1.fq.gz -O seq/"$1"_1.fq.gz")}' \
# <(tail -n+2 result/metadata.txt)
# awk '{system("wget -c ftp://download.big.ac.cn/gsa/"$5"/"$6"/"$6"_r2.fq.gz -O seq/"$1"_2.fq.gz")}' \
# <(tail -n+2 result/metadata.txt)
# 公司返回的测序结果，通常为一个样品一对fq/fastq.gz格式压缩文件
# 文件名与样品名务必对应：不一致时手工修改，批量改名见"常见问题6"
# 如果测序数据是.gz的压缩文件，有时需要使用gunzip解压后使用，vsearch通常可直接读取压缩文件
# gunzip seq/*.gz
# zless按页查看压缩文件，空格翻页、q退出；head默认查看前10行，-n指定行
ls -sh seq/
zless seq/KO1_1.fq.gz | head -n4
# 每行太长，指定查看每行的1-60个字符
zless seq/KO1_1.fq | head | cut -c 1-60
# 统计测序数据，依赖seqkit程序
seqkit stat seq/KO1_1.fq.gz
# 批量统计测序数据并汇总表
seqkit stat seq/*.fq.gz > result/seqkit.txt
head result/seqkit.txt
### 1.3. 流程和数据库 pipeline & database
# 数据库第一次使用必须解压，以后可跳过此段
# usearchs可用16S/18S/ITS数据库：RDP, SILVA和UNITE，本地文件位置 ${db}/usearch/
# usearch数据库database下载页: http://www.drive5.com/usearch/manual/sintax_downloads.html
# 解压16S RDP数据库，gunzip解压缩，seqkit stat统计
gunzip -c ${db}/usearch/rdp_16s_v18.fa.gz > ${db}/usearch/rdp_16s_v18.fa
seqkit stat ${db}/usearch/rdp_16s_v18.fa # 2.1万条序列
# 解压ITS UNITE数据库，需自行从官网或网盘db/amplicon/usearch中下载
# gunzip -c ${db}/usearch/utax_reference_dataset_all_25.07.2023.fasta.gz >${db}/usearch/unite.fa
# seqkit stat ${db}/usearch/unite.fa # 32.6万
# Greengene数据库用于功能注释: ftp://greengenes.microbio.me/greengenes_release/gg_13_5/gg_13_8_otus.tar.gz
# 默认解压会删除原文件，-c指定输出至屏幕，> 写入新文件(可改名)
gunzip -c ${db}/gg/97_otus.fasta.gz > ${db}/gg/97_otus.fa
seqkit stat ${db}/gg/97_otus.fa
## 2. 序列合并和重命名 reads merge and rename
### 2.1 合并双端序列并按样品重命名 Merge pair-end reads and rename
#测试.以WT1单样品合并为例
#time统计计算时间，real是物理时间，user是计算时间，sys是硬件等待时间
# time vsearch --fastq_mergepairs seq/WT1_1.fq.gz \
# --reverse seq/WT1_2.fq.gz \
# --fastqout temp/WT1.merged.fq \
# --relabel WT1.
# head temp/WT1.merged.fq
#依照实验设计批处理并合并
#tail -n+2去表头，cut -f1取第一列，获得样本列表；18个样本x1.5万对序列合并8s
#Win下复制Ctrl+C为Linux下中止，为防止异常中断，结尾添加&转后台，无显示后按回车继续
# 一部分电脑 rush 不支持，运行时调度失败，请使用 for 循环部分
# for 循环部分是放入后台运行的，点完 run 之后，看上去程序已运行完，实际没运行完，而是正在运行中。
# 不要急于运行后面的程序。
# 之前课程，有发现每次运行结果都不一样，就是因为 for 循环部分没运行完，只生成了部分数据，导致后面
# 每个样品 reads 数每次运行都会不一致。
#方法1.for循环顺序处理
# time for i in `tail -n+2 result/metadata.txt|cut -f1`;do
# vsearch --fastq_mergepairs seq/${i}_1.fq.gz --reverse seq/${i}_2.fq.gz \
# --fastqout temp/${i}.merged.fq --relabel ${i}.
# done &
#方法2.rush并行处理，任务数jobs(j),2可加速1倍4s；建议设置2-4
# 部分 mac 可能不支持 gzip 版本的文件，可通过下面的命令查看
# vsearch --version
## vsearch v2.22.1_macos_x86_64, 8.0GB RAM, 8 cores
## https://github.com/torognes/vsearch
##
## Rognes T, Flouri T, Nichols B, Quince C, Mahe F (2016)
## VSEARCH: a versatile open source tool for metagenomics
## PeerJ 4:e2584 doi: 10.7717/peerj.2584 https://doi.org/10.7717/peerj.2584
##
## Compiled with support for gzip-compressed files, but the library was not found.
## Compiled with support for bzip2-compressed files, but the library was not found.
# 如果不支持 zlib，可以用这个方式运行。
time tail -n+2 result/metadata.txt | cut -f 1 | \
rush -j 2 "vsearch --fastq_mergepairs <(zcat seq/{}_1.fq.gz) --reverse <(zcat seq/{}_2.fq.gz) \
--fastqout temp/{}.merged.fq --relabel {}."
# time tail -n+2 result/metadata.txt | cut -f 1 | \
# rush -j 3 "vsearch --fastq_mergepairs seq/{}_1.fq.gz --reverse seq/{}_2.fq.gz \
# --fastqout temp/{}.merged.fq --relabel {}."
# 检查最后一个文件前10行中样本名
head temp/`tail -n+2 result/metadata.txt | cut -f 1 | tail -n1`.merged.fq | grep ^@
##方法3.不支持压缩文件时解压再双端合并
# time tail -n+2 result/metadata.txt | cut -f 1 | \
# rush -j 1 "vsearch --fastq_mergepairs <(zcat seq/{}_1.fq.gz) --reverse <(zcat seq/{}_2.fq.gz) \
# --fastqout temp/{}.merged.fq --relabel {}."
#
# time for i in `tail -n+2 result/metadata.txt|cut -f1`;do
# vsearch --fastq_mergepairs seq/${i}_1.fq --reverse seq/${i}_2.fq \
# --fastqout temp/${i}.merged.fq --relabel ${i}.
# done &
### 2.2 (可选)单端文件改名 Single-end reads rename
# # 单个序列改名示例
# i=WT1
# gunzip -c seq/${i}_1.fq.gz > seq/${i}.fq
# usearch -fastx_relabel seq/${i}.fq -fastqout temp/${i}.merged.fq -prefix ${i}.
#
# # 批量改名，需要有单端fastq文件，且解压(usearch不支持压缩格式)
# gunzip seq/*.gz
# time for i in `tail -n+2 result/metadata.txt|cut -f1`;do
# usearch -fastx_relabel seq/${i}.fq -fastqout temp/${i}.merged.fq -prefix ${i}.
# done &
# # vsearch大数据方法参考“常见问题2”
### 2.3 改名后序列整合 integrate renamed reads
#合并所有样品至同一文件
cat temp/*.merged.fq > temp/all.fq
#查看文件大小223M，软件不同版本结果略有差异
ls -lsh temp/all.fq
# 查看序列名，“.”之前是否为样本名，样本名绝不允许有点 (".")
# 样本名有点的一个显著特征是生成的特征表会很大，特征表里面列很多，导致后面分析出现内存不足。
# 后面分析获得特征表后要看一眼有没有问题，遇到内存不足问题，也要回头来排查。
head -n 6 temp/all.fq | cut -c1-60
## 3. 切除引物与质控 Cut primers and quality filter
# 左端10bp标签+19bp上游引物V5共为29，右端V7为18bp下游引物
# Cut barcode 10bp + V5 19bp in left and V7 18bp in right
# 务必清楚实验设计和引物长度，引物已经去除可填0，27万条序列14s
time vsearch --fastx_filter temp/all.fq \
--fastq_stripleft 29 --fastq_stripright 18 \
--fastq_maxee_rate 0.01 \
--fastaout temp/filtered.fa
# 查看文件了解fa文件格式
head temp/filtered.fa
## 4. 去冗余挑选OTU/ASV Dereplicate and cluster/denoise
### 4.1 序列去冗余 Dereplicate
# 并添加miniuniqusize最小为8或1/1M，去除低丰度噪音并增加计算速度
# -sizeout输出丰度, --relabel必须加序列前缀更规范, 1s
vsearch --derep_fulllength temp/filtered.fa \
--minuniquesize 10 --sizeout --relabel Uni_ \
--output temp/uniques.fa
#高丰度非冗余序列非常小(500K~5M较适合)，名称后有size和频率
ls -lsh temp/uniques.fa
# Uni_1;size=6423 - 去冗余后序列的名字 Uni_1；该序列在所有样品测序数据中出现 6423 次
# 为出现最多的序列。
head -n 2 temp/uniques.fa
### 4.2 聚类OTU/去噪ASV Cluster or denoise
#有两种方法：推荐unoise3去噪获得单碱基精度ASV，备选传统的97%聚类OTU(属水平精度)
#usearch两种特征挑选方法均自带de novo去嵌合体
#-minsize二次过滤，控制OTU/ASV数量至1-5千，方便下游统计分析
#方法1. 97%聚类OTU，适合大数据/ASV规律不明显/reviewer要求
#结果耗时1s, 产生508 OTUs, 去除126 chimeras
# usearch -cluster_otus temp/uniques.fa -minsize 10 \
# -otus temp/otus.fa \
# -relabel OTU_
#方法2. ASV去噪 Denoise: predict biological sequences and filter chimeras
#6s, 1530 good, 41 chimeras, 序列百万条可能需要几天/几周
# usearch -unoise3 temp/uniques.fa -minsize 10 \
# -zotus temp/zotus.fa
#修改序列名：Zotu为改为ASV方便识别
# sed 's/Zotu/ASV_/g' temp/zotus.fa > temp/otus.fa
# head -n 2 temp/otus.fa
#方法3. 数据过大无法使用usearch时，备选vsearch方法见"常见问题3"
vsearch --cluster_unoise temp/uniques.fa --centroids temp/otus1.fa --relabel ASV_
vsearch --uchime3_denovo temp/otus1.fa --nonchimeras temp/otus.fa
head -n 2 temp/otus.fa
### 4.3 基于参考去嵌合 Reference-based chimera detect
# 不推荐，容易引起假阴性，因为参考数据库无丰度信息
# 而de novo时要求亲本丰度为嵌合体16倍以上防止假阴性
# 因为已知序列不会被去除，数据库选择越大越合理，假阴性率最低
mkdir -p result/raw
# 方法1. vsearch+rdp去嵌合(快但容易假阴性)
# 可自行下载silva并解压(替换rdp_16s_v18.fa为silva_16s_v123.fa)，极慢但理论上更好
vsearch --uchime_ref temp/otus.fa \
-db ${db}/usearch/rdp_16s_v18.fa \
--nonchimeras result/raw/otus.fa
# RDP: 7s, 143 (9.3%) chimeras; SILVA：9m, 151 (8.4%) chimeras
# Win vsearch结果添加了windows换行符^M需删除，mac不要执行此命令
sed -i 's/\r//g' result/raw/otus.fa
# 方法2. 不去嵌合
# cp -f temp/otus.fa result/raw/otus.fa
## 5. 特征表构建和筛选 Feature table create and filter
# OTU和ASV统称为特征(Feature)，它们的区别是：
# OTU通常按97%聚类后挑选最高丰度或中心的代表性序列；
# ASV是基于序列进行去噪(排除或校正错误序列，并挑选丰度较高的可信序列)作为代表性序列
### 5.1 生成特征表
# 方法1. usearch生成特征表，小样本(<30)快；但大样本受限且多线程效率低，83.2%, 4核17s
# time usearch -otutab temp/filtered.fa \
# -otus result/raw/otus.fa \
# -threads 4 \
# -otutabout result/raw/otutab.txt
# 方法2. vsearch生成特征表
# id(1)：100%相似度比对49.45%序列，1m50s
# id(0.97)：97%相似度比对83.66%序列，1m10s(更高数据使用率，更快)
time vsearch --usearch_global temp/filtered.fa \
--db result/raw/otus.fa \
--id 0.97 --threads 4 \
--otutabout result/raw/otutab.txt
#212862 of 268019 (79.42%)可比对
# vsearch结果windows用户删除换行符^M校正为标准Linux格式
sed -i 's/\r//' result/raw/otutab.txt
head -n6 result/raw/otutab.txt | cut -f 1-6 |cat -A
# csvtk统计表行列
# 这里一定看好列数，是不是等于你的样品数；如果不等，一般是样品命名存在问题，具体看上面解释
csvtk -t stat result/raw/otutab.txt
### 5.2 物种注释，且/或去除质体和非细菌 Remove plastid and non-Bacteria
# 物种注释-去除质体和非细菌/古菌并统计比例(可选)
# RDP物种注释(rdp_16s_v18)更快，但缺少完整真核来源数据,可能不完整，耗时15s;
# SILVA数据库(silva_16s_v123.fa)更好注释真核、质体序列，极慢耗时3h起
# 置信阈值通常0.6/0.8，vserch最低0.1/usearch可选0输出最相似物种注释用于观察潜在分类
vsearch --sintax result/raw/otus.fa \
--db ${db}/usearch/rdp_16s_v18.fa \
--sintax_cutoff 0.1 \
--tabbedout result/raw/otus.sintax
head result/raw/otus.sintax | cat -A
sed -i 's/\r//' result/raw/otus.sintax
# 方法1. 原始特征表行数
wc -l result/raw/otutab.txt
#R脚本选择细菌古菌(真核)、去除叶绿体、线粒体并统计比例；输出筛选并排序的OTU表
#输入为OTU表result/raw/otutab.txt和物种注释result/raw/otus.sintax
#输出筛选并排序的特征表result/otutab.txt和
#统计污染比例文件result/raw/otutab_nonBac.txt和过滤细节otus.sintax.discard
#真菌ITS数据，请改用otutab_filter_nonFungi.R脚本，只筛选真菌
# Rscript ${db}/script/otutab_filter_nonBac.R -h # 显示参数说明
Rscript ${db}/script/otutab_filter_nonBac.R \
--input result/raw/otutab.txt \
--taxonomy result/raw/otus.sintax \
--output result/otutab.txt\
--stat result/raw/otutab_nonBac.stat \
--discard result/raw/otus.sintax.discard
# 筛选后特征表行数
wc -l result/otutab.txt
#过滤特征表对应序列
cut -f 1 result/otutab.txt | tail -n+2 > result/otutab.id
vsearch -fastx_getseqs result/raw/otus.fa \
-labels result/otutab.id \
-fastaout result/otus.fa
#过滤特征表对应序列注释
awk 'NR==FNR{a[$1]=$0}NR>FNR{print a[$1]}'\
result/raw/otus.sintax result/otutab.id \
> result/otus.sintax
#补齐末尾列
sed -i 's/\t$/\td:Unassigned/' result/otus.sintax
head -n2 result/otus.sintax
# 方法2. 觉得筛选不合理可以不筛选
# cp result/raw/otu* result/
#可选统计方法：OTU表简单统计 Summary OTUs table
# 因为 Mac 升级后用不了 usearch，额外写了个小程序处理
summary.sh -f result/otutab.txt -c CTCTCT
#注意最小值、分位数，或查看result/raw/otutab_nonBac.stat中样本详细数据量，用于重采样
### 5.3 等量抽样标准化
# Normlize by subsample
#使用vegan包进行等量重抽样，输入reads count格式Feature表result/otutab.txt
#可指定输入文件、抽样量和随机数，输出抽平表result/otutab_rare.txt和多样性alpha/vegan.txt
mkdir -p result/alpha
Rscript ${db}/script/otutab_rare.R --input result/otutab.txt \
--depth 10000 --seed 1 \
--normalize result/otutab_rare.txt \
--output result/alpha/vegan.txt
# usearch -otutab_stats result/otutab_rare.txt \
# -output result/otutab_rare.stat
# cat result/otutab_rare.stat
summary.sh -f result/otutab_rare.txt -c CTCTCT
## 6. α多样性 alpha diversity
### 6.1. 计算α多样性 calculate alpha diversity
# 使用USEARCH计算14种alpha多样性指数(Chao1有错勿用)
#details in http://www.drive5.com/usearch/manual/alpha_metrics.html
# usearch -alpha_div result/otutab_rare.txt \
# -output result/alpha/alpha.txt
# 或用下面的R脚本计算6种多样性值
compute_alpha.R --input result/otutab_rare.txt \
--depth 0 \
--output result/alpha/alpha.txt
### 6.2. 计算稀释丰富度 calculate rarefaction richness
#稀释曲线：取1%-100%的序列中OTUs数量，每次无放回抽样
#Rarefaction from 1%, 2% .. 100% in richness (observed OTUs)-method without_replacement https://drive5.com/usearch/manual/cmd_otutab_subsample.html
# 这步部分 mac 用不了，也意义不大；如果需要，自行购买收费版 usearch
# usearch -alpha_div_rare result/otutab_rare.txt \
# -output result/alpha/alpha_rare.txt \
# -method without_replacement
#预览结果
head -n2 result/alpha/alpha_rare.txt
#样本测序量低出现非数值"-"的处理，详见常见问题8
sed -i "s/-/\t0.0/g" result/alpha/alpha_rare.txt
#预览结果
head -n2 result/alpha/alpha_rare.txt
#样本测序量低出现非数值"-"的处理，详见常见问题8
sed -i "s/-/\t0.0/g" result/alpha/alpha_rare.txt
### 6.3. 筛选高丰度菌 Filter by abundance
#计算各特征的均值，有组再求分组均值，需根据实验设计metadata.txt修改组列名
#输入文件为feautre表result/otutab.txt，实验设计metadata.txt
#输出为特征表按组的均值-一个实验可能有多种分组方式
#-h显示脚本帮助(参数说明)
Rscript ${db}/script/otu_mean.R -h
#scale是否标准化，zoom标准化总和，all输出全部样本均值，type计算类型mean或sum
Rscript ${db}/script/otu_mean.R --input result/otutab.txt \
--metadata result/metadata.txt \
--group Group --thre 0 \
--scale TRUE --zoom 100 --all TRUE --type mean \
--output result/otutab_mean.txt
# 结果为全部和各组均值
head -n3 result/otutab_mean.txt
#如以平均丰度>0.1%筛选，可选0.5或0.05，得到每个组的OTU组合
awk 'BEGIN{OFS=FS="\t"}{if(FNR==1) {for(i=3;i<=NF;i++) a[i]=$i; print "OTU","Group";} \
else {for(i=3;i<=NF;i++) if($i>0.1) print $1, a[i];}}' \
result/otutab_mean.txt > result/alpha/otu_group_exist.txt
head result/alpha/otu_group_exist.txt
cut -f 2 result/alpha/otu_group_exist.txt | sort | uniq -c
# 试一试：不同丰度下各组有多少OTU/ASV
# 可在 http://ehbio.com/test/venn/ 中绘图并显示各组共有和特有维恩或网络图
# 也可在 https://www.bic.ac.cn/ImageGP 或 https://www.bic.ac.cn/BIC/#/ 绘制Venn、upSetView和Sanky
## 7. β多样性 Beta diversity
# 这一步 mac 用不了，可以用在线网站https://www.bic.ac.cn/BIC/
# 计算 beta 多样性+PCoA 可视化+统计检验一起
#结果有多个文件，需要目录
mkdir -p result/beta/
#基于OTU构建进化树 Make OTU tree, 4s
usearch -cluster_agg result/otus.fa -treeout result/otus.tree
#生成5种距离矩阵：bray_curtis, euclidean, jaccard, manhatten, unifrac
usearch -beta_div result/otutab_rare.txt -tree result/otus.tree \
-filename_prefix result/beta/
## 8. 物种注释分类汇总
#OTU对应物种注释2列格式：去除sintax中置信值，只保留物种注释，替换:为_，删除引号
cut -f 1,4 result/otus.sintax \
|sed 's/\td/\tk/;s/:/__/g;s/,/;/g;s/"//g' \
> result/taxonomy2.txt
head -n3 result/taxonomy2.txt
#OTU对应物种8列格式：注意注释是非整齐
#生成物种表格OTU/ASV中空白补齐为Unassigned
awk 'BEGIN{OFS=FS="\t"}{delete a; a["k"]="Unassigned";a["p"]="Unassigned";a["c"]="Unassigned";a["o"]="Unassigned";a["f"]="Unassigned";a["g"]="Unassigned";a["s"]="Unassigned";\
split($2,x,";");for(i in x){split(x[i],b,"__");a[b[1]]=b[2];} \
print $1,a["k"],a["p"],a["c"],a["o"],a["f"],a["g"],a["s"];}' \
result/taxonomy2.txt > temp/otus.tax
sed 's/;/\t/g;s/.__//g;' temp/otus.tax|cut -f 1-8 | \
sed '1 s/^/OTUID\tKingdom\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies\n/' \
> result/taxonomy.txt
head -n3 result/taxonomy.txt
#统计门纲目科属，使用 rank参数 p c o f g，为phylum, class, order, family, genus缩写
mkdir -p result/tax
for i in p c o f g;do
sintax_summary.R -t result/taxonomy.txt -i result/otutab_rare.txt --rank ${i} \
--output result/tax/sum_${i}.txt
done
sed -i 's/(//g;s/)//g;s/\"//g;s/\#//g;s/\/Chloroplast//g' result/tax/sum_*.txt
# 列出所有文件
wc -l result/tax/sum_*.txt
head -n3 result/tax/sum_g.txt
## 9. 有参定量特征表
# 比对Greengenes97% OTUs比对，用于PICRUSt/Bugbase功能预测
mkdir -p result/gg/
#方法1. usearch比对更快，但文件超限报错选方法2
# 默认10核以下使用1核，10核以上使用10核
# usearch -otutab temp/filtered.fa -otus ${db}/gg/97_otus.fa \
# -otutabout result/gg/otutab.txt -threads 4
# 比对率80.0%, 1核11m，4核3m，10核2m，内存使用743Mb
#head -n3 result/gg/otutab.txt
# #方法2. vsearch比对，更准更慢，但并行24-96线程更强
vsearch --usearch_global temp/filtered.fa --db ${db}/gg/97_otus.fa \
--otutabout result/gg/otutab.txt --id 0.97 --threads 8
# 比对率81.04%, 1核30m, 12核7m
head -n3 result/gg/otutab.txt
#统计
#usearch -otutab_stats result/gg/otutab.txt -output result/gg/otutab.stat
#cat result/gg/otutab.stat
summary.sh -f result/gg/otutab.txt -c CTCTCT
## 10. 空间清理及数据提交
#删除中间大文件
rm -rf temp/*.fq
# 分双端统计md5值，用于数据提交
cd seq
md5sum *_1.fq.gz > md5sum1.txt
md5sum *_2.fq.gz > md5sum2.txt
paste md5sum1.txt md5sum2.txt | awk '{print $2"\t"$1"\t"$4"\t"$3}' | sed 's/*//g' > ../result/md5sum.txt
rm md5sum*
cd ..
cat result/md5sum.txt
# R语言多样性和物种组成分析
## 1. Alpha多样性
### 1.1 Alpha多样性箱线图
# 在线绘图平台 https://www.bic.ac.cn/BIC/#/ 提供更多定制参数和绘制灵活性
# 查看帮助
Rscript ${db}/script/alpha_boxplot.R -h
# 完整参数，多样性指数可选richness chao1 ACE shannon simpson invsimpson
Rscript ${db}/script/alpha_boxplot.R --alpha_index richness \
--input result/alpha/vegan.txt --design result/metadata.txt \
--group Group --output result/alpha/ \
--width 89 --height 59
# 使用循环绘制6种常用指数
for i in `head -n1 result/alpha/vegan.txt|cut -f 2-`;do
Rscript ${db}/script/alpha_boxplot.R --alpha_index ${i} \
--input result/alpha/vegan.txt --design result/metadata.txt \
--group Group --output result/alpha/ \
--width 89 --height 59
done
mv alpha_boxplot_TukeyHSD.txt result/alpha/
# Alpha多样性柱状图+标准差
Rscript ${db}/script/alpha_barplot.R --alpha_index richness \
--input result/alpha/vegan.txt --design result/metadata.txt \
--group Group --output result/alpha/ \
--width 89 --height 59
### 1.2 稀释曲线
Rscript ${db}/script/alpha_rare_curve.R \
--input result/alpha/alpha_rare.txt --design result/metadata.txt \
--group Group --output result/alpha/ \
--width 120 --height 59
### 1.3 多样性维恩图
# 交互式 Venn 图 https://www.bic.ac.cn/test/venn
# 三组比较:-f输入文件,-a/b/c/d/g分组名,-w/u为宽高英寸,-p输出文件名后缀
bash ${db}/script/sp_vennDiagram.sh \
-f result/alpha/otu_group_exist.txt \
-a WT -b KO -c OE \
-w 3 -u 3 \
-p WT_KO_OE
# 四组比较，图和代码见输入文件目录，运行目录为当前项目根目录
bash ${db}/script/sp_vennDiagram.sh \
-f result/alpha/otu_group_exist.txt \
-a WT -b KO -c OE -d All \
-w 3 -u 3 \
-p WT_KO_OE_All
## 2. Beta多样性
### 2.1 距离矩阵热图pheatmap
# 添加分组注释，如2，4列的基因型和地点
cut -f 1-2 result/metadata.txt > temp/group.txt
# 以bray_curtis为例，-f输入文件,-h是否聚类TRUE/FALSE,-u/v为宽高英寸
# -P添加行注释文件，-Q添加列注释
bash ${db}/script/sp_pheatmap.sh \
-f result/beta/bray_curtis.txt \
-H 'TRUE' -u 6.9 -v 5.6 \
-P temp/group.txt -Q temp/group.txt
### 2.2 主坐标分析PCoA
# 在线绘图平台 https://www.bic.ac.cn/BIC/#/ 提供更多定制参数和绘制灵活性
# 输入文件，选择分组，输出文件，图片尺寸mm，统计见beta_pcoa_stat.txt
Rscript ${db}/script/beta_pcoa.R \
--input result/beta/bray_curtis.txt --design result/metadata.txt \
--group Group --label FALSE --width 89 --height 59 \
--output result/beta/bray_curtis.pcoa.pdf
# 添加样本标签 --label TRUE
Rscript ${db}/script/beta_pcoa.R \
--input result/beta/bray_curtis.txt --design result/metadata.txt \
--group Group --label TRUE --width 89 --height 59 \
--output result/beta/bray_curtis.pcoa.label.pdf
mv beta_pcoa_stat.txt result/beta/
### 2.3 限制性主坐标分析CPCoA
Rscript ${db}/script/beta_cpcoa.R \
--input result/beta/bray_curtis.txt --design result/metadata.txt \
--group Group --output result/beta/bray_curtis.cpcoa.pdf \
--width 89 --height 59
# 添加样本标签 --label TRUE
Rscript ${db}/script/beta_cpcoa.R \
--input result/beta/bray_curtis.txt --design result/metadata.txt \
--group Group --label TRUE --width 89 --height 59 \
--output result/beta/bray_curtis.cpcoa.label.pdf
## 3. 物种组成Taxonomy
### 3.1 堆叠柱状图Stackplot
# 在线绘图平台 https://www.bic.ac.cn/BIC/#/ 提供更多定制参数和绘制灵活性
# 以门(p)水平为例，结果包括output.sample/group.pdf两个文件
Rscript ${db}/script/tax_stackplot.R \
--input result/tax/sum_p.txt --design result/metadata.txt \
--group Group --color ggplot --legend 7 --width 89 --height 59 \
--output result/tax/sum_p.stackplot
# 修改颜色--color ggplot, manual1(22), Paired(12) or Set3(12)
Rscript ${db}/script/tax_stackplot.R \
--input result/tax/sum_p.txt --design result/metadata.txt \
--group Group --color Paired --legend 12 --width 181 --height 119 \
--output result/tax/sum_p.stackplotPaired
# 批量绘制输入包括p/c/o/f/g共5级
# https://www.bic.ac.cn/BIC 提供更多 top 物种筛选方式
for i in p c o f g; do
Rscript ${db}/script/tax_stackplot.R \
--input result/tax/sum_${i}.txt --design result/metadata.txt \
--group Group --output result/tax/sum_${i}.stackplot \
--legend 8 --width 89 --height 59; done
### 3.2 弦/圈图circlize
# 以纲(class,c)为例，绘制前5组
i=c
Rscript ${db}/script/tax_circlize.R \
--input result/tax/sum_${i}.txt --design result/metadata.txt \
--group Group --legend 5
# 结果位于当前目录circlize.pdf(随机颜色)，circlize_legend.pdf(指定颜色+图例)
# 移动并改名与分类级一致
mv circlize.pdf result/tax/sum_${i}.circlize.pdf
mv circlize_legend.pdf result/tax/sum_${i}.circlize_legend.pdf
### 3.3 树图treemap(参考)
# 多层级包含物种关系，输入特征表和物种注释，输出树图
# 指定包含特征数量和图片宽高，100个ASV耗时12s
Rscript ${db}/script/tax_maptree.R \
--input result/otutab.txt --taxonomy result/taxonomy.txt \
--output result/tax/tax_maptree.pdf \
--topN 100 --width 183 --height 118
# 24、差异比较
## 1. R语言差异分析
### 1.1 差异比较
# Error in file(file, ifelse(append, "a", "w")),输出目录不存在，创建目录即可
mkdir -p result/compare/
# 输入特征表、元数据；指定分组列名、比较组和丰度
# 选择方法 wilcox/t.test/edgeR、pvalue和fdr和输出目录
compare="KO-WT"
Rscript ${db}/script/compare.R \
--input result/otutab.txt --design result/metadata.txt \
--group Group --compare ${compare} --threshold 0.1 \
--method edgeR --pvalue 0.05 --fdr 0.2 \
--output result/compare/
### 1.2 火山图
# 在线绘图平台 https://www.bic.ac.cn/BIC/#/ 提供更多定制参数和绘制灵活性
# 输入compare.R的结果，输出火山图带数据标签，可指定图片大小
Rscript ${db}/script/compare_volcano.R \
--input result/compare/${compare}.txt \
--output result/compare/${compare}.volcano.pdf \
--width 89 --height 59
### 1.3 热图
# 输入compare.R的结果，筛选列数，指定元数据和分组、物种注释，图大小英寸和字号
bash ${db}/script/compare_heatmap.sh -i result/compare/${compare}.txt -l 7 \
-d result/metadata.txt -A Group \
-t result/taxonomy.txt \
-w 8 -h 5 -s 7 \
-o result/compare/${compare}
### 1.4 曼哈顿图
# i差异比较结果,t物种注释,p图例,w宽,v高,s字号,l图例最大值
# 图例显示不图，可增加高度v为119+即可，后期用AI拼图为KO-WT.heatmap.emf
bash ${db}/script/compare_manhattan.sh -i result/compare/${compare}.txt \
-t result/taxonomy.txt \
-p result/tax/sum_p.txt \
-w 183 -v 59 -s 7 -l 10 \
-o result/compare/${compare}.manhattan.p.pdf
# 上图只有6个门，切换为纲c和-L Class展示细节
bash ${db}/script/compare_manhattan.sh -i result/compare/${compare}.txt \
-t result/taxonomy.txt \
-p result/tax/sum_c.txt \
-w 183 -v 59 -s 7 -l 10 -L Class \
-o result/compare/${compare}.manhattan.c.pdf
# 显示完整图例，再用AI拼图
bash ${db}/script/compare_manhattan.sh -i result/compare/${compare}.txt \
-t result/taxonomy.txt \
-p result/tax/sum_c.txt \
-w 183 -v 149 -s 7 -l 10 -L Class \
-o result/compare/${compare}.manhattan.c.legend.pdf
### 1.5 单个特征的绘制
# 筛选显示差异ASV，按KO组丰度降序列，取ID展示前10
awk '$4<0.05' result/compare/KO-WT.txt | sort -k7,7nr | cut -f1 | head
# 差异OTU细节展示
Rscript ${db}/script/alpha_boxplot.R --alpha_index ASV_2 \
--input result/otutab.txt --design result/metadata.txt \
--transpose TRUE --scale TRUE \
--width 89 --height 59 \
--group Group --output result/compare/feature_
# ID不存在会报错： Error in data.frame(..., check.names = FALSE) : 参数值意味着不同的行数: 0, 18 Calls: alpha_boxplot -> cbind -> cbind -> data.frame
# 指定某列排序：按属丰度均值All降序
# Mac 下的输出是已经排序过的
# csvtk -t sort -k All:nr result/tax/sum_g.txt | head
# 差属细节展示
Rscript ${db}/script/alpha_boxplot.R --alpha_index Lysobacter \
--input result/tax/sum_g.txt --design result/metadata.txt \
--transpose TRUE \
--width 89 --height 59 \
--group Group --output result/compare/feature_
### 1.5 三元图
#参考示例见：result\compare\ternary\ternary.Rmd 文档
#备选教程[246.三元图的应用与绘图实战](https://mp.weixin.qq.com/s/3w3ncpwjQaMRtmIOtr2Jvw)
## 2. STAMP输入文件准备
### 2.1 生成输入文件
Rscript ${db}/script/format2stamp.R -h
mkdir -p result/stamp
Rscript ${db}/script/format2stamp.R --input result/otutab.txt \
--taxonomy result/taxonomy.txt --threshold 0.01 \
--output result/stamp/tax
# 可选Rmd文档见result/format2stamp.Rmd
### 2.2 绘制扩展柱状图和表
compare="KO-WT"
# 替换ASV(result/otutab.txt)为属(result/tax/sum_g.txt)
Rscript ${db}/script/compare_stamp.R \
--input result/stamp/tax_5Family.txt --metadata result/metadata.txt \
--group Group --compare ${compare} --threshold 0.1 \
--method "t.test" --pvalue 0.05 --fdr "none" \
--width 189 --height 159 \
--output result/stamp/${compare}
# 可选Rmd文档见result/CompareStamp.Rmd
## 3. LEfSe输入文件准备
### 3.1. 命令行生成文件
# 可选命令行生成输入文件
Rscript ${db}/script/format2lefse.R -h
mkdir -p result/lefse
# threshold控制丰度筛选以控制作图中的枝数量
Rscript ${db}/script/format2lefse.R --input result/otutab.txt \
--taxonomy result/taxonomy.txt --design result/metadata.txt \
--group Group --threshold 0.4 \
--output result/lefse/LEfSe
### 3.2 Rmd生成输入文件(可选)
#1. result目录中存在otutab.txt, metadata.txt, taxonomy.txt三个文件；
#2. Rstudio打开EasyAmplicon中format2lefse.Rmd，另存至result目录并Knit生成输入文件和可重复计算网页；
### 3.3 LEfSe分析
#方法1. 打开LEfSe.txt并在线提交 https://www.bic.ac.cn/BIC/#/analysis?tool_type=tool&page=b%27MzY%3D%27
#方法2. LEfSe本地分析(限Linux系统、选学)，参考代码见附录
#方法3. LEfSe官网在线使用
# 25、QIIME 2分析流程
# 代码详见 qiime2/pipeline_qiime2.sh
# 31、功能预测
## 1. PICRUSt功能预测
# PICRUSt 1.0
# 方法1. 使用 https://www.bic.ac.cn/ImageGP 在线分析 gg/otutab.txt
# 方法2. Linux服务器用户可参考"附录2. PICRUSt功能预测"实现软件安装和分析
# 然后结果使用STAMP/R进行差异比较
# R语言绘图
# 输入文件格式调整
l=L2
sed '/# Const/d;s/OTU //' result/picrust/all_level.ko.${l}.txt > result/picrust/${l}.txt
num=`head -n1 result/picrust/${l}.txt|wc -w`
cut -f $num result/picrust/${l}.txt >a1
cut -f 1-$[num-1] result/picrust/${l}.txt >a2
paste a1 a2 > result/picrust/${l}.spf
cut -f 2- result/picrust/${l}.spf > result/picrust/${l}.mat.txt
awk 'BEGIN{FS=OFS="\t"} {print $2,$1}' result/picrust/${l}.spf | \
sed 's/;/\t/' | sed '1 s/ID/Pathway\tCategory/' \
> result/picrust/${l}.anno.txt
head result/picrust/${l}.anno.txt
# 差异比较
compare="KO-WT"
Rscript ${db}/script/compare.R \
--input result/picrust/${l}.mat.txt --design result/metadata.txt \
--group Group --compare ${compare} --threshold 0 \
--method wilcox --pvalue 0.05 --fdr 0.2 \
--output result/picrust/
# 可对结果${compare}.txt筛选
# 绘制指定组(A/B)的柱状图，按高分类级着色和分面
Rscript ${db}/script/compare_hierarchy_facet.R \
--input result/picrust/${compare}.txt \
--data MeanA \
--annotation result/picrust/${l}.anno.txt \
--output result/picrust/${compare}.MeanA.bar.pdf
# 绘制两组显著差异柱状图，按高分类级分面
Rscript ${db}/script/compare_hierarchy_facet2.R \
--input result/picrust/${compare}.txt \
--pvalue 0.05 --fdr 0.1 \
--annotation result/picrust/${l}.anno.txt \
--output result/picrust/${compare}.bar.pdf
### PICRUSt 2.0
# 软件安装见附录6. PICRUSt环境导出和导入
# (可选)PICRUSt2(Linux/Windows下Linux子系统，要求>16GB内存)
# 安装参考附录5的方式直接下载安装包并解压即可使用
# Linux中加载conda环境
conda activate picrust2
# 进入工作目录，服务器要修改工作目录
wd=${wd}/result/picrust2
mkdir -p ${wd} && cd ${wd}
# 运行流程，内存15.7GB，耗时12m
picrust2_pipeline.py -s ../otus.fa -i ../otutab.txt -o ./out -p 8
# 添加EC/KO/Pathway注释
cd out
add_descriptions.py -i pathways_out/path_abun_unstrat.tsv.gz -m METACYC \
-o METACYC.tsv
add_descriptions.py -i EC_metagenome_out/pred_metagenome_unstrat.tsv.gz -m EC \
-o EC.tsv
add_descriptions.py -i KO_metagenome_out/pred_metagenome_unstrat.tsv.gz -m KO \
-o KO.tsv
# KEGG按层级合并
db=/mnt/d/EasyMicrobiome/
python3 ${db}/script/summarizeAbundance.py \
-i KO.tsv \
-m ${db}/kegg/KO1-4.txt \
-c 2,3,4 -s ',+,+,' -n raw \
-o KEGG
# 统计各层级特征数量
wc -l KEGG*
# 可视化见picrust2文件夹中ggpicrust2.Rmd
## 2. 元素循环FAPROTAX
## 方法1. 在线分析，推荐使用 https://www.bic.ac.cn/ImageGP 在线分析
## 方法2. Linux下分析、如QIIME 2环境下
# 设置工作目录
wd2=${wd}/result/faprotax/
mkdir -p ${wd2} && cd ${wd2}
# 设置脚本目录
sd=${db}/EasyMicrobiome/script/FAPROTAX_1.2.7
### 1. 软件安装
# 注：软件已经下载至 EasyMicrobiome/script目录，在qiime2环境下运行可满足依赖关系
#(可选)下载软件新版本，以1.2.7版为例， 2023/7/14更新数据库
#wget -c https://pages.uoregon.edu/slouca/LoucaLab/archive/FAPROTAX/SECTION_Download/MODULE_Downloads/CLASS_Latest%20release/UNIT_FAPROTAX_1.2.7/FAPROTAX_1.2.7.zip
#解压
#unzip FAPROTAX_1.2.7.zip
#新建一个python3环境并配置依赖关系，或进入qiime2 python3环境
conda activate qiime2-2023.7
# source /home/silico_biotech/miniconda3/envs/qiime2/bin/activate
#测试是否可运行，弹出帮助即正常工作
python $sd/collapse_table.py
### 2. 制作输入OTU表
#txt转换为biom json格式
biom convert -i ../otutab_rare.txt -o otutab_rare.biom --table-type="OTU table" --to-json
#添加物种注释
biom add-metadata -i otutab_rare.biom --observation-metadata-fp ../taxonomy2.txt \
-o otutab_rare_tax.biom --sc-separated taxonomy \
--observation-header OTUID,taxonomy
#指定输入文件、物种注释、输出文件、注释列名、属性列名
### 3. FAPROTAX功能预测
#python运行collapse_table.py脚本、输入带有物种注释OTU表tax.biom、
#-g指定数据库位置，物种注释列名，输出过程信息，强制覆盖结果，结果文件和细节
#下载faprotax.txt，配合实验设计可进行统计分析
#faprotax_report.txt查看每个类别中具体来源哪些OTUs
python ${sd}/collapse_table.py -i otutab_rare_tax.biom \
-g ${sd}/FAPROTAX.txt \
--collapse_by_metadata 'taxonomy' -v --force \
-o faprotax.txt -r faprotax_report.txt
### 4. 制作OTU对应功能注释有无矩阵
# 对ASV(OTU)注释行，及前一行标题进行筛选
grep 'ASV_' -B 1 faprotax_report.txt | grep -v -P '^--$' > faprotax_report.clean
# faprotax_report_sum.pl脚本将数据整理为表格，位于public/scrit中
perl ${sd}/../faprotax_report_sum.pl -i faprotax_report.clean -o faprotax_report
# 查看功能有无矩阵，-S不换行
less -S faprotax_report.mat
## 3. Bugbase细菌表型预测
### 1. Bugbase命令行分析
cd ${wd}/result
bugbase=${db}/script/BugBase
/bin/rm -rf bugbase/
# 脚本已经优化适合R4.0，biom包更新为biomformat
Rscript ${bugbase}/bin/run.bugbase.r -L ${bugbase} \
-i gg/otutab.txt -m metadata.txt -c Group -o bugbase/
cd ../
### 2. 其它可用分析
# 使用 https://www.bic.ac.cn/ImageGP
# 官网，https://bugbase.cs.umn.edu/ ，有报错，不推荐
# Bugbase细菌表型预测Linux，详见附录3. Bugbase细菌表型预测
# 32、MachineLearning机器学习
# RandomForest包使用的R代码见advanced/RandomForestClassification和RandomForestRegression
## Silme2随机森林/Adaboost使用代码见EasyMicrobiome/script/slime2目录中的slime2.py，详见附录4
# 使用实战(使用QIIME 2的Python3环境，以在Windows中为例)
conda activate qiime2-2023.7
cd /mnt/d/EasyMicrobiome/script/slime2
#使用adaboost计算10000次(16.7s)，推荐千万次
./slime2.py otutab.txt design.txt --normalize --tag ab_e4 ab -n 10000
#使用RandomForest计算10000次(14.5s)，推荐百万次，支持多线程
./slime2.py otutab.txt design.txt --normalize --tag rf_e4 rf -n 10000
# 33、Evolution进化树
cd ${wd}
mkdir -p result/tree
cd ${wd}/result/tree
## 1. 筛选高丰度/指定的特征
#方法1. 按丰度筛选特征，一般选0.001或0.005，且OTU数量在30-150个范围内
#统计特征表中ASV数量，如总计1609个
tail -n+2 ../otutab_rare.txt | wc -l
#按相对丰度0.2%筛选高丰度OTU
# usearch -otutab_trim ../otutab_rare.txt \
# -min_otu_freq 0.002 \
# -output otutab.txt
#统计筛选OTU表特征数量，总计~81个
# tail -n+2 otutab.txt | wc -l
awk 'BEGIN{OFS=FS="\t"}{if(FNR==1) { print "OTUID";} \
else {for(i=3;i<=NF;i++) if($i>0.2 && a[$1]=="") {print $1; a[$1]=1;}}}' \
../otutab_mean.txt > otutab_high.id
#方法2. 按数量筛选
# #按丰度排序，默认由大到小
# usearch -otutab_sortotus ../otutab_rare.txt \
# -output otutab_sort.txt
# #提取高丰度中指定Top数量的OTU ID，如Top100,
# sed '1 s/#OTU ID/OTUID/' otutab_sort.txt \
# | head -n101 > otutab.txt
#修改特征ID列名
# sed -i '1 s/#OTU ID/OTUID/' otutab.txt
#提取ID用于提取序列
# cut -f 1 otutab.txt > otutab_high.id
# 筛选高丰度菌/指定差异菌对应OTU序列
vsearch -fastx_getseqs ../otus.fa -labels otutab_high.id \
-fastaout otus.fa
head -n 2 otus.fa
## 筛选OTU对物种注释
awk 'NR==FNR{a[$1]=$0} NR>FNR{print a[$1]}' ../taxonomy.txt \
otutab_high.id > otutab_high.tax
#获得OTU对应组均值，用于样本热图
#依赖之前otu_mean.R计算过按Group分组的均值
awk 'NR==FNR{a[$1]=$0} NR>FNR{print a[$1]}' ../otutab_mean.txt otutab_high.id \
| sed 's/#OTU ID/OTUID/' > otutab_high.mean
head -n3 otutab_high.mean