-
Notifications
You must be signed in to change notification settings - Fork 116
/
2StatPlot.sh
254 lines (212 loc) · 9.74 KB
/
2StatPlot.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
[TOC]
# 易宏基因组流程EasyMetagenomePipeline
# 版本: 1.21, 2024/5/17
# 测试环境为Windows 10+ / MacOS 10+
# 设置结果目录(通常为项目中的result，此处为12个样本结果result12)
# 可选癌症cancer或年龄age各18个样本数据示例
# mac/linux 需要修改路径格式如下， ~ 代表家目录
# wd=~/meta/result12
# sd=~/EasyMicrobiome/script
wd=/d/EasyMetagenome/result12
# 设置脚本所在目录(Script Directory)，系统为win/mac/linux
sd=/d/EasyMicrobiome/script
PATH=$PATH:$sd/../win:$sd
# 进入结果目录
cd $wd
## 物种Metaphlan4
### 热图
# 显示脚本帮助 help
Rscript ${sd}/metaphlan_hclust_heatmap.R -h
# 按指定分类汇总、排序并取Top25种绘制热图
# -i输入metaphlan4结果转换的spf文件；
# -t指定分类级别，可选Kingdom/Phylum/Class/Order/Family/Genus/Species/Strain(界门纲目科属种株)，推荐门，目，属
# -n 输出物种数量，默认为25，最大值为该类型的数量
# -w、-e指定图片的宽和高，单位为毫米(mm)
# -o输出图pdf、表txt前缀，默认Heatmap+(-t)+(-n)
# 科水平Top25热图
csvtk -t stat metaphlan4/taxonomy.spf
Rscript $sd/metaphlan_hclust_heatmap.R \
-i metaphlan4/taxonomy.spf \
-t Family -n 25 \
-w 183 -e 118 \
-o metaphlan4/HeatmapFamily
# 属水平的Top30
Rscript $sd/metaphlan_hclust_heatmap.R \
-i metaphlan4/taxonomy.spf \
-t Genus -n 30 \
-w 183 -e 118 \
-o metaphlan4/HeatmapGenus
### 维恩图
### 筛选每个样本>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];}}' \
metaphlan4/taxonomy.tsv > metaphlan4/taxonomy_high.tsv
wc -l metaphlan4/taxonomy_high.tsv
# 在线绘制，支持实时查看元素交集 http://www.ehbio.com/test/venn/#/
# 引文：Mei Yang, Tong Chen, Yong-Xin Liu, Luqi Huang. 2024. Visualizing set relationships:
# EVenn's comprehensive approach to Venn diagrams. iMeta 3: e184. https://doi.org/10.1002/imt2.184
# 本地5组比较:-f输入文件,-a/b/c/d/g分组名,-w/u为宽高英寸,-p输出文件名后缀
bash ${sd}/sp_vennDiagram.sh -f metaphlan4/taxonomy_high.tsv \
-a C1 -b C2 -c N1 -d N2 -g N3 \
-w 4 -u 4 \
-p C1_C2_N1_N2_N3
# 求均值再两组比较
sed -i '/^#/d' metaphlan4/taxonomy.tsv
Rscript ${sd}/otu_mean.R --input metaphlan4/taxonomy.tsv \
--metadata metadata.txt \
--group Group --thre 0 \
--scale F --all TRUE --type mean \
--output metaphlan4/group_mean.txt
### 筛选每个组>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];}}' \
metaphlan4/group_mean.txt > metaphlan4/group_high.tsv
bash ${sd}/sp_vennDiagram.sh -f metaphlan4/group_high.tsv \
-a Cancer -b Normal -c All \
-w 4 -u 4 \
-p Cancer_Normal_All
## 功能HUMAnN3
### 分组聚类热图
bash $sd/sp_pheatmap.sh
cut -f 1-2 metadata.txt > group.txt
# -f输入文件，-H水平聚类，u/v图片宽/高，-P添加行注释文件，-Q添加列注释
bash $sd/sp_pheatmap.sh \
-f humann3/pathabundance_relab_unstratified.tsv \
-H 'TRUE' -u 20 -v 50 \
-Q group.txt
# 结果为 输入文件名+pheamap.r/pdf，代码和图片
# 水平标准化，-d row，可选column
bash $sd/sp_pheatmap.sh \
-f humann3/pathabundance_relab_unstratified.tsv \
-H 'TRUE' -u 20 -v 50 \
-Q group.txt -d row
## 物种kraken2
### Alpha多样性
# 提取种级别注释并抽平至最小测序量，计算6种alpha多样性指数
# 查看帮助
Rscript $sd/kraken2alpha.R -h
# -d指定最小样本量，默认0为最小值，抽平文件tax_norm.txt，alpha多样性tax_alpha.txt
Rscript $sd/kraken2alpha.R \
--input kraken2/tax_count.mpa \
--depth 0 \
--species kraken2/tax_count.txt \
--normalize kraken2/tax_count.norm \
--output kraken2/tax_count.alpha
# 绘制Alpha多样性指数，结果为输入文件+类型richness/chao1/ACE/shannon/simpson/invsimpson
# Rscript $sd/alpha_boxplot.R -h # 查看参数
Rscript $sd/alpha_boxplot.R \
-i kraken2/tax_count.alpha \
-a shannon \
-d metadata.txt \
-n Group \
-o kraken2/ \
-w 89 -e 59
# 批量计算6种指数的箱线图+统计
for i in richness chao1 ACE shannon simpson invsimpson;do
Rscript $sd/alpha_boxplot.R -i kraken2/tax_count.alpha -a ${i} \
-d metadata.txt -n Group -w 89 -e 59 \
-o kraken2/
done
### 热图
# 转换为metaphalan2 spf格式，分隔符为下划线“_”
awk 'BEGIN{OFS=FS="\t"}{delete a; a["d"]="unclassified";a["p"]="unclassified";a["c"]="unclassified";a["o"]="unclassified";a["f"]="unclassified";a["g"]="unclassified";a["s"]="unclassified"; \
split($1,x,"|");for(i in x){split(x[i],b,"_");a[b[1]]=b[2];} \
print a["d"],a["p"],a["c"],a["o"],a["f"],a["g"],a["s"],$0;}' \
kraken2/tax_count.txt > temp.txt
cut -f 1-7,9- temp.txt > kraken2/tax_count.spf
sed -i '1 s/unclassified\tunclassified\tunclassified\tunclassified\tunclassified\tunclassified\tunclassified/Domain\tPhylum\tClass\tOrder\tFamily\tGenus\tSpecies/' \
kraken2/tax_count.spf
# 绘制热图，可选域、门、纲、目、科、属、种(Domain Phylum Class Order Family Genus Species)
tax=Genus
Rscript $sd/metaphlan_hclust_heatmap.R \
-i kraken2/tax_count.spf \
-t ${tax} \
-n 25 \
-w 118 -e 118 \
-o kraken2/heatmap_${tax}
### 箱线图
# 绘制属水平Top30箱线图
Rscript $sd/metaphlan_boxplot.R \
-i kraken2/tax_count.spf \
-t Genus \
-n 30 \
-o kraken2/boxplot_Genus
# 绘制门水平Top10箱线图
Rscript $sd/metaphlan_boxplot.R \
-i kraken2/tax_count.spf \
-t Phylum \
-n 10 -w 6 -e 4 \
-o kraken2/boxplot_Phylum
## 物种kraken2-braken2
### Alpha多样性
# 提取种级别注释并抽平至最小测序量，计算6种alpha多样性指数
# 查看帮助
Rscript $sd/otutab_rare.R -h
# -d指定最小样本量，默认0为最小值，抽平文件tax_norm.txt，alpha多样性tax_alpha.
tax=S
Rscript $sd/otutab_rare.R \
--input kraken2/bracken.${tax}.txt \
--depth 0 --seed 1 \
--normalize kraken2/bracken.${tax}.norm \
--output kraken2/bracken.${tax}.alpha
# 绘制Alpha多样性指数，结果为输入文件+类型richness/chao1/ACE/shannon/simpson/invsimpson
# Rscript $sd/alpha_boxplot.R -h # 查看参数
mkdir -p kraken2/${tax}
Rscript $sd/alpha_boxplot.R \
-i kraken2/bracken.${tax}.alpha \
-a shannon \
-d metadata.txt \
-n Group \
-o kraken2/${tax} \
-w 89 -e 59
# 批量计算6种指数的箱线图+统计
for i in richness chao1 ACE shannon simpson invsimpson;do
Rscript $sd/alpha_boxplot.R -i kraken2/bracken.${tax}.alpha -a ${i} \
-d metadata.txt -n Group -w 89 -e 59 \
-o kraken2/${tax}
done
### Beta多样性
# Beta多样性距离矩阵计算
# Mac 用不了面对的 usearch，使用在线平台 https://www.bic.ac.cn/BIC
mkdir -p kraken2/beta/
$sd/../win/usearch -beta_div kraken2/bracken.${tax}.norm \
-filename_prefix kraken2/beta/
# PCoA分析输入文件，选择分组，输出文件，图片尺寸mm，统计见beta_pcoa_stat.txt
# 可选距离有 bray_curtis, euclidean, jaccard, manhattan
dis=bray_curtis
Rscript $sd/beta_pcoa.R \
--input kraken2/beta/${dis}.txt \
--design metadata.txt \
--group Group \
--width 89 --height 59 \
--output kraken2/pcoa.${dis}.pdf
### 堆叠柱状图
# 以门(P)/种(S)水平为例，结果包括output.sample/group.pdf两个文件
tax=P
Rscript ${sd}/tax_stackplot.R \
--input kraken2/bracken.${tax}.txt --design metadata.txt \
--group Group --output kraken2/bracken.${tax}.stackplot \
--legend 10 --width 89 --height 59
## 基因组进化树注释table2itol
cd ../binning/result/itol/
## 方案1. 分类彩带、数值热图、种标签
# -a 找不到输入列将终止运行（默认不执行）-c 将整数列转换为factor或具有小数点的数字，-t 偏离提示标签时转换ID列，-w 颜色带，区域宽度等， -D输出目录，-i OTU列名，-l 种标签替换ID
Rscript ${sd}/table2itol.R -a -c double -D plan1 -i ID -l Species -t %s -w 0.5 annotation.txt
# 生成注释文件中每列为单独一个文件
## 方案2. 数值柱形图，树门背景色，属标签
Rscript ${sd}/table2itol.R -a -d -c none -D plan2 -b Phylum -i ID -l Genus -t %s -w 0.5 annotation.txt
## 方案3.分类彩带、整数为柱、小数为热图
Rscript ${sd}/table2itol.R -c keep -D plan3 -i ID -t %s annotation.txt
## 方案4. 将整数转化成因子生成注释文件
Rscript ${sd}/table2itol.R -a -c factor -D plan4 -i ID -l Genus -t %s -w 0 annotation.txt
# 附录
## windows换行符处理
# 查看行尾是否有^M
cat -A $sd/metaphlan_hclust_heatmap.R | head
# 转换Windows为Linux换行符
dos2unix $sd/metaphlan_hclust_heatmap.R
## 表格合并
# 按末列注释，必须名叫KO
Rscript $sd/mat_gene2ko.R \
-i temp/dbcan2/gene_fam.count \
-o dbcan2/fam_merge.count