基因家族專(zhuān)題(3)：基因家族成員的鑒定

Data preparation

繼續(xù)上次的內(nèi)容，下載好數(shù)據(jù)后就可以正式開(kāi)始鑒定了。首先回顧一下，下載好的數(shù)據(jù)。

  1. 基因組序列信息，存儲(chǔ)基因組序列信息的.fasta文件。還有其蛋白質(zhì)序列，也是以.fasta結(jié)尾的文件。一般來(lái)說(shuō)注釋的比較好的基因組都會(huì)含有這些文件。
  2. 基因組基因結(jié)構(gòu)注釋信息。儲(chǔ)存基因的intron，exon，CDS,gene等坐標(biāo)信息的.gff3或.gtf文件。
  3. 所感興趣的基因家族隱馬可夫模型，hmm文件
-rw-r----- 1 hhu pawsey0149  9738306 Oct 25 12:24 Arabidopsis_thaliana.TAIR10.41.gff3.gz
-rw-r----- 1 hhu pawsey0149 14623390 Oct 25 12:23 Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
-rw-r----- 1 hhu pawsey0149 36462703 Oct 25 12:23 Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
-rw-r----- 1 hhu pawsey0149  9776319 Oct 25 12:22 Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
-rw-r----- 1 hhu pawsey0149   118379 Oct 25 12:26 NBS-ARC.hmm

基因家族鑒定的工具h(yuǎn)mmer

一般尋找基因家族，都可以通過(guò)保守結(jié)構(gòu)域來(lái)預(yù)測(cè)，從而找到物種的某一基因家族，從而進(jìn)行之后的分析。這里就需要用到HMMER，來(lái)鑒定物種某一基因家族。

HMMER3.1下載地址：http://hmmer.org/download.html
HMMER3.1 manual：http://eddylab.org/software/hmmer3/3.1b2/Userguide.pdf

hmmbuild/hmmsearch/hmmscan/hmmalign 這幾個(gè)功能是主要用于蛋白質(zhì)結(jié)構(gòu)與分析和注釋的hmmer中小工具

在鑒定基因家族時(shí)，常用到的工具是hmmsearch，里面常用的算法有三種。一般我們使用--cut_tc算法對(duì)隱馬可夫模型進(jìn)行搜索，tc算法是使用pfam提供的hmm文件中trusted cutoof的值進(jìn)行篩選，相對(duì)比較可靠。

  --cut_ga : use profile's GA gathering cutoffs to set all thresholding
  --cut_nc : use profile's NC noise cutoffs to set all thresholding
  --cut_tc : use profile's TC trusted cutoffs to set all thresholding

具體實(shí)戰(zhàn)操作

下面會(huì)根據(jù)一篇經(jīng)典文獻(xiàn)中的方法，對(duì)擬南芥進(jìn)行NBS-LRR基因組的探索。首先，回顧一下文獻(xiàn)看看整體它的思路和方法。

Identification of NBS-LRR genes

Predicted proteins from the cassava genome were scanned using HMMER v3 [39] using the Hidden Markov Model (HMM) corresponding to the Pfam [40] NBS (NB-ARC) family (PF00931; http://pfam.sanger.ac.uk/). From the proteins obtained using the raw NBS HMM, a high-quality protein set (E-value?<?1 × 10?20 and manual verification of an intact NBS domain) was aligned and used to construct a cassava-specific NBS HMM using hmmbuild from the HMMER v3 suite. This new cassava-specific HMM was used, and all proteins with an E-value lower than 0.01 were selected. NBS-LRR genes were further filtered based on manual curation and functional annotation against both the closest homolog from Arabidopsis and the UNIREF100 sequence database. Most of the proteins that were removed had at least a partial kinase domain, but no relationship to NBS-LRR genes; this result was expected because the NBS domain has smaller kinase subdomains

這副圖就是對(duì)應(yīng)了該文章的基因家族鑒定思路，首先在全基因組的范圍內(nèi)使用hmmersearch和NBS-ARC基因家族的隱馬可夫模型進(jìn)行基因家族的進(jìn)行初步搜索，接著把質(zhì)量比較高的基因家族候選基因篩選出來(lái)E-value?<?1 × 10?20， 然后使用clustalw2對(duì)高質(zhì)量的序列進(jìn)行多序列比對(duì)，多序列比對(duì)后，對(duì)這些置信的序列進(jìn)行隱馬可夫模型的構(gòu)建（使用hmmbuild），最后使用該新建的隱馬可夫模型，進(jìn)一步篩選完整的NSB基因家族序列（需再次過(guò)濾，找到基因家族的成員數(shù)量一般比第一步初步篩選的多）。

把該流程用到我測(cè)試數(shù)據(jù)中：

##目標(biāo)基因家族搜索
hmmsearch --cut_tc --domtblout NBS-ABC.out NBS-ARC.hmm Arabidopsis_thaliana.TAIR10.pep.all.fa

##簡(jiǎn)單看看運(yùn)算的初步輸出結(jié)果
head NBS-ABC.out
##output 
#                                                                            --- full sequence --- -------------- this domain -------------   hmm coord   ali coord   env coord
# target name        accession   tlen query name           accession   qlen   E-value  score  bias   #  of  c-Evalue  i-Evalue  score  bias  from    to  from    to  from    to  acc description of target
#
AT1G61180.1          -            889 NB-ARC               PF00931.22   252   1.4e-90  304.3   0.6   1   1   2.2e-92   2.5e-90  303.5   0.6     1   251   156   397   156   398 0.99 pep chromosome:TAIR10:1:22551271:22554684:1 gene:AT1G61180 transcript:AT1G61180.1 gene_biotype:protein_coding transcript_biotype:protein_coding description:LRR and NB-ARC domains-containing disease resistance protein [Source:UniProtKB/TrEMBL;Acc:Q2V4G0]
AT1G61180.2          -            899 NB-ARC               PF00931.22   252   1.5e-90  304.2   0.6   1   1   2.2e-92   2.5e-90  303.5   0.6     1   251   156   397   156   398 0.99 pep chromosome:TAIR10:1:22551271:22554684:1 gene:AT1G61180 transcript:AT1G61180.2 gene_biotype:protein_coding transcript_biotype:protein_coding description:LRR and NB-ARC domains-containing disease resistance protein [Source:UniProtKB/TrEMBL;Acc:Q2V4G0]



###對(duì)其e-value進(jìn)行篩選，篩選出高質(zhì)量的NBS-LRR蛋白質(zhì)序列。
grep -v "#" NBS-ABC.out|awk '($7 + 0) < 1E-20'|cut -f1 -d  " "|sort -u > NBS-ARC_qua_id.txt

~/biosoft/seqtk/seqtk subseq Arabidopsis_thaliana.TAIR10.pep.all.fa NBS-ARC_qua_id.txt >NBS-ARC_qua.fa

對(duì)篩選出來(lái)的序列,使用clustalw2進(jìn)行多序列的比較

clustalw2
 **************************************************************
 ******** CLUSTAL 2.1 Multiple Sequence Alignments  ********
 **************************************************************

     1. Sequence Input From Disc
     2. Multiple Alignments
     3. Profile / Structure Alignments
     4. Phylogenetic trees

     S. Execute a system command
     H. HELP
     X. EXIT (leave program)

Your choice: 1
Sequences should all be in 1 file.
7 formats accepted: 
NBRF/PIR, EMBL/SwissProt, Pearson (Fasta), GDE, Clustal, GCG/MSF,                  RSF.
Enter the name of the sequence file : NBS-ARC_qua.fa
 **************************************************************
 ******** CLUSTAL 2.1 Multiple Sequence Alignments  ********
 **************************************************************

     1. Sequence Input From Disc
     2. Multiple Alignments
     3. Profile / Structure Alignments
     4. Phylogenetic trees

     S. Execute a system command
     H. HELP
     X. EXIT (leave program)

Your choice: 1

Sequences should all be in 1 file.

7 formats accepted: 
NBRF/PIR, EMBL/SwissProt, Pearson (Fasta), GDE, Clustal, GCG/MSF,                  RSF.

Enter the name of the sequence file : new_NBS.aln

對(duì)這些置信的序列進(jìn)行隱馬可夫模型的構(gòu)建（使用hmmbuild），構(gòu)建更加準(zhǔn)確地盡可能預(yù)測(cè)所有的基因家族成員。

hmmbuild NBS-ARC.second.out  NBS-ARC_qua.aln 

hmmsearch --cut_tc --domtblout NBS-ARC.second.out NBS-ARC_qua.hmm ../Arabidopsis_thaliana.TAIR10.pep.all.fa

最后對(duì)再次對(duì)這些基因進(jìn)行過(guò)濾與提取

grep -v "#" NBS-ABC.second.out|awk '($7 + 0) < 1E-03' | cut -f1 -d " "|sort -u >final.NBS.list

~/biosoft/seqtk/seqtk subseq Arabidopsis_thaliana.TAIR10.pep.all.fa final.NBS.list >final_NBS-ARC_qua.fa

BLAST-based method

除了使用隱馬可夫模型和hmmer搜索的方法外，使用同源比對(duì)blast的方法也是鑒定基因家族的其中一種方法。

首先我去了NCBI下載所有植物的存在于Ref-seq（一般認(rèn)為是比較置信的植物基因序列）中的NBS序列

makeblastdb -in ref.nbs.plant.fa -dbtype prot

cat blastp.out |awk '$3>75' |cut -f1 |sort -u > blastp_result_id.list

最后我們還可將上述兩種方法重合的gene id，找出兩種方法共有的基因家族，這樣結(jié)果就比較置信了。

comm -12 blastp_result_id.list hmm_out_id.list > common.list

~/biosoft/seqtk/seqtk subseq Arabidopsis_thaliana.TAIR10.pep.all.fa common.list >final_searh_NBS-ARC_qua.fa

最后可以通過(guò)一些網(wǎng)上的保守結(jié)構(gòu)域搜索網(wǎng)頁(yè)，進(jìn)一步對(duì)所找出的結(jié)果進(jìn)行驗(yàn)證：

比如：NCBI CD-Search tool
https://www.ncbi.nlm.nih.gov/Structure/bwrpsb/bwrpsb.cgi

Pfam的搜索：
https://pfam.xfam.org/search#tabview=tab1

又或者：InterProScan sequence search
https://www.ebi.ac.uk/interpro/search/sequence-search

這些工具都可再次驗(yàn)證所搜尋的蛋白質(zhì)序列是不是含有基因家族對(duì)應(yīng)的domain。在查看保守結(jié)構(gòu)域后，如果該區(qū)域含有NBS所對(duì)應(yīng)的保守結(jié)構(gòu)域，例如LRR區(qū)域等，該蛋白質(zhì)序列可以保留進(jìn)行后續(xù)的分析。如果在該區(qū)域沒(méi)有找到對(duì)應(yīng)的保守區(qū)域，為了分析的嚴(yán)謹(jǐn)性，需進(jìn)行進(jìn)一步的排查來(lái)確定是否要去掉該序列。這種情況一般分為兩種情況，第一種就是注釋無(wú)誤，該序列確實(shí)丟失了對(duì)應(yīng)的保守結(jié)構(gòu)域，需要去掉。第二種情況就是注釋有誤，該序列的結(jié)構(gòu)域可能沒(méi)有被完整的保留下來(lái)，這種情況應(yīng)該截取該序列的上下游重新注釋分析。

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成，瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明：文章內(nèi)容（如有圖片或視頻亦包括在內(nèi)）由作者上傳并發(fā)布，文章內(nèi)容僅代表作者本人觀點(diǎn)，簡(jiǎn)書(shū)系信息發(fā)布平臺(tái)，僅提供信息存儲(chǔ)服務(wù)。

相關(guān)閱讀更多精彩內(nèi)容

友情鏈接更多精彩內(nèi)容