扩增子分析——usearch+vsearch+qiime1
2020-12-22 本文已影响0人
wanghaihua888
参考文章：
1.https://www.jianshu.com/p/c72bb359f050
2.http://blog.sciencenet.cn/blog-3334560-1071618.html
usearch下载地址：https://drive5.com/software.html
usearch安装：
1. 解压缩
2. chmod +x /apps/users/user01/wanghhh/software/amplicon_example_workflow/usearch
3. export PATH=/apps/users/user01/wanghhh/software/amplicon_example_workflow:$PATH
vsearch安装：
conda install -c bioconda vsearch
#修改名字
ls *.fq |cut -d"_" -f 1 |sort -u |while read id; do
mv ${id}_1R.fq ${id}_R1.fq
mv ${id}_2R.fq ${id}_R2.fq
done
数据质控
fastqc
运行：
- 双端reads的合并
usearch -fastq_mergepairs rawdata/*R1*.fq -fastqout usearch_analysis/all_samples_merged.fq -relabel @
- 引物剔除和数据质控
2.1设置引物文件
head primers.fa
>forward_primer
GTGCCAGCMGCCGCGGTAA
>reverse_primer
GGACTACHVGGGTWTCTAAT
2.2#生成随机抽样5000条序列文件 all_sub_for_primer_check.fq
usearch -fastx_subsample all_samples_merged.fq -sample_size 5000 -fastqout all_sub_for_primer_check.fq
2.3 #在随机抽样文件中搜索检测引物序列，根据生成文件确定引物位置
usearch -search_oligodb all_sub_for_primer_check.fq -db primers.fa -strand both -userout primer_hits.txt -userfields query+qlo+qhi+qstrand
2.4 质控，切除引物
bacterial: vsearch -fastq_filter all_samples_merged.fq --fastq_stripleft 30 --fastq_stripright 25 -fastq_maxee 1 --fastq_maxlen 440 --fastq_minlen 150 --fastaout QCd_merged.fa --fastq_qmax 42
vsearch -fastq_filter all_samples_merged.fq --fastq_stripleft 30 --fastq_stripright 150 -fastq_maxee 1 --fastq_maxlen 440 --fastq_minlen 110 --fastaout QCd_merged.fa --fastq_qmax 42
fungal: vsearch -fastq_filter all_samples_merged.fq --fastq_stripleft 35 --fastq_stripright 150 -fastq_maxee 1 --fastq_maxlen 500 --fastq_minlen 150 --fastaout QCd_merged.fa --fastq_qmax 42
- 序列去重复
vsearch --derep_fulllength QCd_merged.fa -sizeout -relabel Uniq -output unique_seqs.fa
4.聚类OTU或者产生ASVs
4.1方法1
bacterial: usearch -unoise3 unique_seqs.fa -zotus ASVs.fa -minsize 5
fungal: usearch -unoise3 unique_seqs.fa -zotus ASVs.fa -minsize 27
4.2方法2
usearch -cluster_otus unique_seqs.fa \
-otus otus.fa \
-relabel OTU_
5.生成ASVs的统计表
5.1 #将文件内的每个序列的名字Zotu改为ASV开头。
sed -i.tmp 's/Zotu/OTU_/' ASVs.fa
rm ASVs.fa.tmp
5.2 #将每个样品的质控合并后的序列map到 ASVs.fa上。
bacterial：vsearch -usearch_global QCd_merged.fa --db ASVs.fa --id 0.970 --otutabout ASV_counts.txt
fungal: vsearch -usearch_global QCd_merged.fa --db ASVs.fa --id 0.97 --otutabout ASV_counts.txt
用awk命令计算文件中某一列的总和
awk 'BEGIN{sum=0}{sum+=$2}END{print sum}' ASV_counts.txt
- ASVs注释 -物种注释（通过改变参数调节注释效果）
6.1 操作环境及文件安装
操作环境:qiime1：http://qiime.org/1.9.0/install/index.html#
rdp-classifier安装：http://qiime.org/1.9.0/install/install.html#rdp-install
echo "export RDP_JAR_PATH=/apps/users/user01/wanghhh/software/packages/RDP_classifier/qiime-master-rdp_classifier_2.2/rdp_classifier_2.2/rdp_classifier-2.2.jar " >> $HOME/.bashrc
source $HOME/.bashrc
6.2 运行注释
（1）vsearch运行 （生成biom文件，不可读）
vsearch --usearch_global ASVs.fa --db /apps/users/user01/wanghhh/qiime2/usearch_database/rdp_16s_v16.fa --biomout out_tax.txt --id 0.97
（2）qiime1运行（主流）
#bacterial
assign_taxonomy.py \
-i usearch_analysis/ASVs.fa \
-r /apps/users/user01/wanghhh/qiime2/qiime_database/gg_13_8_otus/rep_set/99_otus.fasta \
-t /apps/users/user01/wanghhh/qiime2/qiime_database/gg_13_8_otus/taxonomy/99_otu_taxonomy.txt \
-m rdp \
--confidence=0.5 --min_consensus_fraction=0.51 --similarity=0.8 \
-o usearch_analysis/
--blast_e_value / #使用blast方法时
# fungal:
assign_taxonomy.py \
-i usearch_analysis/ASVs.fa \
-r /apps/users/user01/wanghhh/qiime2/qiime_database/its_12_11_otus/rep_set/99_otus.fasta \
-t /apps/users/user01/wanghhh/qiime2/qiime_database/its_12_11_otus/taxonomy/99_otu_taxonomy.txt \
-m rdp \
--confidence=0.5 --min_consensus_fraction=0.51 --similarity=0.9 \
-o usearch_analysis/
--blast_e_value=0.001\
（3）usearch运行
usearch -sintax ASVs.fa -db /apps/users/user01/wanghhh/qiime2/usearch_database/rdp_16s_v16.fa -tabbedout ASV_tax_raw.txt -strand both -sintax_cutoff 0.4
fungal:
usearch -sintax ASVs.fa -db /apps/users/user01/wanghhh/qiime2/qiime_database/silva_18s_v123.fa -tabbedout ASVs_tax_assignments.txt -strand both -sintax_cutoff 0.6
- OTU表统计、格式转换、添加信息
将OTU表转换为Biom格式，这样便于其它软件对其操作。可添加上面获得的物种信息，这样表格的信息就更丰富了，再转换为文本，便于人类可读，同时使用summarize-table查看OTU表的基本信息。
7.1# 文本OTU表转换为BIOM：方便操作
biom使用：http://biom-format.org/documentation/biom_conversion.html
biom convert \
-i ASV_counts.txt \
-o ASV_counts.biom \
--table-type="OTU table" --to-json
7.2 # 添加物种信息至OTU表最后一列，命名为taxonomy
biom add-metadata -i ASV_counts.biom \
--observation-metadata-fp ASVs_tax_assignments.txt \
-o otu_table_tax.biom \
--sc-separated taxonomy --observation-header OTUID,taxonomy
7.3 # 转换biom为txt格式，带有物种注释：人类可读
biom convert -i otu_table_tax.biom -o otu_table_tax.txt --to-tsv --header-key taxonomy
查看OTU表的基本信息：样品，OUT数量统计
biom summarize-table -i otu_table_tax.biom -o otu_table_tax.sum
Num samples: 27 # 样品数据
Num observations: 975 # OTU数据
Total count: 409647 # 总数据量
Table density (fraction of non-zero values): 0.464 # 非零的单元格
Counts/sample summary:
Min: 2352.0 # 样品数据量最小值
Max: 35955.0 # 样品数据量最大值
Median: 14851.000 # 样品数据量中位数
Mean: 15172.111 # 样品数据量平均数
Std. dev.: 10691.823 # 样品数据量标准变异
Sample Metadata Categories: None provided # 样品分类信息：末提供
Observation Metadata Categories: taxonomy # 观察值分类：物种信息
Counts/sample detail: # 每个样品的数据量
8 OTU表筛选
过滤条件是根据经验、相关文献设计的，如果不清楚，也不要随便过滤，容易引起假阴性。得到的最终结果，还要转换为文本格式，和提取OTU表对应的序列，用于下游分析。
8.1 按样品数据量过滤：选择counts>3000的样品
filter_samples_from_otu_table.py -i result/otu_table_tax.biom -o result/otu_table2.biom -n 3000
查看过滤后结果：
biom summarize-table -i result/otu_table2.biom
8.2 # 按OTU丰度过滤：选择相对丰度均值大于万分之一的OTU
同时还要过滤低丰度的OTU，一般低于万分之一丰度的菌，在功能研究可能还是比较困难的(早期文章454测序数据量少，通常只关注丰度千分之五以上的OTU)。
filter_otus_from_otu_table.py --min_count_fraction 0.0001 -i result/otu_table2.biom -o result/otu_table3.biom
查看过滤后结果：
biom summarize-table -i result/otu_table3.biom
8.3 # 按物种筛选OTU表：去除p__Chloroflexi菌门
有些研究手段在特定有实验中存在偏差，如2012Nature报导V5-V7在植物中扩增会偏好扩增Chloroflexi菌门，建议去除。
filter_taxa_from_otu_table.py -i result/otu_table3.biom -o result/otu_table4.biom -n p__Chloroflexi
查看过滤后结果：
biom summarize-table -i result/otu_table4.biom
8.4# 转换最终biom格式OTU表为文本OTU表格
biom convert -i result/otu_table4.biom -o result/otu_table4.txt --table-type="OTU table" --to-tsv
OTU表格式调整方便R读取
sed -i '/# Const/d;s/#OTU //g;s/ID.//g' result/otu_table4.txt
筛选最终OTU表中对应的OTU序列
filter_fasta.py -f result/rep_seqs.fa -b result/otu_table4.biom -o result/rep_seqs4.fa
- 进化树构建
进化树是基于多序列比对的结果，可展示丰富的信息，我们将在R绘图中详细解读。此处只是建树，用于Alpha, Beta多样性分析的输入文件。
clustalo多序列比对，如果没有请安装Clustal Omega
clustalo -i rep_seqs.fa -o rep_seqs_align.fa --seqtype=DNA --full --force --threads=10
筛选结果中保守序列和保守区
filter_alignment.py -i rep_seqs_align.fa -o temp/ # rep_seqs_align_pfiltered.fa, only very short conserved region saved
基于fasttree建树
make_phylogeny.py -i temp/rep_seqs_align_pfiltered.fasta -o rep_seqs.tree # generate tree by FastTree
10.Alpha多样性计算
Alpha多样性计算前需要对OTU表进行标准化，因为不同测序深度，检测到的物种数量会不同。我们将OTU表重抽样至相同数据量，以公平比较各样品的物种数量。方法如下：
# 查看样品的数据量最小值
biom summarize-table -i otu_table_tax.biom
# 基于最小值进行重抽样标准化
single_rarefaction.py -i otu_table_tax.biom -o temp/otu_table_rare.biom -d 8685
# 计算常用的四种Alpha多样性指数
alpha_diversity.py -i temp/otu_table_rare.biom -o alpha.txt -t rep_seqs.tree -m shannon,chao1,observed_otus,PD_whole_tree
- Beta多样性
Beta多样性是计算各样品间的相同或不同，OTU表也需要标准化。采用重抽样方法丢失的信息太多，不利于统计。此步我们选择CSS标准化方法。
# CSS标准化OTU表
normalize_table.py -i otu_table_tax.biom -o temp/otu_table_css.biom -a CSS
# 转换标准化OTU表为文本，用于后期绘图
biom convert -i temp/otu_table_css.biom -o result/otu_table_css.txt --table-type="OTU table" --to-tsv
# 删除表格多余信息，方便R读取
sed -i '/# Const/d;s/#OTU //g;s/ID.//g' otu_table_tax.txt > otu_table_tax_1.txt
# 计算Beta多样性
beta_diversity.py -i otu_table_tax.biom -o beta/ -t rep_seqs.tree -m bray_curtis,weighted_unifrac,unweighted_unifrac
# Beta多样性距离文件整理，方便R读取
sed -i 's/^\t//g' beta/*
- 按物种分类级别分类汇总
OTU表中最重要的注释信息是物种注释信息。通常的物种注释信息分为7个级别：界、门、纲、目、科、属、种。种是最小的级别，和OTU类似但有不相同。
我们除了可以比较样品和组间OTU水平差异外，还可以研究不同类似级别上的差异，它们是否存在那些共同的变化规律。
按照注释的级别进行分类汇总，无论是Excel还R操作起来，都是很麻烦的过程。这里我们使用QIIME自带 的脚本summarize_taxa.py。
# 结果按门、纲、目、科、属五个级别进行分类汇总，对应结果的L2-L6
summarize_taxa.py -i otu_table_tax.biom -o sum_taxa # summary each level percentage
# 修改一下文本表头，适合R读取的表格格式
sed -i '/# Const/d;s/#OTU ID.//g' sum_taxa/* # format for R read
# 以门为例查看结果
less -S sum_taxa/otu_table4_L2.txt
转入R、统计软件以及在线分析。
