EasyAmplicon (易扩增子)-扩增子高通量序列分析软件流程及脚本-详细使用方法——来自刘永鑫团队的秘籍_easyamplicon能分析上百个样本吗

网上学习资料一大堆，但如果学到的知识不成体系，遇到问题时只是浅尝辄止，不再深入研究，那么很难做到真正的技术提升。

需要这份系统化的资料的朋友，可以点击这里获取！

一个人可以走的很快，但一群人才能走的更远！不论你是正从事IT行业的老鸟或是对IT行业感兴趣的新人，都欢迎加入我们的的圈子（技术交流、学习资源、职场吐槽、大厂内推、面试辅导），让我们一起学习成长！

或者一般下载链接


[https://github.com/YongxinLiu/EasyAmplicon/archive/refs/heads/master.zip]( )


百度网盘（这个真贴心，国内不担心上不了github了：


[百度网盘 请输入提取码]( )


这里还有相关的近500个R包直接下载，省去安装和调试R了，不过不同环境下可能出现小错误，需要有一定的调试技能：[百度网盘 请输入提取码]( )


## 2、开始使用


这里这个主要是扩增子的流程，这里刘大神收集了大家分析常用的软件，缺少软件的先看看这个仓库吧，[GitHub - YongxinLiu/EasyMicrobiome]( )


这里有EasyMicrobiome的介绍 ：[EasyMicrobiome-易扩增子、易宏基因组等分析流程依赖常用软件、脚本文件和数据库注释文件-CSDN博客]( )



习惯用linux了，所以按照linux方法做， 先打开EasyAmplicon下的pipeline.sh文件，建议使用notepad或者sublime txt


运行前可以先将EasyAmplicon文件加入系统环境，这样使用脚本的时候不需要带EasyAmplicon文件夹的绝对路径了。按个人喜好来定吧。


先进入全新目录：

#进入工作目录，工作目录名字可以自己修改
cd /work

1. 分析流程pipeline.sh

2. 样本元信息metadata.txt，保存于result目录

3. 测序数据fastq文件保存于seq目录，通常以.fq.gz结尾，每个样品一对文件

4. 创建临时文件存储目录，分析结束可删除

mkdir -p seq result temp

 先按照大神自己的实验设计描述整理一个全信息的表格，表头设置为一样，这里就不解释字段意思了，会实验设计的基本都了解，样例文件在EasyAmplicon/result下，大家同样整理一个名为metadata\_raw.txt的文件放入新的work文件夹


![](https://img-blog.csdnimg.cn/direct/61a8c620c23e42879f323a49c291a2f7.png)


在work文件夹下，将metadata\_raw.txt进一步去除结尾符号，主要是windows系统和linux系统差异，并命名为metadata.txt

# windows用户结尾有^M，运行sed命令去除，再用cat -A检查
sed 's/\r//' result/metadata_raw.txt > result/metadata.txt
cat -A result/metadata.txt | head -n3

如果直接在linux下建立的文件大家不需要这个步骤了，直接按实验设计整理一个metadata.txt文件即可，如果是windows下写的设计文件建议在ftp上传的时候选择使用二进制方式上传就不用这个步骤处理了。


进入seq文件夹

cd seq

下载样品序列，放入seq目录下


已有样品数据在EasyAmplicon/seq下，这里是每个样品一对序列文件，如果样品是混合样可以先拆样后按相同的方式将序列文件命名后放入下面目录：


![](https://img-blog.csdnimg.cn/direct/d950c3bd75a946d89456f9b66fd4018c.png)


统计测序数据，依赖seqtk程序，没有的自己搜索找一个吧：

#seqkit统计测序数据
seqkit stat seq/*.fq.gz > result/seqkit.txt

#查看统计信息
head result/seqkit.txt

![](https://img-blog.csdnimg.cn/direct/2a860dfe3bdc49c89ac0182277052e5a.png)


单个文件的查询方法在pipeline里面已有说明，大家按照里面的命令执行即可查看，不再细说 


### 序列合并并按样品重命名，仍然为单个样品

#测试.以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

可以使用for循环或者rush进行批量处理

#方法1.for循环顺序处理
time for i in tail -n+2 result/metadata.txt|cut -f1;do
vsearch --fastq_mergepairs seq/$i_1.fq.gz--reverseseq/${i}_2.fq.gz
–fastqout temp/${i}.merged.fq --relabel ${i}.
done &

一部分电脑 rush 不支持，运行时调度失败，请使用 for 循环部分

#方法2.rush并行处理，任务数jobs(j),2可加速1倍4s；建议设置2-4
time tail -n+2 result/metadata.txt | cut -f 1 |
rush -j 2 “vsearch --fastq_mergepairs seq/{}_1.fq.gz --reverse seq/{}_2.fq.gz
–fastqout temp/{}.merged.fq --relabel {}.”

序列合并处理完成后，所有的序列名称都加上了label，所以为提高效率，后面将合并所有样品的序列在一起


### 合并所有样品序列


从这里开始就与混合样品测序的样品在merge后进行relabel之后的状态一致了。

#合并所有样品至同一文件
cat temp/*.merged.fq > temp/all.fq

#查看文件大小223M，软件不同版本结果略有差异
ls -lsh temp/all.fq

查看序列名，“.”之前是否为样本名，样本名绝不允许有点 (“.”)

样本名有点 (.) 的一个显著特征是生成的特征表会很大，特征表里面列很多，导致后面分析出现内存不足。

后面分析获得特征表后要看一眼有没有问题，遇到内存不足问题，也要回头来排查。

head -n 6 temp/all.fq|cut -c1-60

### 切除引物和质控

左端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

### 序列去冗余 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和频率

### 聚类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

查看otus.fa，序列名称都更改了。


![](https://img-blog.csdnimg.cn/direct/f51f20e0c7d7459c8749dc8fd690d726.png)


### 基于参考去嵌合 Reference-based chimera detect


这里如果去除嵌合体的话就需要有rdp的数据库了，需提前下载好

不推荐，容易引起假阴性，因为参考数据库无丰度信息

而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不要执行此命令，linux下应该也不需要

sed -i ‘s/\r//g’ result/raw/otus.fa

方法2. 不去嵌合，直接cp到指定工作目录

cp -f temp/otus.fa result/raw/otus.fa

### 特征表构建和筛选 Feature table create and filter


# OTU和ASV统称为特征(Feature)，它们的区别是：  
 # OTU通常按97%聚类后挑选最高丰度或中心的代表性序列；  
 # ASV是基于序列进行去噪(排除或校正错误序列，并挑选丰度较高的可信序列)作为代表性序列


#### 生成特征表

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

#### 物种注释，且/或去除质体和非细菌 Remove plastid and non-Bacteria


 先统计otu表的特征数，包括每个样品的otu总数


其他pipeline的筛选可根据需求和喜好来吧，不确定或感觉没必要的掠过即可，直接统计otu表信息

物种注释-去除质体和非细菌/古菌并统计比例(可选)

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

sed -i ‘s/\r//’ result/raw/otus.sintax
########################
#物种注释结果：
ASV_7 d:Bacteria(1.00),p:Bacteroidetes(1.00),c:Flavobacteriia(1.00),o:Flavobacteriales(1.00),f:Flavobacteriaceae(1.00),g:Flavobacterium(1.00) + d:Bacteria,p:Bacteroidetes,c:Flavobacteriia,o:Flavobacteriales,f:Flavobacteriaceae,g:Flavobacterium
ASV_5 d:Bacteria(1.00),p:Cyanobacteria/Chloroplast(1.00),c:Chloroplast(1.00),f:Chloroplast(1.00),g:Streptophyta(1.00) + d:Bacteria,p:Cyanobacteria/Chloroplast,c:Chloroplast,f:Chloroplast,g:Streptophyta
ASV_4 d:Bacteria(1.00),p:Proteobacteria(1.00),c:Betaproteobacteria(1.00),o:Burkholderiales(1.00),f:Comamonadaceae(1.00),g:Rhizobacter(0.98) + d:Bacteria,p:Proteobacteria,c:Betaproteobacteria,o:Burkholderiales,f:Comamonadaceae,g:Rhizobacter
ASV_3 d:Bacteria(1.00),p:Proteobacteria(1.00),c:Betaproteobacteria(1.00),o:Burkholderiales(1.00),f:Comamonadaceae(1.00),g:Rhizobacter(0.93) + d:Bacteria,p:Proteobacteria,c:Betaproteobacteria,o:Burkholderiales,f:Comamonadaceae,g:Rhizobacter
ASV_2 d:Bacteria(1.00),p:Proteobacteria(1.00),c:Betaproteobacteria(1.00),o:Burkholderiales(1.00),f:Comamonadaceae(1.00),g:Pelomonas(1.00) + d:Bacteria,p:Proteobacteria,c:Betaproteobacteria,o:Burkholderiales,f:Comamonadaceae,g:Pelomonas

#OTU表简单统计 Summary OTUs table
usearch -otutab_stats result/otutab.txt
-output result/otutab.stat

cat result/otutab.stat

207572  Reads (207.6k)
    18  Samples
  1363  OTUs

 24534  Counts
  4507  Count  =0  (18.4%)
  4229  Count  =1  (17.2%)
  3999  Count >=10 (16.3%)

   413  OTUs found in all samples (30.3%)
   576  OTUs found in 90% of samples (42.3%)
  1276  OTUs found in 50% of samples (93.6%)

Sample sizes: min 10425, lo 10953, med 11549, mean 11531.8, hi 11924, max 13209

####  等量抽样标准化 Normlize by subsample


这里主要是依赖R包了，使用R运行，可以脚本运行，或者直接将otutab\_rare.R文件打开在R终端或studio的交互窗口直接运行

#使用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

#### 计算α多样性 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程序来进行alpha多样性的计算

#### 计算稀释丰富度 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
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

#### 筛选高丰度菌 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]=KaTeX parse error: Expected 'EOF', got '}' at position 24: … "OTU","Group";}̲ \     else {fo…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/ 中绘图并显示各组共有和特有维恩或网络图

也可在 http://www.ehbio.com/ImageGP 绘制Venn、upSetView和Sanky

#### β多样性 Beta diversity

#结果有多个文件，需要目录
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/

#### 物种注释分类汇总

#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
  usearch -sintax_summary result/otus.sintax
  -otutabin result/otutab_rare.txt -rank KaTeX parse error: Expected group after '_' at position 31: … 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

#### 有参定量特征表

比对Greengenes97% OTUs比对，用于PICRUSt/Bugbase功能预测

mkdir -p result/gg/

usearch比对更快，但文件超限报错选附录14 vsearch比对

默认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

#统计
usearch -otutab_stats result/gg/otutab.txt -output result/gg/otutab.stat
cat result/gg/otutab.stat

#### 空间清理及数据提交


这个步骤应该放在最后文章撰写完成后再处理，因为有些中间文件可能需要调整参数重新处理，有些大文件不用再重复cp了，另外，提交数据也是在文章接收后才会用到。

#删除中间大文件
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

#### 统计分析及出图


##### Alpha多样性

#Alpha多样性箱线图，对于不同样品数的作图，建议单个运行后再逐个调整参数以使出图布局最佳

查看帮助

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 多样性维恩图

三组比较:-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

EVenn在线绘制维恩图 https://www.ehbio.com/test/venn

##### 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

##### 主坐标分析PCoA

输入文件，选择分组，输出文件，图片尺寸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/

#####  限制性主坐标分析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

#####  物种组成Taxonomy

3.1 堆叠柱状图Stackplot

以门§水平为例，结果包括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级

for i in p c o f g; do
Rscript KaTeX parse error: Expected group after '_' at position 55: … 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 KaTeX parse error: Expected group after '_' at position 54: … result/tax/sum_̲{i}.txt --design result/metadata.txt
  --group Group --legend 5

结果位于当前目录circlize.pdf(随机颜色)，circlize_legend.pdf(指定颜色+图例)

移动并改名与分类级一致

mv circlize.pdf result/tax/sum_KaTeX parse error: Expected group after '_' at position 55: … 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

##### 差异比较

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 火山图

输入compare.R的结果，输出火山图带数据标签，可指定图片大小

Rscript $db/script/compar{e}_{v}olcano.R--inputresult/compare/${compare}.txt
  --output result/compare/${compare}.volcano.pdf
  --width 89 --height 59

1.3 热图

输入compare.R的结果，筛选列数，指定元数据和分组、物种注释，图大小英寸和字号

bash $db/script/compar{e}_{h}eatmap.sh-iresult/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/compar{e}_{m}anhattan.sh-iresult/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/compar{e}_{m}anhattan.sh-iresult/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/compar{e}_{m}anhattan.sh-iresult/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降序

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.三元图的应用与绘图实战

##### 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--threshold0.1--method"t.test"--pvalue0.05--fdr"none"--width189--height159--outputresult/stamp/${compare}

可选Rmd文档见result/CompareStamp.Rmd

##### 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?page=b%27MzY%3D%27
#方法2. LEfSe本地分析(限Linux系统、选学)，参考代码见附录
#方法3. LEfSe官网在线使用

####  QIIME 2分析流程


代码详见 qiime2/pipeline\_qiime2.sh


#####  功能预测


## 1. PICRUSt功能预测


### PICRUSt 1.0

方法1. 使用 http://www.ehbio.com/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
paste <(cut -f $numresult/picrust/${l}.txt) <(cut -f 1-$[num-1]result/picrust/${l}.txt)
  > result/picrust/$l.spfcut-f2-result/picrust/${l}.spf > result/picrust/${l}.mat.txt
awk 'BEGIN{FS=OFS=“\t”} {print $2,KaTeX parse error: Expected 'EOF', got '}' at position 2: 1}̲' result/picrus…{l}.spf | sed ‘s/;/\t/’ | sed ‘1 s/ID/Pathway\tCategory/’
  > result/picrust/${l}.anno.txt

差异比较

compare=“KO-WT”
Rscript $db/script/compare.R--inputresult/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/compar{e}_{h}ierarch{y}_{f}acet.R--inputresult/picrust/${compare}.txt
  --data MeanA
  --annotation result/picrust/$l.anno.txt--outputresult/picrust/${compare}.MeanA.bar.pdf

绘制两组显著差异柱状图，按高分类级分面

Rscript $db/script/compar{e}_{h}ierarch{y}_{f}acet2.R--inputresult/picrust/${compare}.txt
  --pvalue 0.05 --fdr 0.1
  --annotation result/picrust/$l.anno.txt--outputresult/picrust/${compare}.bar.pdf

 ### PICRUSt 2.0

软件安装见附录6. PICRUSt环境导出和导入

(可选)PICRUSt2(Linux/Windows下Linux子系统，要求>16GB内存)

安装参考附录5的方式直接下载安装包并解压即可使用

Linux中加载conda环境

conda activate picrust2

进入工作目录，服务器要修改工作目录

wd=/mnt/d/amplicon/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

##### 元素循环FAPROTAX

方法1. 在线分析，推荐使用 http://www.bic.ac.cn/ImageGP/index.php/Home/Index/FAPROTAX.html 在线分析

方法2. Linux下分析、如QIIME 2环境下

设置工作目录

wd=/mnt/d/amplicon/result/faprotax/
mkdir -p ${wd} && cd ${wd}

设置脚本目录

sd=/mnt/d/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

##### Bugbase细菌表型预测

1. Bugbase命令行分析

cd $wd/resultbugbase=${db}/script/BugBase
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/

2. 其它可用分析

使用 http://www.bic.ac.cn/ImageGP/index.php/Home/Index/BugBase.html

官网，https://bugbase.cs.umn.edu/ ，有报错，不推荐

Bugbase细菌表型预测Linux，详见附录3. Bugbase细菌表型预测

##### 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

####  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

#方法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序列

usearch -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

#合并物种注释和丰度为注释文件
cut -f 2- otutab_high.mean > temp
paste otutab_high.tax temp > annotation.txt
head -n 3 annotation.txt

为了做好运维面试路上的助攻手，特整理了上百道 【运维技术栈面试题集锦】 ，让你面试不慌心不跳，高薪offer怀里抱！

这次整理的面试题，小到shell、MySQL，大到K8s等云原生技术栈，不仅适合运维新人入行面试需要，还适用于想提升进阶跳槽加薪的运维朋友。

本份面试集锦涵盖了

• 174 道运维工程师面试题
• 128道k8s面试题
• 108道shell脚本面试题
• 200道Linux面试题
• 51道docker面试题
• 35道Jenkis面试题
• 78道MongoDB面试题
• 17道ansible面试题
• 60道dubbo面试题
• 53道kafka面试
• 18道mysql面试题
• 40道nginx面试题
• 77道redis面试题
• 28道zookeeper

总计 1000+ 道面试题， 内容 又全含金量又高

• 174道运维工程师面试题

1、什么是运维?

2、在工作中，运维人员经常需要跟运营人员打交道，请问运营人员是做什么工作的?

3、现在给你三百台服务器，你怎么对他们进行管理?

4、简述raid0 raid1raid5二种工作模式的工作原理及特点

5、LVS、Nginx、HAproxy有什么区别?工作中你怎么选择?

6、Squid、Varinsh和Nginx有什么区别，工作中你怎么选择?

7、Tomcat和Resin有什么区别，工作中你怎么选择?

8、什么是中间件?什么是jdk?

9、讲述一下Tomcat8005、8009、8080三个端口的含义？

10、什么叫CDN?

11、什么叫网站灰度发布?

12、简述DNS进行域名解析的过程?

13、RabbitMQ是什么东西?

14、讲一下Keepalived的工作原理?

15、讲述一下LVS三种模式的工作过程?

16、mysql的innodb如何定位锁问题，mysql如何减少主从复制延迟?

17、如何重置mysql root密码?

网上学习资料一大堆，但如果学到的知识不成体系，遇到问题时只是浅尝辄止，不再深入研究，那么很难做到真正的技术提升。

需要这份系统化的资料的朋友，可以点击这里获取！

一个人可以走的很快，但一群人才能走的更远！不论你是正从事IT行业的老鸟或是对IT行业感兴趣的新人，都欢迎加入我们的的圈子（技术交流、学习资源、职场吐槽、大厂内推、面试辅导），让我们一起学习成长！

锦涵盖了

• 174 道运维工程师面试题
• 128道k8s面试题
• 108道shell脚本面试题
• 200道Linux面试题
• 51道docker面试题
• 35道Jenkis面试题
• 78道MongoDB面试题
• 17道ansible面试题
• 60道dubbo面试题
• 53道kafka面试
• 18道mysql面试题
• 40道nginx面试题
• 77道redis面试题
• 28道zookeeper

总计 1000+ 道面试题， 内容 又全含金量又高

• 174道运维工程师面试题

1、什么是运维?

2、在工作中，运维人员经常需要跟运营人员打交道，请问运营人员是做什么工作的?

3、现在给你三百台服务器，你怎么对他们进行管理?

4、简述raid0 raid1raid5二种工作模式的工作原理及特点

5、LVS、Nginx、HAproxy有什么区别?工作中你怎么选择?

6、Squid、Varinsh和Nginx有什么区别，工作中你怎么选择?

7、Tomcat和Resin有什么区别，工作中你怎么选择?

8、什么是中间件?什么是jdk?

9、讲述一下Tomcat8005、8009、8080三个端口的含义？

10、什么叫CDN?

11、什么叫网站灰度发布?

12、简述DNS进行域名解析的过程?

13、RabbitMQ是什么东西?

14、讲一下Keepalived的工作原理?

15、讲述一下LVS三种模式的工作过程?

16、mysql的innodb如何定位锁问题，mysql如何减少主从复制延迟?

17、如何重置mysql root密码?

网上学习资料一大堆，但如果学到的知识不成体系，遇到问题时只是浅尝辄止，不再深入研究，那么很难做到真正的技术提升。

需要这份系统化的资料的朋友，可以点击这里获取！

一个人可以走的很快，但一群人才能走的更远！不论你是正从事IT行业的老鸟或是对IT行业感兴趣的新人，都欢迎加入我们的的圈子（技术交流、学习资源、职场吐槽、大厂内推、面试辅导），让我们一起学习成长！