基因家族生信分析

1、基因家族生信分析一、 什么是基因家族概念：是來源于同一個祖先，有一個基因通過基因重復(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)錄因子，它們廣泛地調(diào)控著植物的生長、發(fā)育和生

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

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

4、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)收到積極選擇，有利于適應(yīng)性遺傳，這些受正向選擇的基因?qū)⒆鳛橐院蟮难芯恐攸c。 軟件的安裝從圖片中獲得進(jìn)入NCBI-blast官

5、網(wǎng)復(fù)制blast-linux版本的鏈接在Linux終端1. blast的安裝#wget blast鏈接#tar xvfz 文件名 解壓縮文件# mv 解壓縮文件 /root/local/app# mv 解壓縮文件 blast# vi .bashrc#在最后一行添加export $PATH=/root/local/app/blast/bin:$PATH 并保存退出#source .bashrc 運行 #blastp -version 查看是否安裝成功。2hummer的安裝#yum install -y wget /安裝wget#wget hmmer源碼鏈接 #tar -zxvf hmmer-3.

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

7、w install-cpanm#/root/perl5/perlbrew/bin/cpanm Bio:Perl三、 具體操作：1.保守域結(jié)構(gòu)分析下載蛋白保守結(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)域的蛋白及核酸序列安裝hummsearchyum ins

8、tall -y wget /安裝wget#wget hmmer源碼鏈接 #tar -zxvf hmmer-3.2.1# vi .bashrc#(在最末端添加的語句) PATH=$PATH:/biosoft/ /binaries#yum install -y gcc #./configure#make #make check#make install#which hmmsearch 查看是否安裝成功。解壓文件移動到APP目錄下面在app目錄下面新建文件夾 mkdir hmmer將hmmer-3.2.1移到hmmer目錄下面 mmove -v c:/hmmer-3.2.1 c:/hmmer刪除安裝

9、包打開文字編輯器vi /.bashrc 在文字編輯器里最后一行添加以上內(nèi)容安裝好wget如果make check出現(xiàn)錯誤百度用以下方法解決出現(xiàn)complete安裝完成#source /.bashrc #which hmmsearch至此hmmer安裝完成。虛擬機(jī)操作：1.導(dǎo)入下載好的文件；2.hmmsearch -cut_tc domtblout NB-ARC.txt NB-ARC.hmm Arabidopsis_thaliana.TAIR10 pep.all.fa 可以用editplus打開.txt文件3. perl domain_xulie.pl 結(jié)果文件.txt 蛋白序列文件 domai

10、n.fa 1e-204.clustalw進(jìn)行多序列比對，得到aln文件和dnd文件。5.hmmbuild 擬南芥特異的hmm模型文件 domain.aln6.hmmsearch cut_tc domtblout result.txt newhmm文件 蛋白質(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.搜索基因家族成員

11、的MOTIF2.1 需要準(zhǔn)備的文件 1. 擬南芥NBS基因蛋白質(zhì)序列 2. 蛋白保守結(jié)構(gòu)域的隱馬爾科夫模型矩陣文件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 需

12、要準(zhǔn)備的文件1. 擬南芥NBS基因id2. 擬南芥的注釋文件（gff3文件）3. 擬南芥基因組長度4.1在線繪圖工具：MapGene2Chrom4.2 samtools faidx 擬南芥.dna.toplevel.fa 可得到擬南芥.文件 該文件包括 各個染色體，染色體長度。4.3 對基因的id文件在Excel中進(jìn)行分列，去重復(fù)處理。4.4 使用處理過的id文件，對擬南芥的注釋文件進(jìn)行篩選使用perl腳本得到基因在染色體上的位置。命令： perl get_gene_gff.pl -in1 基因的id文件 -in2 擬南芥gff3文件 -out 新文件名稱 4.5 新文件存放的是基因在染色體上

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

14、ript:”3.使用perl腳本提取擬南芥轉(zhuǎn)錄本結(jié)構(gòu)信息：命令： perl get_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 gene,and(b) similarity of aligned regions >75%】參考文獻(xiàn)：Extent of gene duplication in t

15、he genomes of Drosophila, nematode, and yeast.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序列：命令： perl get_fa_by_id.pl id文件 擬南芥cds序列文件 > cds.fa6.使用blast軟件篩選串聯(lián)重復(fù)基因6.1建立目標(biāo)序列的數(shù)據(jù)庫：

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

17、t.txt6.6 用editplus打開，將內(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中，利用perl腳本提取序列：perl get_fa_by_id.p

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

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

20、rt|uniq > idd.txt去重復(fù)Perl get_fa_by_id.pl idd.txt去重復(fù) 蛋白序列wenjain >結(jié)果輸出文件在Pfam或者NCBI的cdd中搜索檢查是否有相關(guān)蛋白結(jié)構(gòu)域。之后再進(jìn)行motif分析一些命令及軟件應(yīng)用說明（參考一些視頻資料）hmmsearch使用說明用途：利用蛋白保守結(jié)構(gòu)域的隱馬爾科夫模型搜索蛋白序列中具有該保守結(jié)構(gòu)域的蛋白用法：說明：result.txt是輸出的結(jié)果文件，*.hmm在pfam數(shù)據(jù)庫下載的模型， pep.all.fa物種基因組所有蛋白序列hmmbuild使用說明用途：利用clustalw比對生成的aln文件構(gòu)建蛋白保守結(jié)

21、構(gòu)域的隱馬爾科夫模型用法：hmmbuild new.hmm domain.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)域的隱馬爾科夫模型用法：perl domain_xulie.pl（腳本不在使用目錄下要寫全路徑） hmmoutfile pep.all.fa domain.fasta E-value說明：hmmoutfile 是hmmsearch搜索結(jié)果文件

22、 pep.all.fa物種基因組所有蛋白序列 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.txt cds.fastat >id_cds.fasta 說明：id.txt是包含你的ID的文件，cds.fasta是你叢數(shù)據(jù)庫中下載的包含所有cds序列的文件，id_cds.fasta是輸出文件內(nèi)容是ID對應(yīng)的序列samtools faidx 用途：提取fasta文件信息用法：samtool

23、s faidx *.fa說明：輸入文件是fasta文件，自動生成輸出目錄*.fa.fai,結(jié)果的fai文件第一列是你輸入的fasta文件的ID第二列是其序列長度Gffread 使用說明用途：將基因組注釋文件gff3轉(zhuǎn)化成基因的注釋文件gtf用法：gffread my.gff3 -T -o my.gtf說明：my.gff3是輸入文件基因組注釋文件，my.gtf是輸出文件是基因的注釋文件Get_gene_gff.pl用途：想要繪制基因的染色體位置圖必須要拿到基因在染色體上的具體信息，該腳本就是從總的gff文件抽去你需要的基因的信息如：所在染色體，起始終止位置等信息用法：Perl Get_gene_

24、gff.pl -in1 gene_id.txt -in2 my.gff3 -out gene_location.txt說明：gene_id.txt 是第一個輸入文件基因的ID文件，my.gff3是第二個輸入文件是物種基因組所有蛋白序列 ，gene_location.txt 是結(jié)果輸出文件Get_gtf.pl腳本使用說明用途：從基因注釋文件gtf文件中提取轉(zhuǎn)錄本的結(jié)構(gòu)信息用法：perl get_gtf.pl -in1 id.txt -in2 gene.gtf -out structure.txt說明：id.txt 是第一個輸入文件是id文件 ，gene.gtf是第二個輸入文件是基因注釋文件gtf

25、文件 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作用，生成的文件，result.txt 是拿cds進(jìn)行多序列比對得到的結(jié)果，shaixuan.txt 是篩選后的結(jié)果存放的文本aln文件轉(zhuǎn)化 axt文件命令：/home/manag