kraken2bracken build db

作者: 胡童远 | 来源:发表于2023-08-22 09:50 被阅读0次

database is locked Possibly ther
私有API检测工具使用
IOS 新建项目忽略文件配置
mybatis-generator使用配置
uvm_component---super.xxx_phase
python: bracken merge samples
db/db小鼠
dB与dB※
CC2640R2F学习笔记（12）——修改发射功率
mysql新建db、复制db、删除db

kraken2 github:
https://github.com/DerrickWood/kraken2/blob/master/docs/MANUAL.markdown

环境

import argparse
import gzip
import itertools
import os
import pathlib
import textwrap
import sys
from Bio import Phylo

source /hwfsxx1/ST_HN/P18Z10200N0423/huty/software/miniconda3_2/etc/profile.d/conda.sh
conda activate py38
conda install -c bioconda biopython

GTDB注释

source /hwfssz1/**/hutongyuan/softwares/miniconda3/etc/profile.d/conda.sh
conda activate gtdbtk
# 配置数据库
echo "export GTDBTK_DATA_PATH=/hwfssz1/**/hutongyuan/databases/release207_v2/" > /hwfssz1/**/hutongyuan/softwares/miniconda3/envs/gtdbtk/etc/conda/activate.d/gtdbtk.sh
# 注释
gtdbtk classify_wf \
--cpus 32 \
-x fa \
--genome_dir bins_drep_comp70/dereplicated_genomes/ \
--out_dir anno_drep_comp70/

处理注释文件
很多潜在新物种没有注释信息，用MAG号作为补充
存在同种注释，种水平可全部加上MAG号区分

cp /hwfssz5/**/anno_drep_comp70/classify/gtdbtk.bac120.summary.tsv ./
cp /hwfssz5/**/anno_drep_comp70/classify/gtdbtk.ar53.summary.tsv ./
# 提取合并细菌，古菌
touch anno.bac.ar.tsv
cat gtdbtk.bac120.summary.tsv | sed '1d' | awk -F"\t" '{OFS="\t"}{print $1, $2}' >> anno.bac.ar.tsv
cat gtdbtk.ar53.summary.tsv | sed '1d' | awk -F"\t" '{OFS="\t"}{print $1, $2}' >> anno.bac.ar.tsv
# 补充未注释的部分
#!/usr/bin/env python3
import re,sys,os
with open("anno_edit.txt", 'w') as o:
    with open("anno.bac.ar.tsv", 'r') as infile:
        for line in infile:
            line = line.strip()
            line = re.split(r'[\t|;]', line)
            #if line[7] == "s__":
            line[7] = "s__" + line[0]
            if line[6] == "g__":
                line[6] = "g__" + line[0]
            if line[5] == "f__":
                line[5] = "f__" + line[0]
            if line[4] == "o__":
                line[4] = "o__" + line[0]
            if line[3] == "c__":
                line[3] = "c__" + line[0]
            o.write("{}\t{};{};{};{};{};{};{}\n".format(line[0],line[1],line[2],line[3],line[4],line[5],line[6],line[7]))

gtdb注释文件要求：
1 中间是制表符
2 末尾不留空行

kraken2建库准备

source /hwfssz1/ST_HEALTH/P18Z10200N0423/hutongyuan/softwares/miniconda3/etc/profile.d/conda.sh
conda activate python3.7

bins="/hwfssz5/ST_HEALTH/P18Z10200N0423/hutongyuan/analysis/fish/02_bin/bins_drep_comp70/dereplicated_genomes"
python ./tax_gtdb_final.py \
--gtdb anno_edit.txt \
--assemblies $bins/ \
--nodes nodes.dmp \
--names names.dmp \
--kraken_dir kraken_library

bash ./sc_buildk2db.sh

过程

Loading GTDB taxonomy:
    336 taxa found

Assigning tax IDs:
    domains:       2
    phyla:        13
    classes:      19
    orders:       38
    families:     56
    genera:       88
    species:     120

Writing taxonomy nodes to nodes.dmp
Writing taxonomy names to names.dmp

Searching for an assembly for each accession:
    120 / 120 assemblies found

Making Kraken assembly directory:
    120 / 120 assemblies

kraken_library里就是drep的MAG

建库

conda activate kraken2
for file in kraken_library/*.fa
do
    kraken2-build --add-to-library $file --db kraken_database
done

mkdir kraken_database/taxonomy
mv names.dmp kraken_database/taxonomy
mv nodes.dmp kraken_database/taxonomy

## kraken2-build --build --db kraken_database
vi build_kk2.sh

source /hwfssz1/ST_HEALTH/P18Z10200N0423/hutongyuan/softwares/miniconda3/etc/profile.d/conda.sh
conda activate kraken2
kraken2-build --build --db kraken_database --threads 4

qsub -cwd -l vf=10G,p=4 -q st.q -P P18Z10200N0423 -binding linear:4 build_kk2.sh

过程

Creating sequence ID to taxonomy ID map (step 1)...
Sequence ID to taxonomy ID map complete. [0.762s]
Estimating required capacity (step 2)...
Estimated hash table requirement: 1491144408 bytes
Capacity estimation complete. [3m54.774s]
Building database files (step 3)...
Taxonomy parsed and converted.
CHT created with 11 bits reserved for taxid.
Completed processing of 282544 sequences, 1399196925 bp
Writing data to disk...  complete.
Database files completed. [9m25.653s]
Database construction complete. [Total: 13m35.437s]

generate bracken distribution

# Let's continue to generate bracken distribution for profile correction
cd kraken_database
kraken2 --db=./ --threads=4 <( find -L library \( -name "*.fna" -o -name "*.fasta" -o -name "*.fa" \) -exec cat {} + ) > database.kraken

过程

Loading database information... done.
41646 sequences (388.85 Mbp) processed in 37.970s (65.8 Kseq/m, 614.45 Mbp/m).
  41646 sequences classified (100.00%)
  0 sequences unclassified (0.00%)

readlen=50、100、200、300（这一步是后面bracken可能会用到的）

conda activate bracken
kmer2read_distr --seqid2taxid seqid2taxid.map --taxonomy ./taxonomy --kraken database.kraken --output database50mers.kraken -l 50 -t 1
kmer2read_distr --seqid2taxid seqid2taxid.map --taxonomy ./taxonomy --kraken database.kraken --output database100mers.kraken -l 100 -t 1
kmer2read_distr --seqid2taxid seqid2taxid.map --taxonomy ./taxonomy --kraken database.kraken --output database200mers.kraken -l 200 -t 1
kmer2read_distr --seqid2taxid seqid2taxid.map --taxonomy ./taxonomy --kraken database.kraken --output database300mers.kraken -l 300 -t 1

generate_kmer_distribution.py -i database50mers.kraken -o database50mers.kmer_distrib
generate_kmer_distribution.py -i database100mers.kraken -o database100mers.kmer_distrib
generate_kmer_distribution.py -i database200mers.kraken -o database200mers.kmer_distrib
generate_kmer_distribution.py -i database300mers.kraken -o database300mers.kmer_distrib

# rm -r ../kraken_library/

过程

(bracken) [hutongyuan@cngb-xcompute-0-18 kraken_database]$ kmer2read_distr --seqid2taxid seqid2taxid.map --taxonomy ./taxonomy --kraken database.kraken --output database50mers.kraken -l 50 -t 1
        >>STEP 0: PARSING COMMAND LINE ARGUMENTS
                Taxonomy nodes file: ./taxonomy/nodes.dmp
                Seqid file:          seqid2taxid.map
                Num Threads:         1
                Kmer Length:         31
                Read Length:         50
        >>STEP 1: READING SEQID2TAXID MAP
                282544 total sequences read
        >>STEP 2: READING NODES.DMP FILE
                1056 total nodes read
        >>STEP 3: CONVERTING KMER MAPPINGS INTO READ CLASSIFICATIONS:
                50mers, with a database built using 31mers
                282546 sequences converted (finished: )oncoct_E9_bin.21_k119_17748|kraken:taxid|595)97)))
        Time Elaped: 1 minutes, 23 seconds, 0.00000 microseconds
        =============================
(bracken) [hutongyuan@cngb-xcompute-0-18 kraken_database]$ kmer2read_distr --seqid2taxid seqid2taxid.map --taxonomy ./taxonomy --kraken database.kraken --output database100mers.kraken -l 100 -t 1
        >>STEP 0: PARSING COMMAND LINE ARGUMENTS
                Taxonomy nodes file: ./taxonomy/nodes.dmp
                Seqid file:          seqid2taxid.map
                Num Threads:         1
                Kmer Length:         31
                Read Length:         100
        >>STEP 1: READING SEQID2TAXID MAP
                282544 total sequences read
        >>STEP 2: READING NODES.DMP FILE
                1056 total nodes read
        >>STEP 3: CONVERTING KMER MAPPINGS INTO READ CLASSIFICATIONS:
                100mers, with a database built using 31mers
                282546 sequences converted (finished: )oncoct_E9_bin.21_k119_17748|kraken:taxid|595)97)))
        Time Elaped: 2 minutes, 29 seconds, 0.00000 microseconds
        =============================
(bracken) [hutongyuan@cngb-xcompute-0-18 kraken_database]$ kmer2read_distr --seqid2taxid seqid2taxid.map --taxonomy ./taxonomy --kraken database.kraken --output database200mers.kraken -l 200 -t 1
        >>STEP 0: PARSING COMMAND LINE ARGUMENTS
                Taxonomy nodes file: ./taxonomy/nodes.dmp
                Seqid file:          seqid2taxid.map
                Num Threads:         1
                Kmer Length:         31
                Read Length:         200
        >>STEP 1: READING SEQID2TAXID MAP
                282544 total sequences read
        >>STEP 2: READING NODES.DMP FILE
                1056 total nodes read
        >>STEP 3: CONVERTING KMER MAPPINGS INTO READ CLASSIFICATIONS:
                200mers, with a database built using 31mers
                282546 sequences converted (finished: )oncoct_E9_bin.21_k119_17748|kraken:taxid|595)97)))
        Time Elaped: 4 minutes, 12 seconds, 0.00000 microseconds
        =============================
(bracken) [hutongyuan@cngb-xcompute-0-18 kraken_database]$ kmer2read_distr --seqid2taxid seqid2taxid.map --taxonomy ./taxonomy --kraken database.kraken --output database300mers.kraken -l 300 -t 1
        >>STEP 0: PARSING COMMAND LINE ARGUMENTS
                Taxonomy nodes file: ./taxonomy/nodes.dmp
                Seqid file:          seqid2taxid.map
                Num Threads:         1
                Kmer Length:         31
                Read Length:         300
        >>STEP 1: READING SEQID2TAXID MAP
                282544 total sequences read
        >>STEP 2: READING NODES.DMP FILE
                1056 total nodes read
        >>STEP 3: CONVERTING KMER MAPPINGS INTO READ CLASSIFICATIONS:
                300mers, with a database built using 31mers
                282546 sequences converted (finished: )oncoct_E9_bin.21_k119_17748|kraken:taxid|595)97)))
        Time Elaped: 5 minutes, 39 seconds, 0.00000 microseconds
        =============================

(bracken) [hutongyuan@cngb-login-0-8 kraken_database]$ generate_kmer_distribution.py -i database50mers.kraken -o database50mers.kmer_distrib
PROGRAM START TIME: 12-04-2022 14:10:33
...473 total genomes read from kraken output file
...creating kmer counts file -- lists the number of kmers of each classification per genome
...creating kmer distribution file -- lists genomes and kmer counts contributing to each genome
PROGRAM END TIME: 12-04-2022 14:10:35
(bracken) [hutongyuan@cngb-login-0-8 kraken_database]$ generate_kmer_distribution.py -i database100mers.kraken -o database100mers.kmer_distrib
PROGRAM START TIME: 12-04-2022 14:10:35
...473 total genomes read from kraken output file
...creating kmer counts file -- lists the number of kmers of each classification per genome
...creating kmer distribution file -- lists genomes and kmer counts contributing to each genome
PROGRAM END TIME: 12-04-2022 14:10:37
(bracken) [hutongyuan@cngb-login-0-8 kraken_database]$ generate_kmer_distribution.py -i database200mers.kraken -o database200mers.kmer_distrib
PROGRAM START TIME: 12-04-2022 14:10:37
...473 total genomes read from kraken output file
...creating kmer counts file -- lists the number of kmers of each classification per genome
...creating kmer distribution file -- lists genomes and kmer counts contributing to each genome
PROGRAM END TIME: 12-04-2022 14:10:38
(bracken) [hutongyuan@cngb-login-0-8 kraken_database]$ generate_kmer_distribution.py -i database300mers.kraken -o database300mers.kmer_distrib
PROGRAM START TIME: 12-04-2022 14:10:39
...473 total genomes read from kraken output file
...creating kmer counts file -- lists the number of kmers of each classification per genome
...creating kmer distribution file -- lists genomes and kmer counts contributing to each genome
PROGRAM END TIME: 12-04-2022 14:10:40

database is locked Possibly ther
报错： unable to attach DB: error: accessing build database ...
私有API检测工具使用
1、下载已经build好的ios_private.db库，将 ios_private.db 放入到项目的根目录，主...
IOS 新建项目忽略文件配置
.DS_Store/build/*.pbxuser.db!default.pbxuser*.mode1v3!def...
mybatis-generator使用配置
添加依赖，配置mybatis-generator build.gradle 再resource下添加配置文件 db...
uvm_component---super.xxx_phase
对于build_phase来说，uvm_component对其做的最重要的事情就是自动获取通过config_db：...
python: bracken merge samples
用了类似构造泛基因组表格的方法 1 kraken2bracken分析 kraken2数据处理 bracken数据处...
db/db小鼠
db/db小鼠为Leptin receptor基因缺陷鼠。模型简介：Leptin receptor（LEPR）是 ...
dB与dB※
一言以蔽之是相对值的单位，是绝对值的单位。 dB 是表达两个功率（或者场量）比值关系的单位。功率比的关系为只知道...
CC2640R2F学习笔记（12）——修改发射功率
一、背景 1.1 发射功率范围 (最大)5dB、4dB、3dB、2dB、1dB、（默认）0dB、-3dB、-6dB...
mysql新建db、复制db、删除db
新建数据库，修改类型为utf8 create database abc;use abc;ALTER DATABAS...