最近偶然從一篇PNAs的文章，看到作者不僅進(jìn)行了16S rRNA測序，而且通過marker基因和宏基因組bins同時計算了群落結(jié)構(gòu)，對比三者發(fā)現(xiàn)，菌群結(jié)構(gòu)相似。

看了一下marker基因?qū)昊蚪Mreads進(jìn)行分類的方法，作者提到了singlem（見下圖）。

于是就想試試singlem好不好用，按照作者的說法，為了與宏基因組bins的分類盡可能一致，所以要使用GTDB的分類標(biāo)準(zhǔn)。其實很有道理，因為對于細(xì)菌來說，GTDB的分類和NCBI里面分類好多都不一樣。而后作者又提到，他download了所有GTDB的rplp序列，這個我是沒找到從GTDB下載單個基因的方法， GTDB難道不是只能下載基因組和某個物種串聯(lián)后的蛋白文件嗎？

所以就干脆自己翻譯，自己用hmmsearch找rplp基因。

先說singlem的安裝，其實有不少坑，首先它依賴于GraftM。GraftM用conda安裝的時候，經(jīng)常會漏掉一些python的包，其實問題不大，耐心點用pip一一安裝完畢。

下面說說操作步驟

1 ?使用prodigal預(yù)測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循環(huán)，畢竟6萬多個基因組，我這只是貼了腳本的一部分。

2 pfam下載rplp基因hmm文件，即Ribosomal_L16.hmm。而后hmmsearch繼續(xù)循環(huán)查找

3 繼續(xù)寫個小shell循環(huán)提取所有的，能找到的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