最近偶然從一篇PNAs的文章，看到作者不僅進(jìn)行了16S rRNA測(cè)序逗柴，而且通過(guò)marker基因和宏基因組bins同時(shí)計(jì)算了群落結(jié)構(gòu)蛹头，對(duì)比三者發(fā)現(xiàn)，菌群結(jié)構(gòu)相似戏溺。
看了一下marker基因?qū)昊蚪Mreads進(jìn)行分類(lèi)的方法渣蜗，作者提到了singlem（見(jiàn)下圖）。
于是就想試試singlem好不好用于购，按照作者的說(shuō)法袍睡，為了與宏基因組bins的分類(lèi)盡可能一致，所以要使用GTDB的分類(lèi)標(biāo)準(zhǔn)肋僧。其實(shí)很有道理斑胜，因?yàn)閷?duì)于細(xì)菌來(lái)說(shuō)控淡，GTDB的分類(lèi)和NCBI里面分類(lèi)好多都不一樣。而后作者又提到止潘，他download了所有GTDB的rplp序列掺炭，這個(gè)我是沒(méi)找到從GTDB下載單個(gè)基因的方法， GTDB難道不是只能下載基因組和某個(gè)物種串聯(lián)后的蛋白文件嗎凭戴？
所以就干脆自己翻譯涧狮，自己用hmmsearch找rplp基因。
先說(shuō)singlem的安裝么夫，其實(shí)有不少坑者冤，首先它依賴于GraftM。GraftM用conda安裝的時(shí)候档痪，經(jīng)常會(huì)漏掉一些python的包涉枫，其實(shí)問(wèn)題不大，耐心點(diǎn)用pip一一安裝完畢腐螟。
下面說(shuō)說(shuō)操作步驟
1 ?使用prodigal預(yù)測(cè)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
建議寫(xiě)個(gè)shell循環(huán)愿汰，畢竟6萬(wàn)多個(gè)基因組，我這只是貼了腳本的一部分乐纸。
2 pfam下載rplp基因hmm文件衬廷，即Ribosomal_L16.hmm。而后hmmsearch繼續(xù)循環(huán)查找
3 繼續(xù)寫(xiě)個(gè)小shell循環(huán)提取所有的汽绢，能找到的rplp基因序列
grep -v "^#" rplp.out | awk '{print $1}' | seqkit grep -f - ${line%/}.faa > ${line%_genomic/}.rplp
里面還涉及到修改序列名字等吗跋，自行寫(xiě)腳本，此處不表宁昭。
4 合并rplp基因小腊，以及GTDB庫(kù)的taxonomy文件，使用GraftM建庫(kù)
graftM create --output ./my.gpkg --sequences gtdb.rplp.fasta --taxonomy taxonomy.txt
耐心等待完成后久窟，將my.gpkg移動(dòng)到miniconda3/envs/singlem/lib/python3.6/site-packages/singlem/data/
這里還有個(gè)坑，注意json文件本缠，修改成自定義庫(kù)的名字斥扛，否則python找不到。
此時(shí)建庫(kù)完成丹锹，參考singlem說(shuō)明書(shū)稀颁，跑一下試試
singlem pipe --forward 58_clean_1.fq --reverse 58_clean_2.fq --otu_table 2.tsv --threads 20
