单菌基因预测、去冗余和定量
nrMAGs 的基因组注释: MAGs 的基因组注释首先使用 Prokka （v1.13） 使用 metaWRAP 的 annotate_bins 模块进行。然后使用 MicrobeAnnotator （v2.0.5） 处理注释的基因组，以进行功能注释并计算 KEGG 模块完整性。使用 Kofamsca在精选的 KEGG 直系同源物 （KO） 数据库中搜索所有蛋白质;根据 Kofamscan 的自适应分数阈值选择最佳匹配项。提取没有 KO 标识符（或匹配）的蛋白质，并针对其他数据库（例如 Swissprot、精选的 RefSeq 数据库或非精选的 trEMBL 数据库）进行搜索51。提取与每个基因组（或蛋白质集）中所有蛋白质相关的 KO 标识符，并根据模块中的总步骤、每个步骤所需的蛋白质 （KO） 以及每个基因组中存在的 KO 来计算 KEGG 模块完整性。最后，将结果编译成所有基因组的单个矩阵样模块完整性表。
为了闭合 EGF1-FE1 的基因组，使用 Unicycler34 （v.0.4.8） 进行了 Illumina 短读长和 Oxford Nanopore 长读长的杂交组装。使用 FastQC35 （v.0.11.9） 对 Illumina 短读长进行质量控制检查，使用 Trimmomatic36 （v.0.36） 进行修整和过滤，以去除 Nextera 接头和低质量读长。对于 ONT 读数，使用 NanoStat37 （v.1.1.2） 检查原始读数的质量，然后使用 Filtlong （v.0.2.0） 过滤，以仅提取最佳读数，直到总数碱基为 1.85 亿（50× ONT 覆盖率）。在 Unicycler 中，默认设置用于混合组装。E. gallinarum 的封闭基因组在 RAST server38,39 （v.2.0） 和 Prokka40 （v.1.14.5） 中进行了注释。
# 输入文件：拼装好的序列文件 result/megahit/final.contigs.fa # 输出文件：prodigal预测的基因序列 temp/prodigal/gene.fa # 基因文件大，可参考附录prodigal拆分基因文件，并行计算 mkdir -p temp/prodigal # prodigal的meta模式预测基因，35s，>和2>&1记录分析过程至gene.log prodigal -i result/megahit/final.contigs.fa \ -d temp/prodigal/gene.fa \ -o temp/prodigal/gene.gff \ -p meta -f gff > temp/prodigal/gene.log 2>&1 # 查看日志是否运行完成，有无错误 tail temp/prodigal/gene.log # 统计基因数量 seqkit stat temp/prodigal/gene.fa # 统计完整基因数量，数据量大可只用完整基因部分 grep -c 'partial=00' temp/prodigal/gene.fa # 提取完整基因(完整片段获得的基因全为完整，如成环的细菌基因组) grep 'partial=00' temp/prodigal/gene.fa | cut -f1 -d ' '| sed 's/>//' > temp/prodigal/full_length.id seqkit grep -f temp/prodigal/full_length.id temp/prodigal/gene.fa > temp/prodigal/full_length.fa seqkit stat temp/prodigal/full_length.fa
# 输入文件：prodigal预测的基因序列 temp/prodigal/gene.fa # 输出文件：去冗余后的基因和蛋白序列：result/NR/nucleotide.fa, result/NR/protein.fa mkdir -p result/NR # aS覆盖度，c相似度，G局部比对，g最优解，T多线程，M内存0不限制 # 2万基因2m，2千万需要2000h，多线程可加速 cd-hit-est -i temp/prodigal/gene.fa \ -o result/NR/nucleotide.fa \ -aS 0.9 -c 0.95 -G 0 -g 0 -T 0 -M 0 # 统计非冗余基因数量，单次拼接结果数量下降不大，多批拼接冗余度高 grep -c '>' result/NR/nucleotide.fa # 翻译核酸为对应蛋白序列, --trim去除结尾的* seqkit translate --trim result/NR/nucleotide.fa \ > result/NR/protein.fa # 两批数据去冗余使用cd-hit-est-2d加速，见附录
# 输入文件：去冗余后的基因和蛋白序列：result/NR/nucleotide.fa # 输出文件：Salmon定量后的结果：result/salmon/gene.count, gene.TPM mkdir -p temp/salmon salmon -v # 1.4.0 # 建索引, -t序列, -i 索引，10s salmon index \ -t result/NR/nucleotide.fa \ -p 9 \ -i temp/salmon/index # 定量，l文库类型自动选择，p线程，--meta宏基因组模式, 2个任务并行2个样 # 注意parallel中待并行的命令必须是双引号，内部变量需要使用原始绝对路径 tail -n+2 result/metadata.txt|cut -f1|rush -j 2 \ "salmon quant \ -i temp/salmon/index -l A -p 3 --meta \ -1 temp/qc/{1}_1.fastq \ -2 temp/qc/{1}_2.fastq \ -o temp/salmon/{1}.quant" # 合并 mkdir -p result/salmon salmon quantmerge --quants temp/salmon/*.quant \ -o result/salmon/gene.TPM salmon quantmerge --quants temp/salmon/*.quant \ --column NumReads -o result/salmon/gene.count sed -i '1 s/.quant//g' result/salmon/gene.* # 预览结果表格 head -n3 result/salmon/gene.*
# 输入数据：上一步预测的蛋白序列 result/NR/protein.fa # 中间结果：temp/eggnog/protein.emapper.seed_orthologs # temp/eggnog/output.emapper.annotations # temp/eggnog/output # COG定量表：result/eggnog/cogtab.count # result/eggnog/cogtab.count.spf (用于STAMP) # KO定量表：result/eggnog/kotab.count # result/eggnog/kotab.count.spf (用于STAMP) # CAZy碳水化合物注释和定量：result/dbcan2/cazytab.count # result/dbcan2/cazytab.count.spf (用于STAMP) # 抗生素抗性：result/resfam/resfam.count # result/resfam/resfam.count.spf (用于STAMP) # 这部分可以拓展到其它数据库
对prokka预测基因生成的faa文件使用eggNOG注释
emapper.py \ --data_dir /data/wangyang/databases/eggnog-mapper-data \ --itype proteins \ -m diamond \ --output out_eggnog \ --cpu 8 \ -i /ssd1/wy/workspace2/leipu/leipu_workspace2/output/prokka/OCF-1-PACBIO-HIFI/OCF-1-PACBIO-HIFI.faa
# https://github.com/eggnogdb/eggnog-mapper/wiki/eggNOG-mapper-v2 emapper.py --version # 2.1.6 # diamond比对基因至eggNOG 5.0数据库, 9p11m, 1~9h，默认diamond 1e-3 mkdir -p temp/eggnog time emapper.py --no_annot --no_file_comments --override \ --data_dir ${db}/eggnog \ -i result/NR/protein.fa \ --cpu 9 -m diamond \ -o temp/eggnog/protein # 比对结果功能注释, 1h # sqlite3.OperationalError: no such table: prots是数据库不配套，重新下载即可 emapper.py \ --annotate_hits_table temp/eggnog/protein.emapper.seed_orthologs \ --data_dir ${db}/eggnog \ --cpu 9 --no_file_comments --override \ -o temp/eggnog/output # 2.1较2.0结果又有新变化，添加了#号表头，减少了列 sed '1 s/^#//' temp/eggnog/output.emapper.annotations > temp/eggnog/output csvtk -t headers -v temp/eggnog/output summarizeAbundance生成COG/KO/CAZy丰度汇总表 mkdir -p result/eggnog # 显示帮助，需要Python3环境，可修改软件第一行指定python位置，如指定某Python执行脚本 /mnt/bai/yongxin/miniconda2/envs/humann3/bin/python3 /db/EasyMicrobiome/script/summarizeAbundance.py summarizeAbundance.py -h # 汇总，7列COG_category按字母分隔，12列KEGG_ko和19列CAZy按逗号分隔，原始值累加 # 指定humann3中的Python 3.7.6运行正常，qiime2中的Python 3.6.13报错 summarizeAbundance.py \ -i result/salmon/gene.TPM \ -m temp/eggnog/output \ -c '7,12,19' -s '*+,+,' -n raw \ -o result/eggnog/eggnog sed -i 's/^ko://' result/eggnog/eggnog.KEGG_ko.raw.txt sed -i '/^-/d' result/eggnog/eggnog* # eggnog.CAZy.raw.txt eggnog.COG_category.raw.txt eggnog.KEGG_ko.raw.txt # 添加注释生成STAMP的spf格式，结合metadata.txt进行差异比较 awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \ /db/EasyMicrobiome/kegg/KO_description.txt \ result/eggnog/eggnog.KEGG_ko.raw.txt | \ sed 's/^\t/Unannotated\t/' \ > result/eggnog/eggnog.KEGG_ko.TPM.spf # KO to level 1/2/3 summarizeAbundance.py \ -i result/eggnog/eggnog.KEGG_ko.raw.txt \ -m /db/EasyMicrobiome/kegg/KO1-4.txt \ -c 2,3,4 -s ',+,+,' -n raw \ -o result/eggnog/KEGG # CAZy awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \ /db/EasyMicrobiome/dbcan2/CAZy_description.txt result/eggnog/eggnog.CAZy.raw.txt | \ sed 's/^\t/Unannotated\t/' > result/eggnog/eggnog.CAZy.TPM.spf # COG awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2"\t"$3} NR>FNR{print a[$1],$0}' \ /db/EasyMicrobiome/eggnog/COG.anno result/eggnog/eggnog.COG_category.raw.txt > \ result/eggnog/eggnog.COG_category.TPM.spf
# 比对CAZy数据库, 用时2~18m mkdir -p temp/dbcan2 # --sensitive慢10倍，dbCAN2推荐e值为1e-102，此处结果3条太少，以1e-3为例演示 diamond blastp \ --db /db/dbcan2/CAZyDB.09242021 \ --query result/NR/protein.fa \ --threads 9 -e 1e-3 --outfmt 6 --max-target-seqs 1 --quiet \ --out temp/dbcan2/gene_diamond.f6 wc -l temp/dbcan2/gene_diamond.f6 # 整理比对数据为表格 mkdir -p result/dbcan2 # 提取基因与dbcan分类对应表 format_dbcan2list.pl \ -i temp/dbcan2/gene_diamond.f6 \ -o temp/dbcan2/gene.list # 按对应表累计丰度，依赖 summarizeAbundance.py \ -i result/salmon/gene.TPM \ -m temp/dbcan2/gene.list \ -c 2 -s ',' -n raw \ -o result/dbcan2/TPM # 添加注释生成STAMP的spf格式，结合metadata.txt进行差异比较 awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$2} NR>FNR{print a[$1],$0}' \ /db/EasyMicrobiome/dbcan2/CAZy_description.txt result/dbcan2/TPM.CAZy.raw.txt | \ sed 's/^\t/Unannotated\t/' > result/dbcan2/TPM.CAZy.raw.spf # 检查未注释数量，有则需要检查原因 # grep 'Unannotated' result/dbcan2/TPM.CAZy.raw.spf|wc -l
数据库：https://card.mcmaster.ca/ ，有在线分析平台，本地代码供参考
# 参考文献：http://doi.org/10.1093/nar/gkz935 # 软件使用Github: https://github.com/arpcard/rgi # 启动rgi环境 conda activate rgi rgi -h # 5.2.1 # 蛋白注释 mkdir -p result/card cut -f 1 -d ' ' result/NR/protein.fa > temp/protein.fa rgi main -i temp/protein.fa -t protein \ -n 9 -a DIAMOND --include_loose --clean \ -o result/card/protein
结果说明：
# Generate report in default taxid output conda activate meta memusg -t kraken2 --db /db/kraken2/mini \ result/NR/nucleotide.fa \ --threads 3 \ --report temp/NRgene.report \ --output temp/NRgene.output # Genes & taxid list grep '^C' temp/NRgene.output|cut -f 2,3|sed '1 i Name\ttaxid' \ > temp/NRgene.taxid # Add taxonomy awk 'BEGIN{FS=OFS="\t"} NR==FNR{a[$1]=$0} NR>FNR{print $1,a[$2]}' \ /db/EasyMicrobiome/kraken2/taxonomy.txt \ temp/NRgene.taxid \ > result/NR/nucleotide.tax memusg -t /conda2/envs/humann3/bin/python3 /db/EasyMicrobiome/script/summarizeAbundance.py \ -i result/salmon/gene.TPM \ -m result/NR/nucleotide.tax \ -c '2,3,4,5,6,7,8,9' -s ',+,+,+,+,+,+,+,' -n raw \ -o result/NR/tax wc -l result/NR/tax*|sort -n