基因家族生信分析

1、頁眉內(nèi)容基因家族生信分析一、什么是基因家族概念：是來源于同一個祖先，有一個基因通過基因重復(fù)而產(chǎn)生兩個或更多的拷貝而構(gòu)成的一組基因， 他們在結(jié)構(gòu)和功能上具有明顯的相似性， 編碼相似的蛋白質(zhì)產(chǎn)物。劃分：按功能劃分 ：把一些功能類似的基因聚類，形成一個家族。按照序列相似程度劃分: 一般將同源的基因放在一起認(rèn)為是一個家族。1. 常見基因家族：WRKY基因家族 ：是植物前十大蛋白質(zhì)基因家族之一，大量研究表明， WRKY 基因家族的許多成員參與調(diào)控植物的生長發(fā)育，形態(tài)建成與抗病蟲。NBS-LRR抗病基因家族 ：是植物中最大類抗病基因家族之一。MADS-BOX基因家族： 是植物體內(nèi)的重要轉(zhuǎn)錄因子，它們廣泛地

2、調(diào)控著植物的生長、發(fā)育和生殖等過程。在植物中參與花器官的發(fā)育，開花時間的調(diào)節(jié)，在果實，根，莖，葉的發(fā)育中都起著重要的作用。熱激蛋白 70 家族（HSP70）是一類在植物中高度保守的分子伴侶蛋白， 在細(xì)胞中協(xié)助蛋白質(zhì)正確折疊。二、基因家族分析流程 ：利用蛋白保守域結(jié)構(gòu)提取號在 Pfam 數(shù)據(jù)庫提取其隱馬爾科夫模型矩陣文件（ *.hmm）在數(shù)據(jù)庫（ Ensemble 、JGI、NVBI）下載你所需要的物種的基因組數(shù)據(jù)（ *.fa,*.gff ）在虛擬機(jī)中 Bio-Linux 中的 hummsearch程序，用隱馬爾科夫模型矩陣文件在蛋白序列文件中搜索含有該保守結(jié)構(gòu)域的蛋白將蛋白序列導(dǎo)入 MEGA軟

3、件構(gòu)建進(jìn)化樹（可以闡明成員之間系統(tǒng)進(jìn)化關(guān)系，從進(jìn)化關(guān)系上揭示其多樣性）利用 MEME搜索蛋白質(zhì)的保守結(jié)構(gòu)域利用 MEME搜索基因家族成員的 motif 可以揭示基因家族在物種內(nèi)的多樣化及其功能， 如果他們都含有相同的 motif 表明其功能具有相似性，如果部分家族成員含有其他不同的 motif, 很可能這些成員有其他特異功能，或者可以歸分為一個亞族繪制基因染色體位置圖從*.gff文件中抽取我們搜索到的基因位置信息，_v2.0/ 在線繪制基因染色體位置圖通過染色體位置分布，可以了解基因主要分布字哪條染色體上，及是頁眉內(nèi)容否能形成基因簇（被認(rèn)為是通過重組與錯配促進(jìn)基因交流）基因結(jié)構(gòu)分析從 gff

4、文件中抽取基因的結(jié)構(gòu)信息，繪制轉(zhuǎn)錄本結(jié)構(gòu)圖。計算串聯(lián)重復(fù)基因的Ka,Ks1. 首先 將篩 選到 的基 因 的 cds 序列 進(jìn)行 多序 列對比， 篩選 identity > 75% ，tength 大于對比的兩條序列中較長的那條的長度的 75%，將篩選到的基因分別用 clustalw 進(jìn)行比對，比對結(jié)果導(dǎo)入 KsKs_Calculster 計算 Ka，Ks、Ka/ks 比，計算核苷酸的非同義替代（ ka）與核苷酸的同義替代（ks）的平均速率。2. Ka/ks 比值 <1 表明：通過純化選擇降低了氨基酸變化的速率；比值=1 表示中性選擇；比值 >1, 表明這些基因可能已經(jīng)收到積

5、極選擇，有利于適應(yīng)性遺傳，這些受正向選擇的基因?qū)⒆鳛橐院蟮难芯恐攸c。軟件的安裝從圖片中獲得進(jìn)入NCBI-blast 官網(wǎng)復(fù)制 blast-linux 版本的鏈接在 Linux 終端1. blast的安裝#wget blast 鏈接#tar xvfz 文件名 解壓縮文件# mv 解壓縮文件 /root/local/app# mv 解壓縮文件 blast頁眉內(nèi)容# vi .bashrc#在最后一行添加 export $PATH=/root/local/app/blast/bin:$PATH 并保存退出#source .bashrc 運行#blastp -version 查看是否安裝成功。2humm

6、er的安裝#yum install -y wget/安裝 wget#wgethmmer 源碼鏈接# vi .bashrc#(在最末端添加的語句 ) PATH=$PATH:/biosoft/ hmmer-#yum install -y gcc#./configure#make#make check#make install#which hmmsearch 查看是否安裝成功。3.perl的安裝#wget源代碼鏈接# tar xvfz perl- 解壓縮#./configure#make#make install 安裝完成。3.bioperl的安裝#wget -O -| bash#perlbrew

7、install-cpanm#/root/perl5/perlbrew/bin/cpanm Bio:Perl三、具體操作 ：1.保守域結(jié)構(gòu)分析頁眉內(nèi)容下載蛋白保守結(jié)構(gòu)域文件、cds、cDNA、gff注釋文件和隱馬爾科夫矩陣模型。以擬南芥為例：下載完成后，需要將文件傳到Linux 系統(tǒng)上進(jìn)行分析：打開虛擬機(jī)輸入 ip a 將虛擬機(jī) IP 連接到 Xshell 上，在 Xshell 上進(jìn)行操作，將文件通過 xftp （同樣需要連接 IP ）傳到 Linux 系統(tǒng)上，然后進(jìn)行解壓。( 一個命令解壓多個文件 :gunzip*.gz)接下來用 hummsearch尋找含有該蛋白保守結(jié)構(gòu)域的蛋白及核酸序列安

8、裝 hummsearchyum install -y wget#wgethmmer 源碼鏈接/安裝wget# vi .bashrc#(在最末端添加的語句 ) PATH=$PATH:/biosoft/ hmmer-#yum install -y gcc#./configure#make#make check#make install#which hmmsearch 查看是否安裝成功。解壓文件移動到 APP目錄下面在 app 目錄下面新建文件夾mkdir hmmer將刪除安裝包打開文字編輯器vi /.bashrc在文字編輯器里最后一行添加以上內(nèi)容安裝好 wget如果 make check 出現(xiàn)錯誤

9、百度用以下方法解決出現(xiàn) complete 安裝完成#source /.bashrc#which hmmsearch至此 hmmer安裝完成。頁眉內(nèi)容虛擬機(jī)操作：1.導(dǎo)入下載好的文件；2.hmmsearch -cut_tc domtblout NB-ARC.txt NB-ARC.hmm Arabidopsis_thaliana.TAIR10可以用 editplus 打開 .txt 文件3. perl domain_xulie.pl 結(jié)果文件 .txt 蛋白序列文件domain.fa 1e-204.clustalw 進(jìn)行多序列比對，得到aln 文件和 dnd 文件。5.hmmbuild 擬南芥特異

10、的hmm 模型文件domain.aln6.hmmsearchcut_tcdomtbloutresult.txtnewhmm 文件 蛋白質(zhì)序列文件7.在 Excel 中，根據(jù)特定的evalue 進(jìn)行篩選，并對第一列進(jìn)行去重復(fù)，得到第一列去重復(fù)的 id，保存為 id.txt8. 用 perl 腳本根據(jù) id 提取序列Per get_fa_by_id.pl id.txt 蛋白序列 wenjain >結(jié)果輸出文件可以根據(jù)的得到的序列文件進(jìn)行后續(xù)的構(gòu)建進(jìn)化樹、motif 分析等。2.搜索基因家族成員的MOTIF2.1 需要準(zhǔn)備的文件1. 擬南芥 NBS 基因蛋白質(zhì)序列2. 蛋白保守結(jié)構(gòu)域的隱馬爾科

11、夫模型矩陣文件2.2 MOTIF 的搜索使用 meme 軟件命令：meme nbs_pep.fa -protein -oc nbs_motif -nostatus -maxsize 600000 -mor anr-nmotifs 10 -minw 6 -maxw 50搜索結(jié)果存放在nbs_motif 文件夾中。文件夾中的eps 文件可以用AI打開編輯，可以另存為png或 jpg 格式，也可打開網(wǎng)頁版，也可用tbtools軟件打開，下載motif在基因上的位置信息。3. 繪制基因在染色體上的位置圖3.1 需要準(zhǔn)備的文件1. 擬南芥 NBS 基因 id2. 擬南芥的注釋文件（ gff3 文件）頁眉

12、內(nèi)容3. 擬南芥基因組長度4.1 在線繪圖工具： MapGene2Chrom4.2 samtools faidx 擬南芥 . 可得到擬南芥 . 該文件包括各個染色體，染色體長度。4.3 對基因的 id 文件在 Excel 中進(jìn)行分列，去重復(fù)處理。4.4 使用處理過的 id 文件，對擬南芥的注釋文件進(jìn)行篩選使用 perl 腳本得到基因在染色體上的位置。命令： perlget_gene_gff.pl -in1 基因的 id 文件-in2 擬南芥 gff3 文件 -out 新文件名稱4.5 新文件存放的是基因在染色體上的位置4.6 在在線文件 MapGene2chrom 中，將基因在染色體上的位置信

13、息文件復(fù)制到，input1 框中，在 input2 中粘入 samtools得到的 fai 文件。4.繪制轉(zhuǎn)錄本的結(jié)構(gòu)圖4.1 需要準(zhǔn)備的文件1. 擬南芥 NBS 基因轉(zhuǎn)錄本 id（通過家族成員鑒定得到的蛋白2.擬南芥基因的注釋文件（gtf 文件）id 文件）3.在線繪圖工具： Gene Structure Display Server 2.0http： /4.2 具體方法1. 準(zhǔn)備 gtf 文件：輸入命令： gffreadgff3 注釋文件-T-o輸出文件（gtf文件）2.editplus 打開 gtf 文件，去除 ” transcript:”3. 使用 perl 腳本提取擬南芥轉(zhuǎn)錄本結(jié)構(gòu)信

14、息：命令：perlget_gtf.pl-in1擬南芥轉(zhuǎn)錄本 id 文件 -in2 gtf 文件 -out 輸出文件（ nbs_gtf.txt）4.通過在線繪圖工具，進(jìn)行繪圖。5.篩選出串聯(lián)重復(fù)基因5.1 準(zhǔn)備文件1.擬南芥 NBS 基因 CDS 序列串聯(lián)重復(fù)基因篩選標(biāo)準(zhǔn)【 （a）length of alignable sequence covers>75% of longer頁眉內(nèi)容gene,and(b) similarity of aligned regions >75%】參考文獻(xiàn)： Extent of gene duplication in the genomes of Dro

15、sophila, nematode, andyeast.2.由于篩選時產(chǎn)生的文件較多，因此創(chuàng)建新的目錄：mkdir 新目錄3.用 editplus 打開家族成員的id 文件，對轉(zhuǎn)錄本id 進(jìn)行處理，使一個基因只拿一個轉(zhuǎn)錄本。4.把 id 復(fù)制到 Excel，首先排序處理， 然后進(jìn)行分列， 然后以第一列刪除重復(fù)值。最后將第一列和第二列進(jìn)行合并。將處理好的id 導(dǎo)入 Linux 。5.使用 perl 腳本提取 cds 序列：命令： perlget_fa_by_id.plid文件 擬南芥cds序列文件> cds.fa6.使用 blast 軟件篩選串聯(lián)重復(fù)基因6.1建立目標(biāo)序列的數(shù)據(jù)庫： mak

16、eblastdb-in cds.fa -dbtype nucl -title cds.fa6.2進(jìn)行多序列比對： blastn -query cds.fa-db cds.fa -evalue 1e-20 -outfmt9-out result.txt6.3用 editplus 打開6.4得到 cds 序列的長度，使用 samtools工具建立索引：命令 ：samtools faidxcds.fa6.5 用 perl 腳本對result.txt 進(jìn)行篩選， perl KAKS_SHAIXUAN.pl-in1-in2result.txt-outcleanresult.txt6.6 用 editpl

17、us 打開，將內(nèi)容復(fù)制到Excel ，在id 后插入一列用公式：if（ A1>B1,A1&B1,B1&A1 ）。然后全選，以第 C 列刪除重復(fù)值。并保存到新的文件中，并導(dǎo)入到 Linux 中7. 計算串聯(lián)重復(fù)基因的 KaKs 。7.1 準(zhǔn)備文件1.串聯(lián)重復(fù)基因的CDS 序列文件7.2 將成對的串聯(lián)重復(fù)序列保存在一個文件中。、方法一：復(fù)制需要找到的序列的id，在 editplus 中按 ctrl+F 搜索，找到后復(fù)制粘貼到一個文件中。方法二：首先將成對的id 保存在同一個文件，導(dǎo)入到Linux 中，在 Linux 中，利頁眉內(nèi)容用 perl 腳本提取序列： perlget_

18、fa_by_id.pl新的 id 文件 cds.fa 文件 >id1.fa7.3 計算 KaKS1.計算之前需要使用CLUSTAW 對序列進(jìn)行比較?？色@得id1.aln。2.使用 KaKs_calculator 工具將 id1.aln 文件轉(zhuǎn)換成 id1.axt 文件命令： axtvenvertorid1.alnid1.axt3.計算 KaKs，輸入命令：KaKs_calculstor-iid.axt-oid1_kaks.txt4.如果報錯，則把兩條序列長度保持一致。依此將所有的串聯(lián)重復(fù)基因?qū)?，進(jìn)行計算。四基因家族成員的鑒定（未知隱馬爾科夫模型）1. 鑒定測略?在 NCBI 數(shù)據(jù)庫中盡量

19、多下載幾個物種的需要鑒定的蛋白保守結(jié)構(gòu)域序列，以及所要研究物種的所有蛋白序列?在虛擬機(jī)中本地建庫，并進(jìn)行blast建庫命令： makeblastdb -in 研究物種的蛋白序列文件-dbtype prot -title 庫名稱?進(jìn)行序列比對命令： blastp -query 下載的多個物種序列文件-db 庫名稱 -evalue 1e-10 -outfmt6 -out 結(jié)構(gòu)域 .blast?使用 sed 命令去除表頭和結(jié)尾得到新的new 結(jié)構(gòu)域 . blast?Awk print$1結(jié)new構(gòu)域 . Blast |less 可查看打印的結(jié)果? Awk print$1結(jié)構(gòu)new域 . Blast

20、> id.txt?cat id.txt|sort|uniq > idd.txt去重復(fù)?Perlget_fa_by_id.plidd.txt 去重復(fù) 蛋白序列 wenjain >結(jié)果輸出文件?在 Pfam 或者 NCBI 的 cdd 中搜索檢查是否有相關(guān)蛋白結(jié)構(gòu)域。之后再進(jìn)行 motif分析一些命令及軟件應(yīng)用說明（參考一些視頻資料）hmmsearch使用說明用途：利用蛋白保守結(jié)構(gòu)域的隱馬爾科夫模型搜索蛋白序列中具有該保守結(jié)構(gòu)域的蛋白用法： hmmsearch -cut_tc -domtblout result.txt*.hmm頁眉內(nèi)容說明： result.txt是輸出的結(jié)果文件

21、，*.hmm在 pfam數(shù)據(jù)庫下載的模型，hmmbuild使用說明用途：利用clustalw比對生成的aln文件構(gòu)建蛋白保守結(jié)構(gòu)域的隱馬爾科夫模型用法：hmmbuildnew.hmmdomain.aln說 明 ： new.hmm是結(jié)果文件也就是構(gòu)建的蛋白保守結(jié)構(gòu)域的隱馬爾科夫模型，domain.aln是 clustalw比對生成的aln 文件domain_xulie.pl腳本使用說明用途：提取hmmsearch搜索結(jié)果中蛋白序列中保守結(jié)構(gòu)域的序列，用于構(gòu)建新的物種特異的蛋白保守結(jié)構(gòu)域的隱馬爾科夫模型用 法 ： perldomain_xulie.pl（ 腳 本 不 在 使 用 目 錄 下 要 寫

22、 全 路 徑 ） hmmoutfiledomain.fastaE-value說明： hmmoutfile是 hmmsearch搜索結(jié)果文件domain.fasta是結(jié)果存放文件也就是蛋白序列中保守結(jié)構(gòu)域的序列，E-value是提取序列時設(shè)定的 E 值get_fa_by_id.pl 使用說明用途：通過 ID號獲取其相應(yīng)的基因或蛋白序列用法： perl perl get_fa_by_id.pl id.txtcds.fastat>id_cds.fasta說明： id.txt 是包含你的 ID 的文件， cds.fasta 是你叢數(shù)據(jù)庫中下載的包含所有cds 序列的文件， id_cds.fast

23、a是輸出文件內(nèi)容是 ID 對應(yīng)的序列samtools faidx用途：提取fasta文件信息用法： samtools faidx *.fa說明：輸入文件是 fasta 文件，自動生成輸出目錄結(jié)果的 fai 文件第一列是你輸入的 fasta 文件的 ID 第二列是其序列長度Gffread使用說明用途：將基因組注釋文件 gff3 轉(zhuǎn)化成基因的注釋文件 gtf 用法： gffread my.gff3 -T -o my.gtf說明： my.gff3是輸入文件 基因組注釋文件，my.gtf是輸出文件是基因的注釋文件Get_gene_gff.pl用途：想要繪制基因的染色體位置圖必須要拿到基因在染色體上的具

24、體信息，該腳本就是從總的 gff 文件抽去你需要的基因的信息如：所在染色體，起始終止位置等信息用法：PerlGet_gene_gff.pl-in1 gene_id.txt-in2my.gff3 -out gene_location.txt說明：gene_id.txt是第一個輸入文件基因的ID文件， my.gff3是第二個輸入文件是物種基因組所有蛋白序列， gene_location.txt是結(jié)果輸出文件頁眉內(nèi)容Get_gtf.pl腳本使用說明用途：從基因注釋文件gtf 文件中提取轉(zhuǎn)錄本的結(jié)構(gòu)信息用法： perl get_gtf.pl-in1 id.txt-in2 gene.gtf-outstr

25、ucture.txt說明： id.txt是第一個輸入文件是id 文件 ， gene.gtf是第二個輸入文件是基因注釋文件gtf 文件structure.txt是結(jié)果輸出文件存放著轉(zhuǎn)錄本的結(jié)構(gòu)信息KaKs_shaixuan.pl使用說明用途：在多序列比對結(jié)果文件中篩選identity大于 75% ，比對上的序列長度大于對比的兩條序列中最長序列的長度的75%用法： perl KaKs_shaixuan.pl -in1 cds.fai -in2 result.txt -out shaixuan.txt說明： cds.fai是 samtools faidx對 cds 的序列文件fasta作用，生成的文件， resu