6-跟着science学习宏基因组-从宏基因组中提取16S/18S序列分析1-vsearch分析
[toc]
写在前面
从宏基因组中提取核糖体DNA序列，进行扩增子分析。扩增子数据的分析我们已经熟悉的非常熟悉了，只是从宏基因组中得到，这个过程不够熟悉。其次你以为直接提取出来的序列，直接上vseearch就可以了吗？答案是不行的。
我们细细想来就知道了，宏基因组获取的核糖体rna基因的序列可能是不同片段的，这些片段很可能不是一个区域，所以不能进行无参考聚类。下面这个science给出了他的解决方法，这和我去年发表的isme都是一种方法，只是我使用了gg138数据库，他使用 silva数据库，前者全面且准确，但是序列较少，后者数据库巨大，但长度不一。
bbmap提取宏基因组中扩增子序列
# 安装软件
conda install -c bioconda bbmap
conda install -c agbiome bbtools
#
mkdir extract_16S
bbduk.sh in=./trimmomatic/SUBERR793599_forward_paired.fq.gz in2=./trimmomatic/SUBERR793599_reverse_paired.fq.gz outm=./extract_16S/SUBERR793599.bbduk.fa.gz k=31 ref=./database/ribokmers.fa.gz threads=6 2> ./extract_16S/SUBERR793599.bbduk.log
批量运行bbduk.sh
mkdir extract_16S
# ${base}
for i in ./rawdata_0.01/*_1.fq.gz
do
base=$(basename $i _1.fq.gz)
echo $base
bbduk.sh in=./trimmomatic/${base}_forward_paired.fq.gz in2=./trimmomatic/${base}_reverse_paired.fq.gz outm=./extract_16S/${base}.bbduk.fa.gz k=31 ref=./database/ribokmers.fa.gz threads=6 2> ./extract_16S/${base}.bbduk.log
done
重命名16S序列 awk 命令不能使用要修改
# zcat ./extract_16S/SUBERR793599.bbduk.fa.gz | awk '/^>/ {{$0=">SUBERR793599_\" substr($0,2)}}1' | gzip > ./extract_16S/SUBERR793599.bbduk.rename.fa.gz
zcat ./extract_16S/SUBERR793600.bbduk.fa.gz | awk '/^>/ {{$0=">SUBERR793600_" substr($0,2)}}1' | gzip > ./extract_16S/SUBERR793600.bbduk.rename.fa.gz
zcat ./extract_16S/SUBERR793601.bbduk.fa.gz | awk '/^>/ {{$0=">SUBERR793601_" substr($0,2)}}1' | gzip > ./extract_16S/SUBERR793601.bbduk.rename.fa.gz
zcat ./extract_16S/SUBERR793603.bbduk.fa.gz | awk '/^>/ {{$0=">SUBERR793603_" substr($0,2)}}1' | gzip > ./extract_16S/SUBERR793603.bbduk.rename.fa.gz
zcat ./extract_16S/SUBERR793604.bbduk.fa.gz | awk '/^>/ {{$0=">SUBERR793604_" substr($0,2)}}1' | gzip > ./extract_16S/SUBERR793604.bbduk.rename.fa.gz
zcat ./extract_16S/SUBERR793605.bbduk.fa.gz | awk '/^>/ {{$0=">SUBERR793605_" substr($0,2)}}1' | gzip > ./extract_16S/SUBERR793605.bbduk.rename.fa.gz
zcat ./extract_16S/SUBERR793606.bbduk.fa.gz | awk '/^>/ {{$0=">SUBERR793606_" substr($0,2)}}1' | gzip > ./extract_16S/SUBERR793606.bbduk.rename.fa.gz
zcat ./extract_16S/SUBERR793607.bbduk.fa.gz | awk '/^>/ {{$0=">SUBERR793607_" substr($0,2)}}1' | gzip > ./extract_16S/SUBERR793607.bbduk.rename.fa.gz
zcat ./extract_16S/SUBERR793608.bbduk.fa.gz | awk '/^>/ {{$0=">SUBERR793608_" substr($0,2)}}1' | gzip > ./extract_16S/SUBERR793608.bbduk.rename.fa.gzzcat ./extract_16S/*.bbduk.rename.fa.gz | gzip > ./extract_16S/all.fa.gz
vsearch 跑微生物组 注释是有参哦
vsearch --usearch_global ./extract_16S/all.fa.gz -db ~/db/silva/Silva119_release/rep_set/97/Silva_119_rep_set97.fna -id 0.97 --strand both --uc ./extract_16S/mapping/project.uc --threads 6 --log ./extract_16S/mapping/vsearch.log
mkdir ./extract_16S/otutable/
biom from-uc -i ./extract_16S/mapping/project.uc -o ./extract_16S/otutable/project.biom
#--biom格式的转化 hdf5
biom convert -i ./extract_16S/otutable/project.biom -o ./extract_16S/otutable/table.from_biom.txt --to-tsv
# txt转化为biom
biom convert -i ./extract_16S/otutable/table.from_biom.txt -o ./extract_16S/otutable/project_hdf5.biom --table-type="OTU table" --to-hdf5
# txt转化为biom json
biom convert -i ./extract_16S/otutable/table.from_biom.txt -o ./extract_16S/otutable/project_hdf5.biom --table-type="OTU table" --to-json
biom 添加物种注释
下载数据库；
wget -c https://www.arb-silva.de/fileadmin/silva_databases/qiime/Silva_119_consensus_majority_taxonomy.zip
wget -c https://www.arb-silva.de/fileadmin/silva_databases/qiime/Silva_119_release.zip
使用下面数据库成功注释，注意如果用高版本的代表序列做，确实会增加可注释的比例，119版本的是39%，132版本48%，数据库还是丰富了不少的，但是注意添加物种信息也要用同一个版本的，或者用高版本的。
biom add-metadata -i ./extract_16S/otutable/project.biom -o ./extract_16S/taxonomy/project_tax.biom --observation-metadata-fp ~/db/silva/consensus_majority_taxonomy_SILVA119/consensus_taxonomy/97/taxonomy_97_7_levels_consensus.txt --observation-header OTUID,taxonomy --sc-separated taxonomy --float-fields confidenc
附
biom 添加物种注释问题-数据库问题
mkdir ./extract_16S/taxonomy/
# 物种注释结果匹配出错
biom add-metadata -i ./extract_16S/otutable/project_hdf5.biom -o ./extract_16S/taxonomy/project_tax.biom --observation-metadata-fp ~/db/silva/db/taxmap_ncbi_ssu_ref_132.txt --observation-header OTUID,taxonomy --sc-separated taxonomy --float-fields confidenc# 错误信息
raise BiomParseException("First column values are not unique! "
biom.exception.BiomParseException: First column values are not unique! Cannot be ids.
这主要是之前的注释文件和现在的不一样
之前silva做了qiime版本的注释库，但是现在随着qiime的老化，退出来市场，实际上qiime有许多小工具都是很有用的。
我们下载这个老数据库，就可以了。
wget -c https://www.arb-silva.de/fileadmin/silva_databases/qiime/Silva_119_consensus_majority_taxonomy.zip
那么大家就可能有一个问题，我做一个新版本的不就可以了，其实silva网站上面是是有新版本的，只是比较大，这次我们不做下载,具体查看：https://www.arb-silva.de/no_cache/download/archive/qiime/
最新更新的132版本压缩文件有3G大小，如果要体验新版本，请在这里下载
wget -c https://www.arb-silva.de/fileadmin/silva_databases/qiime/Silva_132_release.zip
根际互作生物学研究室 简介
根际互作生物学研究室是沈其荣教授土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军副教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 isme J， Microbiome， PCE，SBB，Horticulture Research等期刊上发表了多篇文章。欢迎关注 微生信生物 公众号对本研究小组进行了解。