5-跟着science学习宏基因组-kraken物种注释

[toc]

写在前面

kraken基于mini数据库。并且这个序列也比较少，所以，很快就能完成

继续处理

胶水操作：提取序列名称zcat ./trimmomatic/SUBERR793599_forward_paired.fq.gz | awk '(NR%4 == 1){{print $0}}' | cut -d' ' -f 1 | cut -c 2- > ./trimmomatic/SUBERR793599_1.names.txt

# ${base}

for i in ./rawdata_0.01/*_1.fq.gz

do

base=$(basename $i _1.fq.gz)

echo $base

zcat ./trimmomatic/${base}_forward_paired.fq.gz | awk '(NR%4 == 1){{print $0}}' | cut -d' ' -f 1 | cut -c 2- > ./trimmomatic/${base}_1.names.txt

done

胶水操作：合并双端数据

合并测序文件：这部分数据存在于./unpack/文件夹，是未经过质控和barcode去除的最原始的数据。# 用于丰度统计 加州大学戴维斯分校的成果

pip install khmer

interleave-reads.py -o ./unpack/SUBERR793599.fastq ./unpack/SUBERR793599_1.fastq.gz ./unpack/SUBERR793599_2.fastq.gz

# ${base}

for i in ./rawdata_0.01/*_1.fq.gz

do

base=$(basename $i _1.fq.gz)

echo $base

interleave-reads.py -o ./unpack/${base}.fastq ./unpack/${base}_1.fastq.gz ./unpack/${base}_2.fastq.gz

done

胶水操作：合并全部样本的双端序列cat ./unpack/SUBERR793599.fastq > ./interleaved_merged/all.fastq

# `tail -n+2 ./map.txt |cut -f 1|sed 's/^//unpack\//;s/$/_.fastq/'| tr '\n' ','|sed 's/,$//'`

这个到这里就先停下了，又开始了另外一条注释线。

kraken注释物种信息

kraken物种注释

kraken也是基于read的注释，这里呢，我们使用triming文件夹中的数据进行注释

安装krakenconda install kraken -c bioconda

kraken 注释流程mkdir kraken

# sudo apt install pigz

pigz ./trimming/*_1.fq

pigz ./trimming/*_2.fq

kraken --db ~/db/kraken/minikraken_20171019_8GB --threads 6 --fastq-input ./trimming/SUBERR793599_1.fq.gz ./trimming/SUBERR793599_2.fq.gz --paired --check-names --gzip-compressed --output ./kraken/SUBERR793599.kraken 2> ./kraken/SUBERR793599.kraken.log

kraken-translate --db ~/db/kraken/minikraken_20171019_8GB --mpa-format ./kraken/SUBERR793599.kraken > ./kraken/SUBERR793599.kraken.taxonomy

kraken-mpa-report --db ~/db/kraken/minikraken_20171019_8GB ./kraken/SUBERR793599.kraken > ./kraken/SUBERR793599.kraken.report

批量运行# ${base}

for i in ./rawdata_0.01/*_1.fq.gz

do

base=$(basename $i _1.fq.gz)

echo $base

kraken --db ~/db/kraken/minikraken_20171019_8GB --threads 6 --fastq-input ./trimming/${base}_1.fq.gz ./trimming/${base}_2.fq.gz --paired --check-names --gzip-compressed --output ./kraken/${base}.kraken 2> ./kraken/${base}.kraken.log

kraken-translate --db ~/db/kraken/minikraken_20171019_8GB --mpa-format ./kraken/${base}.kraken > ./kraken/${base}.kraken.taxonomy

kraken-mpa-report --db ~/db/kraken/minikraken_20171019_8GB ./kraken/${base}.kraken > ./kraken/${base}.kraken.report

done

后续转化物种注释文件：kraken2phyloseq 这部分数据清晰没什么卵用不过咱们也不需要，直接在R进行然后可视化就可以了。# echo -e \"#OTU ID\tSUBERR793599\tConsensus Lineage\" > ./kraken/SUBERR793599.kraken.qiime.taxonomy

# awk 'BEGIN {{OFS=\"\t\"}} {{ print $1,1,$2 }}' ./kraken/SUBERR793599.kraken.taxonomy | sed 's/|/;/g' >> ./kraken/SUBERR793599.kraken.qiime.taxonomy

作者写了这几行，证明大工作走完了，开始作图了== #TODO: ==

== # Make plots for raw reads tax assignment with kraken/diamond ==

== # Convert to a dummy otutable to load with phyloseq with a fake count column ==

== # Add qiime header (todo) ==

# awk 'BEGIN {OFS="\t"} { print $1,1,$2 }' 1511KMI-0007/kraken/MPR1-1n-EST.kraken.taxonomy | sed 's/|/;/g' > test2.qiime

# make plots with gpplot

到这儿，咱们解决了整个代码的500行内容了。

配置数据库

这将下载NCBI分类信息，以及RefSeq中细菌，古细菌和病毒域的完整基因组。下载所有这些数据之后，构建过程开始；这是最耗时的步骤。如果您具有多个处理核心，则可以使用多个线程来运行此过程，例如：kraken-build --standard --threads 24 --db ~/db/kraken

请注意，如果任何步骤（包括初始下载）失败，则构建过程将中止。但是，kraken-build将在整个安装过程中产生检查点，并且如果您尝试在部分构建的数据库上再次运行同一命令，则会在最后一个未完成的步骤中重新启动构建。

构建数据库之后，要除去任何不必要的文件（包括不再需要的库文件），请运行以下命令：kraken-build --db ~/db/kraken --clean

数据库完成后应该有这几个文件database.kdb：包含k -mer到分类单元的映射

database.idx：在database.kdb中包含最小化程序偏移量位置

taxonomy/nodes.dmp：分类树结构+等级

taxonomy/names.dmp：分类名称

如果要构建自定义数据库，请参考：http://ccb.jhu.edu/software/kraken/MANUAL.html#standard-kraken-database

这里考虑到大家的内存

您可以获取现有的Kraken数据库并从中创建较小的MiniKraken数据库。使用此选项将删除除指定数量的k -mer / taxon对，以创建一个新的较小的数据库。例如：kraken-build --shrink 10000 --db ~/db/minikraken --new-db minikraken

根际互作生物学研究室 简介

根际互作生物学研究室是沈其荣教授土壤微生物与有机肥团队下的一个关注于根际互作的研究小组。本小组由袁军副教授带领，主要关注：1.植物和微生物互作在抗病过程中的作用；2 环境微生物大数据整合研究；3 环境代谢组及其与微生物过程研究体系开发和应用。团队在过去三年中在 isme J， Microbiome， PCE，SBB，Horticulture Research等期刊上发表了多篇文章。欢迎关注 微生信生物 公众号对本研究小组进行了解。

团队工作及其成果 （点击查看）

团队关注

团队文章成果

团队成果-EasyStat专题

ggClusterNet专题

袁老师小小组

了解 交流 合作

团队成员邮箱 袁军：

[e-mail]；

文涛：

[e-mail]

团队公众号：