基因功能的注釋依賴于上一步的基因結構預測,根據預測結果從基因組上提取翻譯后的 蛋白序列 和主流的數據庫進行比對,完成功能注釋。常用數據庫一共有以幾種:
- Nr:NCBI官方非冗余蛋白數據庫,包括PDB, Swiss-Prot, PIR, PRF; 如果要用DNA序列,就是nt庫
- Pfam: 蛋白結構域注釋的分類系統(tǒng)
- Swiss-Prot: 高質量的蛋白數據庫,蛋白序列得到實驗的驗證
- KEGG: 代謝通路注釋數據庫.
- GO: 基因本體論注釋數據庫
除了以上幾個比較通用的數據庫外,其實還有很多小眾數據庫,應該根據課題研究和背景進行選擇。注意,數據庫本身并不能進行注釋,你只是通過序列相似性進行搜索,而返回的結果你稱之為注釋。因此數據庫和搜索工具要進行區(qū)分,所以你需要單獨下載數據庫和搜索工具,或者是同時下載包含數據庫和搜索工具的安裝包。
注意,后續(xù)分析中一定要保證你的蛋白序列中不能有代表氨基酸字符以外的字符,比如說有些軟件會把最后一個終止密碼子翻譯成"."或者"*"
BLASTP
這一部分用到的數據庫都是用BLASTP進行檢索,基本都是四步發(fā):下載數據庫,構建BLASTP索引,數據庫檢索,結果整理。其中結果整理需要根據BLASTP的輸出格式調整。
Nr的NCBI收集的最全的蛋白序列數據庫,但是無論是用NCBI的BLAST還是用速度比較快DIAMOND對nr進行搜索,其實都沒有利用好物種本身的信息。因此在RefSeq上下載對應物種的蛋白序列, 用BLASTP進行注釋即可。
# download
wget -4 -nd -np -r 1 -A *.faa.gz ftp://ftp.ncbi.nlm.nih.gov/refseq/release/plant/
mkdir -p ~/db/RefSeq
zcat *.gz > ~/db/RefSeq/plant.protein.faa
# build index
~/opt/biosoft/ncbi-blast-2.7.1+/bin/makeblastdb -in plant.protein.faa -dbtype prot -parse_seqids -title RefSeq_plant -out plant
# search
~/opt/biosoft/ncbi-blast-2.7.1+/bin/blastp -query protein.fa -out RefSeq_plant_blastp.xml -db ~/db/RefSeq/uniprot_sprot.fasta -evalue 1e-5 -outfmt 5 -num_threads 50 &
Swiss-Prot里收集了目前可信度最高的蛋白序列,一共有55w條記錄,數據量比較小,
# download
wget -4 -q ftp://ftp.uniprot.org/pub/databases/uniprot/current_release/knowledgebase/complete/uniprot_sprot.fasta.gz
gzip -d uniprot_sprot.fasta.gz
# builid index
~/opt/biosoft/ncbi-blast-2.7.1+/bin/makeblastdb -in uniprot_sprot.fasta -dbtype prot -title swiss_prot -parse_seqids
# search
~/opt/biosoft/ncbi-blast-2.7.1+/bin/blastp -query protein.fa -out swiss_prot.xml -db ~/db/swiss_prot/uniprot_sprot.fasta -evalue 1e-5 -outfmt 5 -num_threads 50 &
關于結果整理,已經有很多人寫了腳本,比如說我搜索BLAST XML CSV,就找到了https://github.com/Sunhh/NGS_data_processing/blob/master/annot_tools/blast_xml_parse.py, 所以就不過多介紹。
InterProScan
下面介紹的工具是InterProScan, 從它的9G的體量就可以感受它的強大之處,一次運行同時實現多個信息注釋。
- InterPro注釋
- Pfam數據庫注釋(可以通過hmmscan搜索pfam數據庫完成)
- GO注釋(可以基于NR和Pfam等數據庫,然后BLAST2GO完成,)
- Reactome通路注釋,不同于KEGG
命令如下
./interproscan-5.29-68.0/interproscan.sh -appl Pfam -f TSV -i sample.fa -cpu 50 -b sample -goterms -iprlookup -pa
-appl告訴軟件要執(zhí)行哪些數據分析,勾選的越多,分析速度越慢,Pfam就行。
KEGG
KEGG數據庫目前本地版收費,在線版收費,所以只能將蛋白序列在KEGG服務器上運行。因此你需要在http://www.genome.jp/tools/kaas/選擇合適的工具進行后續(xù)的分析。我上傳的50M大小蛋白序列,在KEGG服務器上只需要運行8個小時,也就是晚上提交任務,白天回來干活。
