扩增子分析——usearch+vsearch+qiime1

参考文章:
1.//www.greatytc.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

运行:

  1. 双端reads的合并
usearch -fastq_mergepairs  rawdata/*R1*.fq -fastqout  usearch_analysis/all_samples_merged.fq -relabel @
  1. 引物剔除和数据质控
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
  1. 序列去重复
 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
  1. 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
  1. 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
  1. 进化树构建
    进化树是基于多序列比对的结果,可展示丰富的信息,我们将在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
  1. 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/*
  1. 按物种分类级别分类汇总
    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、统计软件以及在线分析。

最后编辑于
©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 211,639评论 6 492
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 90,277评论 3 385
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 157,221评论 0 348
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 56,474评论 1 283
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 65,570评论 6 386
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 49,816评论 1 290
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 38,957评论 3 408
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 37,718评论 0 266
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 44,176评论 1 303
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 36,511评论 2 327
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 38,646评论 1 340
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 34,322评论 4 330
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 39,934评论 3 313
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 30,755评论 0 21
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 31,987评论 1 266
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 46,358评论 2 360
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 43,514评论 2 348

推荐阅读更多精彩内容