更多內(nèi)容關(guān)注微信公眾號：“今日之森”? ?后獲取。
直奔主題
基因家族分析（1）——基因家族成員的確定①簡單介紹了一下分析策略和數(shù)據(jù)下載。
今天就隨便找一篇文章，通過實(shí)踐操作來快速重復(fù)出文章的主要結(jié)果。
那就以下面這篇文章為例，事先沒有看，直接進(jìn)行分析。
Genome-wide analysis of the potato Hsp20 gene family: identification, genomic organization and expression profiles in response to heat stress.DOI: 10.1186/s12864-018-4443-1
第一步：由于是對馬鈴薯Hsp20基因家族進(jìn)行鑒定和分析，所以，首先需要下載馬鈴薯基因組的相關(guān)數(shù)據(jù)，下載方法在基因家族分析（1）——基因家族成員的確定①中已經(jīng)提及。
進(jìn)入EnsemblPlants數(shù)據(jù)庫http://plants.ensembl.org/info/data/ftp/index.html搜索Solanum tuberosum，下載這幾項(xiàng)即可。
因?yàn)橹跋螺d過，就不重復(fù)下載了，就是下面的四個(gè)文件。
第二步：從Pfam網(wǎng)站https://pfam.xfam.org/下載Hsp20的隱馬爾可夫模型。
ps.感興趣的可以用記事本打開這個(gè)Hsp20.hmm看看。
第三步：即用下載好的hmm模型作為query，來搜索馬鈴薯蛋白序列數(shù)據(jù)。在這里就需要用到hmmer子程序hmmsearch。
（具體安裝這里先不寫，感興趣的可以自行解決，很簡單，稍微百度就可以，和一般的程序安裝差不多。但其實(shí)很多大佬早就把基于命令行的操作打包在便于Windows下操作的軟件中了，這里不得不提CJ大神開發(fā)的Tbtools，但是我沒有用過這個(gè)軟件里面的hmmsearch功能，但由于Tbtools的操作過于簡單，就不進(jìn)行操作了，感興趣可自行操作。另外還有SPDE，也是一位大佬一直打磨的軟件，應(yīng)該也有很多功能，但是我從來沒用過任何相關(guān)功能，感興趣的都可以去搜一下。）
3.1 先使用基于ubuntu/Linus環(huán)境，在Windows中操作，由于還要安裝ubuntu等，我也忘記怎么安裝了，所以感興趣的自己琢磨安裝即可，很簡單。
由于我安裝過，界面如下：
然后運(yùn)行以下命令行即可。
hmmsearchHSP20.hmmSolanum_tuberosum.SolTub_3.0.pep.all.fa>St_hmmsearch_hsp20.out
很簡單，輸入代碼之后回車即可。沒錯(cuò)，一行命令即可完成搜索，由于命令行過于簡單，簡單解釋一下，hmmsearch：使用hmm模型搜索；HSP20.hmm：剛才下載的模型文件；Solanum_tuberosum.SolTub_3.0.pep.all.fa：馬鈴薯蛋白序列文件；St_hmmsearch_hsp20.out為輸出結(jié)果文件，可以自行命名，以.out結(jié)尾。
運(yùn)行結(jié)束后，即可生成St_hmmsearch_hsp20.out文件。如下：
將該文件用記事本打開（由于這個(gè)文件很小，記事本可以輕松打開，但由于記事本太不友好，因此，這里推薦使用Notepad軟件（安裝很容易，搜索下載并安裝即可）打開文本文檔）。打開后如下：
這里展示了幾列信息，很好辨識，左邊三列是全序列比對結(jié)果（按score從高到低排序），紅框中的三列是hmm模型結(jié)構(gòu)域最佳匹配結(jié)果（排序同上）。
這里我們主要以紅色框中的搜索結(jié)果為準(zhǔn)，由于搜索到的序列較少，故不以e值限制（一般以10-5進(jìn)行限制），即將搜索到的結(jié)果全部保留。
第四步：文章中還通過擬南芥hsp20基因作為query進(jìn)行本地blastp馬鈴薯蛋白序列。通過查看文章附件發(fā)現(xiàn)，附件一包含擬南芥hsp20基因家族的ID信息。
下載后打開附件
紅色方框中即為我們需要的信息。得到了擬南芥hsp20基因ID，接下來肯定要獲取擬南芥hsp20基因家族成員的蛋白序列，這里就不得不提Cj博士的Tbtools軟件了。
操作之前需要說明，由于同一基因在基因組數(shù)據(jù)庫中會有多個(gè)轉(zhuǎn)錄本，例如經(jīng)常見到的基因A.1?? ?A.2? ? A.3，其實(shí)這三個(gè)不同的轉(zhuǎn)錄本是同一個(gè)基因，所以習(xí)慣上經(jīng)常取最長的一個(gè)轉(zhuǎn)錄本作為該基因的代表。都說清楚了，就請大家閱讀Cj博士的文章TBtools: An Integrative Toolkit Developed for Interactive Analyses of Big Biological Data，了解tbtools，并下載安裝該軟件。
馬鈴薯基因組似乎每個(gè)基因都只保留了一個(gè)轉(zhuǎn)錄本，所以就不需要取最長序列作為代表。直接進(jìn)行操作。
4.1打開tbtools，獲取擬南芥每個(gè)基因的最長轉(zhuǎn)錄本。（當(dāng)然，剛才下載馬鈴薯基因組信息的時(shí)候，也要下載擬南芥的基因組信息）
start后即可生成以最長轉(zhuǎn)錄本為代表的擬南芥蛋白序列數(shù)據(jù)。
接下來獲取擬南芥hsp20基因蛋白序列（At_hsp20）
start之后得到At_hsp20.fasta文件，即為獲取到的蛋白序列。
記事本打開后如下：
4.2 本地blastp，這里需要下載本地blast軟件?？雌饋硇枰惭b的軟件還蠻多的，那就這里也先不提軟件安裝了，等有機(jī)會統(tǒng)一寫一下全部軟件的安裝。
首先，輸入如下代碼，先建立馬鈴薯蛋白數(shù)據(jù)庫。
makeblastdb-inSolanum_tuberosum.SolTub_3.0.pep.all.fa-dbtypeprot-outref
出現(xiàn)這三個(gè)文件即表示建庫成功。
然后進(jìn)行本地blastp，輸入如下代碼：
blastp -query At_hsp20.fasta -dbref-evalue1e-5-num_threads4-max_target_seqs5-outblast_hsp20.txt -outfmt6
輸出文件即為比對結(jié)果。打開后如下：
第五步：通過第三步獲得了hmmsearch后得到的St_hsp20基因家族成員，通過第四步獲得了本地blastp后得到的St_hsp20基因家族成員。接下來，就需要將兩種方式獲得的基因進(jìn)行合并，去掉重復(fù)值，注意，這里僅僅是去掉重復(fù)值，并不意味著去掉重復(fù)值后就是馬鈴薯hsp20基因家族的成員，因?yàn)?，去掉重?fù)值后還需要進(jìn)一步研究每個(gè)基因所包含的結(jié)構(gòu)域，這就需要明確hsp20基因包含哪些保守domain。
好，先講兩次獲得的基因ID進(jìn)行合并刪除重復(fù)值。在Excel中操作。
通過合并刪除重復(fù)值后保留了60個(gè)基因，其實(shí)仔細(xì)看會發(fā)現(xiàn)這60個(gè)就是hmmsearch搜索的結(jié)果，包含了blastp得到的全部結(jié)果。這是因?yàn)槠鸪鯖]有限定e值得原因，沒有關(guān)系。
繼續(xù)分析。
第六步：同樣的方式獲取這60個(gè)基因的蛋白序列。
注意看，得到的基因ID后面有很多描述信息，這是我們不需要的，因此還需要借助tbtools進(jìn)行ID簡化。
簡化后就美好了很多。
但是問題來了，我們確定了60個(gè)基因。明顯比文章中的48個(gè)基因要多，那是因?yàn)闆]有進(jìn)行結(jié)構(gòu)域和分子量分析。
因此在進(jìn)行基因家族分析之前很有必要研究一下該基因家族的具體信息。
好，先進(jìn)行分子量計(jì)算，由于hsp20基因分子量在15-42kDa之間，不在區(qū)間之內(nèi)的需要排除。
進(jìn)入https://web.expasy.org/compute_pi/提交剛才簡化ID之后的蛋白序列文件即可。（但需要注意提交格式，這里每個(gè)基因一行）
（最右兩列一列的等電點(diǎn)，一列是分子量。在excel中簡單處理即可排除不在15-42kDa范圍內(nèi)的基因。排除9個(gè)，保留51個(gè)。通過smart查看結(jié)構(gòu)域發(fā)現(xiàn)其中有兩個(gè)基因不包含hsp20結(jié)構(gòu)域，故最終鑒定該基因家族有49個(gè)成員。這段分析可以忽略，由于使用的馬鈴薯基因組版本和文章中使用的不同，故在鑒定過程中會有差異，其實(shí)通過smart結(jié)構(gòu)域查看即可，無需刻意通過分子量進(jìn)行篩選。詳細(xì)的分析如果需要可自行探索，由于文章中后續(xù)涉及的套路式的分析過于簡單，這里不再進(jìn)行），如果朋友需要了解，可在后臺私信。
先寫到這里。
