目前有兩種方式可計算宏基因組基因的豐度昆婿，一種是基于比對生宛，比如bwa,bowtie,soapaligner等主流的比對軟件，另一種是不比對快速估計基因豐度，可以用近倆年來流行的salmon軟件馁痴，這個軟件在轉錄組的數(shù)據(jù)比對中也經(jīng)常用到，可以直接計算出原始的Counts值和標準化的TPM值订雾，此外由于是基于非比對晴玖，計算的速度得到很大的提升，同時也節(jié)省了部分的內(nèi)存(少了龐大的sam/bam文件)谨娜，可以說是多快好省航攒，但是目前的高分文章中也還是不少用基于比對的方法去計算宏基因組的基因豐度的，尤其是Soapaligner比對軟件趴梢，是由華大公司開發(fā)的Soap系列的工具漠畜，首先在速度上具有較大的優(yōu)勢，同時也很好的解決了在小內(nèi)存情況下將短序列比對到參考基因組上的問題坞靶，作用也不容小覷憔狞。當然也有不少研究人員在前期質(zhì)控環(huán)節(jié)進行去宿主也是不錯的選擇~下面我就分別簡單介紹一下基于比對的soapaligner和不比對快速估計的samlon倆個軟件的操作流程！彰阴！
Option1 SOAPaligenr
安裝與比對
1. 安裝
conda install soap
2. 構建索引
2bwt-builder contig.fa
會產(chǎn)生13個索引文件
3.雙端序列進行比對
soap -p 6 -r 2 -m 200 -x 400 -M 4 -a ./cleandata/B1.clean.fq1.gz -b ./sample1/B1.clean.fq2.gz -D contig.index -o B1_PE.soap -2 B1_SE.soap -u B1_UN.fasta
這個軟件的參數(shù)有點多瘾敢，-a,-b分別輸入雙端序列, -D 索引文件的地址，-o,輸出能比對上組裝后的參考序列且reads成對存在的比對文件尿这，-2簇抵，能比對上組裝后的參考序列但reads不成對存在的比對文件；-u射众，不能比對上組裝后的參考序列碟摆。至于-m,-x參數(shù)比較主觀，看過很多文獻都有不同的設定责球，這里貼出一篇文獻供參考(DOI:10.3389/fmicb.2017.01546),文獻中的參數(shù)為?m 4 ?r 2 ?m 100 ?x 1000＝孤摹！
標準化
看過十幾篇文獻在基因豐度標準化這一塊都引用到了一篇2012年發(fā)表在Nature上的一篇文獻（qin2012,doi:10.1038/nature11450）,看過之后在該文章的補充材料里發(fā)現(xiàn)了計算基因豐度的步驟雏逾，如下圖倆步驟嘉裤，看過之后感覺其實和轉錄組中的RPKM/TPM等標準化的思路較為類似，先標準化基因的長度栖博，然后再除以標準化的值的總和屑宠。還是挺簡單的嘛。
我們也按照這個思路根據(jù)soapaligner比對的結果來進行標準化,在結果文件中我們能得到三個文件仇让，分別為：B1_PE.soap:能比對上組裝后的參考序列且reads成對存在的比對文件典奉；B1_SE.soap: 能比對上組裝后的參考序列但reads不成對存在的比對文件躺翻；sample1_UN.fasta：不能比對上組裝后的參考序列，通常為fasta格式卫玖；我們只需要比對上的文件B1_PE.soap和B1_SE.soap即可公你。這兩個文件的第8列中比對上的基因，我們將其取出來sort|uniq -c 就可以計算出基因的原始copy數(shù)目假瞬，也就是比對上特定基因的數(shù)目陕靠。之后再計算出比對基因的長度，也就是你組裝去冗余之后的非冗余基因集的各基因長度脱茉，這步可以寫個小腳本或者用現(xiàn)成的輪子seqkit/seqtk都可以很快的解決剪芥，最后利用上面的公式進行標準化就行了~ 還需要注意的一點，沒有比對上的基因我們就用0代替他的原始豐度就好了琴许。
Option2 Samlon
這種方法就更加簡單了税肪，也是大致兩步走，建立索引榜田，序列比對益兄，只不過時間將大大縮短，這種方法網(wǎng)上的教程有很多串慰，我這里貼出我的代碼大家可以比較一下使用~ 建議大家使用時候都安裝最新版本的軟件偏塞，github地址為：https://github.com/COMBINE-lab/salmon/releases/
#1.建立索引
salmon index -t B1_NR100nl.fasta -p 9 -k 31 -i ./index
#2.比對
salmon quant --validateMappings -i ./index -l A -p 3 --meta -1 ../clean1_$.fastq -2 ../clean2_$邦鲫.fastq -o $.quant #基因定量
注意的就是宏基因組數(shù)據(jù)加上-p meta這個參數(shù)神汹，其他參數(shù)默認就好~庆捺，得到的結果有個quant.sf文件，豐度信息就在里面屁魏。為5列的信息,第4列就是我們想要的標準化后的基因豐度了
OK滔以，拿到基因豐度文件之后，我們就可以用來進行后續(xù)分析氓拼，比如計算eggnog數(shù)據(jù)庫比對后的功能基因豐度~這個留到下次再說吧Ｄ慊！
下面是剛創(chuàng)建個人公眾號桃漾，會定時更新R坏匪、linux、python撬统，組學方面的學習內(nèi)容适滓，請多多支持呦~~
