最近偶然从一篇PNAs的文章，看到作者不仅进行了16S rRNA测序，而且通过marker基因和宏基因组bins同时计算了群落结构，对比三者发现，菌群结构相似。
看了一下marker基因对宏基因组reads进行分类的方法，作者提到了singlem（见下图）。
于是就想试试singlem好不好用，按照作者的说法，为了与宏基因组bins的分类尽可能一致，所以要使用GTDB的分类标准。其实很有道理，因为对于细菌来说，GTDB的分类和NCBI里面分类好多都不一样。而后作者又提到，他download了所有GTDB的rplp序列，这个我是没找到从GTDB下载单个基因的方法， GTDB难道不是只能下载基因组和某个物种串联后的蛋白文件吗？
所以就干脆自己翻译，自己用hmmsearch找rplp基因。
先说singlem的安装，其实有不少坑，首先它依赖于GraftM。GraftM用conda安装的时候，经常会漏掉一些python的包，其实问题不大，耐心点用pip一一安装完毕。
下面说说操作步骤
1 使用prodigal预测GTDB所有基因组
prodigal -i ${j%.gz} -a temp.orfs.faa -d temp.orfs.ffn -m -o temp.txt -p meta
cut -f 1 -d " " temp.orfs.faa >${j%.fna.gz}.faa
cut -f 1 -d " " temp.orfs.ffn >${j%.fna.gz}.ffn
建议写个shell循环，毕竟6万多个基因组，我这只是贴了脚本的一部分。
2 pfam下载rplp基因hmm文件，即Ribosomal_L16.hmm。而后hmmsearch继续循环查找
3 继续写个小shell循环提取所有的，能找到的rplp基因序列
grep -v "^#" rplp.out | awk '{print $1}' | seqkit grep -f - ${line%/}.faa > ${line%_genomic/}.rplp
里面还涉及到修改序列名字等，自行写脚本，此处不表。
4 合并rplp基因，以及GTDB库的taxonomy文件，使用GraftM建库
graftM create --output ./my.gpkg --sequences gtdb.rplp.fasta --taxonomy taxonomy.txt
耐心等待完成后，将my.gpkg移动到miniconda3/envs/singlem/lib/python3.6/site-packages/singlem/data/
这里还有个坑，注意json文件，修改成自定义库的名字，否则python找不到。
此时建库完成，参考singlem说明书，跑一下试试
singlem pipe --forward 58_clean_1.fq --reverse 58_clean_2.fq --otu_table 2.tsv --threads 20
