5-跟着science学习宏基因组-kraken物种注释-微生信生物-科易网技术创新

写在前面

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文件夹中的数据进行注释

安装kraken

conda 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：分类名称

如果要构建自定义数据库，请参考：#standard-kraken-database

这里考虑到大家的内存

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

kraken-build --shrink 10000 --db ~/db/minikraken --new-db minikraken