目錄

• MCScanX下載及安裝
• 數(shù)據(jù)準(zhǔn)備
• MCScanX運(yùn)行及結(jié)果
• Ka & Ks 計(jì)算
• 下游畫圖腳本
• 最后

MCScanX有兩個(gè)主要的功能, 一是, 方便使用者發(fā)現(xiàn)共線性(collinearity)和同線性(synteny)關(guān)系并且可以從共線性區(qū)塊中看到清晰多重比對(duì); 二是, 通過其中眾多的輔助腳本, 更形象的分析同線性和共線性的數(shù)據(jù)
下面就從零開始共線性的分析吧~
注: 稍微解釋一下synteny and collinearity的關(guān)系, 假設(shè)有A B C三個(gè)基因, 在species1中的排列為ACB, 而在species2中排列為CAB, 則這兩個(gè)區(qū)段稱為synteny; 如果在species3中排列為ABC, 在species4中排列也為ABC, 則這兩個(gè)區(qū)段稱為collinearity

MCScanX下載及安裝

$wget http://chibba.pgml.uga.edu/mcscan2/MCScanX.zip
$unzip MCScanX.zip
$cd MCScanX
$make

這部分我是在華為云服務(wù)器上完成的, 因?yàn)楹芏鄸|西沒裝, 所以一直報(bào)錯(cuò), 報(bào)錯(cuò)內(nèi)容及解決辦法記錄如下:

• 報(bào)錯(cuò)1

make: g++: Command not found

• 解決辦法1

$conda install gcc=4.8.5

• 報(bào)錯(cuò)2

g++ struct.cc mcscan.cc read_data.cc out_utils.cc dagchainer.cc msa.cc permutation.cc -o MCScanX
msa.cc: In function ‘void msa_main(const char*)’:
msa.cc:289:22: error: ‘chdir’ was not declared in this scope
     if (chdir(html_fn)<0)
                      ^
makefile:2: recipe for target 'mcscanx' failed
make: *** [mcscanx] Error 1

• 解決辦法2
  錯(cuò)誤的原因是MCScanX 不支持64位系統(tǒng), 如果要在 64位上運(yùn)行, 需要在msa.h, dissect_multiple_alignment.h, and detect_collinear_tandem_arrays.h第一行加上

#include <unistd.h>

• 報(bào)錯(cuò)3

can't find javac

• 解決辦法3
  一般正常服務(wù)器哪有不配java的, 我也是醉了

$yum -y install  java-1.8.0-openjdk-devel.x86_64

數(shù)據(jù)準(zhǔn)備

安裝完成后, 軟件包包括三個(gè)主程序MCScanX, MCScanX_h, duplicate_gene_classifier, 和13個(gè)用于形象化展示的下游分析軟件
本次分析只用到了MCScanX和幾個(gè)下游腳本

BLASTP

• 建庫

$makeblastdb -in pde_ptr_pep.fasta -dbtype prot -parse_seqids -out pep.db/pde_ptr

• 比對(duì)
  以下是 qsub 腳本內(nèi)容, 比對(duì)之前要檢查一下序列名稱是否跟待會(huì)要用的 gff 文件序列名稱是否一致, 不一致的話可以在 vim 中通過:%s/XXX/xxx/g修改

#PBS -N pdeptr_blastp
#PBS -l nodes=1:ppn=62
#PBS -l walltime=1200:00:00
#PBS -q long
cd /home/Dai_XG/Analysis/Pdel_genome/05.evolution_analysis/05.colinearity/out_01
blastp -query \pde_ptr_pep.fasta -out pde_ptr_pep.fasta.blastpout -db pep.db/pdeptr -outfmt 6 -evalue 1e-5 -num_threads 62

特定GFF生成

需要的格式只有四列信息, 分別是染色體ID, 基因ID, 起始和終止, 以下命令行是對(duì)基因注釋文件操作的
如果只想做特定序列的結(jié)果, 在這一步中篩出來就行, 不用動(dòng)上邊的 blast 結(jié)果

#具體文件具體分析
$grep '\smRNA\s' EVM.final.gene.gff3 | awk '{print $1"\t"$4"\t"$5"\t"$9}' |awk -F 'ID=' '{print $1$2}'|awk -F 'Parent=' '{print $1}'|awk '{print $1"\t"$4"\t"$2"\t"$3}' >pde.gff 
$grep '\smRNA\s' Ptrichocarpa_444_v3.1.gene.gff3 | awk '{print $1"\t"$4"\t"$5"\t"$9}'|awk -F '.v3.1;' '{print $1}'|awk -F 'ID=' '{print $1$2}'|awk '{print $1"\t"$4"\t"$2"\t"$3}' > ptr.gff

MCScanX運(yùn)行及結(jié)果

運(yùn)行命令行

$nohup MCScanX pde_ptr > out.log 2>&1 &

結(jié)果

運(yùn)行完成后會(huì)生成三個(gè)結(jié)果文件

• pde_ptr.collinearity
• pde_ptr.html(目錄)
• pde_ptr.tandem

pde_ptr.collinearity

就不粘自己電腦上的圖了, 偷張別人的圖說一下

• 第一部分為執(zhí)行參數(shù)
• 第二部分兩行給出了共線性基因的數(shù)目和比率, 及所有基因的數(shù)目
• 第三部分全部為 block, 也就是比對(duì)出來的共線性區(qū)塊
  1.第一列為 block 編號(hào)
  2.第二列為 block 內(nèi)部, 比對(duì)基因的編號(hào)
  3.第三列為&前邊染色體的基因
  4.第四列為&后邊染色體的基因
  5.第五列 e_value 值, 跟 blast 一樣
  
  XX.collinearity

pde_ptr.html

目錄中包含所有 block 的.html 文件

• 第一列為與參考片段相應(yīng)基因匹配上的基因數(shù)目
• 第二列參考片段, 基因按順序排列

• 第三列匹配上的片段

  
  
  XX.html

pde_ptr.tandem

串聯(lián)重復(fù)的基因, 也不知道有啥用

Ka & Ks 計(jì)算

上邊說了, MCScanX有很多下游分析腳本, 我這幾篇文章都是在記錄之前做的進(jìn)化分析流程, 其中流程中做過KaKs的統(tǒng)計(jì)計(jì)算, 下面就著重介紹一下, 關(guān)于其他下游腳本, 將在下一章節(jié)介紹

定義

• Ka/Ks或者dN/dS表示異義突變率(Ka)和同義突變率(Ks)之間的比例
• 計(jì)算
  Ks = 同義突變SNP數(shù)/同義位點(diǎn)數(shù)
  Ka = 異義突變SNP數(shù)/異義位點(diǎn)數(shù)
  • 關(guān)于同/異義位點(diǎn)數(shù)的理解
    以TTT密碼子為例, 單(雙/三)堿基發(fā)生變異, 會(huì)有9(15*2/63)種變異位點(diǎn)數(shù), 通過比照密碼子表發(fā)現(xiàn), 1(未統(tǒng)計(jì))種為同義, 8(未統(tǒng)計(jì))種為異義, 以此類推, 序列的同/異義位點(diǎn)數(shù)要依次統(tǒng)計(jì)每個(gè)密碼子變異后的情況, 也就是說分母為整條序列所有單堿基變異的可能之和
  • 關(guān)于同/異義突變SNP數(shù)的理解
    以密碼子TTT為例, 突變?yōu)镚AC(以上所說的三堿基突變)自然是有三個(gè)SNP突變數(shù), 但此突變數(shù)非彼突變數(shù), 因?yàn)槲覀儾恢雷儺愡^程究竟是怎么樣的, 所以要遍歷整個(gè)變異流程, 然后統(tǒng)計(jì)整個(gè)變異流程中的同/異義的數(shù)目, 最后取平均值, TTT變?yōu)镚AC有6中流程, 過程中有14中異義突變SNP數(shù), 4種同義突變SNP數(shù), 故同義突變SNP數(shù) = 4/6, 異義突變SNP數(shù) = 14/6
• Ka>>Ks或者Ka/Ks >> 1, 基因受正選擇(positive selection); Ka＝Ks或者Ka/Ks ＝1, 基因中性進(jìn)化; Ka<<Ks或者Ka/Ks << 1, 基因受純化選擇(purify selection)
• 可以參考用以判斷選擇壓力的Ka/Ks的計(jì)算

數(shù)據(jù)準(zhǔn)備

• CDS序列
  合并物種的 cds 序列文件, 序列名稱與 .collinearity中一定要一致, 不然會(huì)發(fā)現(xiàn)結(jié)果中全都是-2 -2 -2 -2 ...
• 上一步生成的 .collinearity

軟件配置

這一步需要下游腳本add_ka_and_ks_to_synteny.pl, 如果你也是普通用戶, 那么就很操蛋了, 要手動(dòng)下載安裝很多個(gè)perl model, 還要裝ClustalX, 而且必須是2.0版本以前的

• 搞定ClustalW

$wget http://www.clustal.org/download/1.X/ftp-igbmc.u-strasbg.fr/pub/ClustalW/clustalw1.83.linux.tar.gz
$tar zxvf clustalw1.83.linux.tar.gz
$vi ~/.bash_profile

export PATH=/home/Dai_XG/software/clustalw1.83.linux:$PATH

• Bio-Tools-Run-Alignment-Clustalw安裝
  報(bào)錯(cuò)提示沒有記錄, 反正意思就是你沒有這個(gè)庫, 我運(yùn)行不了

$wget https://cpan.metacpan.org/authors/id/C/CD/CDRAUG/Bio-Tools-Run-Alignment-Clustalw-1.7.4.tar.gz
$tar zxvf Bio-Tools-Run-Alignment-Clustalw-1.7.4.tar.gz
$cd Bio-Tools-Run-Alignment-Clustalw
$perl Makefile.PL PREFIX=/home/Dai_XG/software/prefix/perl_local_lib
$make
$make test
$make install

• Devel-CheckBin安裝

$wget https://cpan.metacpan.org/authors/id/T/TO/TOKUHIROM/Devel-CheckBin-0.04.tar.gz
$tar zxvf Devel-CheckBin-0.04.tar.gz
$cd Devel-CheckBin-0.04
$perl Makefile.PL PREFIX=/home/Dai_XG/software/prefix/perl_local_lib
$make
$make test
$make install

• BioPerl-Run安裝

$wget https://cpan.metacpan.org/authors/id/C/CJ/CJFIELDS/BioPerl-Run-1.007002.tar.gz
$tar zxvf BioPerl-Run-1.007002.tar.gz
$cd BioPerl-Run-1.007002

這是Bioperl, 手動(dòng)形式跟上邊的不同, 我是在沒找到普通用戶的管理方法, 但發(fā)現(xiàn)目錄下有一個(gè)lib 目錄, 把它寫入環(huán)境竟然也能用, 太好了~~

• 寫入環(huán)境

$vi ~/.bash_profile

export PERL5LIB=/home/Dai_XG/software/prefix/perl_local_lib/lib/site_perl/5.28.0/:/home/Dai_XG/software/BioPerl-Run-1.007002/lib/:$PERL5LIB
export PATH=/home/Dai_XG/software/clustalw1.83.linux:$PATH

運(yùn)行及結(jié)果

終于能用了, 太煩了

• 運(yùn)行

  $nohup perl  ~/software/MCScanX/downstream_analyses/add_ka_and_ks_to_collinearity.pl -i ../../../out_03_1/chr_result/pde_ptr.collinearity -d pde_ptr.fasta -o pde_ptr.kaks > out.log 2>&1&
  

• 結(jié)果
  在.collinearity 結(jié)果基礎(chǔ)上多了兩列, 分別是Ka&Ks

自編腳本處理

兩個(gè)腳本太長了, 備份里找吧, 名稱為extract_peak.kaks.pl, statistics_block.pl

• extract_peak.kaks結(jié)果

Chr&LachesisXXgroup
-19990  2318
0       1052
220     819
200     805
...
Chr&Chr
-19990  810
2320    73
2280    73
...
LachesisXXgroup&LachesisXXgroup
-19990  2430
2290    76
2160    67

• statistics_block結(jié)果1

## BLOCK 1      Chr01   2133    LachesisXXgroup0        1624    Ks
0       Chr01   Potri.001G064150.1      5052138 5053352 -2
1       Chr01   Potri.001G064301.1      5072110 5072703 -2
2       Chr01   Potri.001G065600.1      5191942 5199150 -2
...
0       LachesisXXgroup0        EVM0035433.1    4924636 4927232 -2
1       LachesisXXgroup0        EVM0024093.1    5012566 5014702 -2
2       LachesisXXgroup0        EVM0031968.1    5107185 5112594 -2
...
## BLOCK 2      Chr01   2118    LachesisXXgroup0        1564    Ks
0       Chr01   Potri.001G290600.3      29621454        29622702        -2
1       Chr01   Potri.001G290800.3      29637348        29639177        -2
2       Chr01   Potri.001G291866.1      29721339        29722599        -2
...

• statistics_block結(jié)果2

BLOCK   species blocklength     colinearity genes       block genes with no colinearity
BLOCK1  Chr01   20138601        1344    789     LachesisXXgroup0        20266103        1344    280
BLOCK2  Chr01   20866914        1135    983     LachesisXXgroup0        21403547        1135    429
BLOCK3  Chr01   5024695 434     330     LachesisXXgroup0        5147359 434     157
BLOCK4  Chr01   3688015 309     201     LachesisXXgroup0        4407659 309     86
BLOCK5  Chr01   1523936 21      91      LachesisXXgroup0        1523936 21      62
BLOCK6  Chr01   226943  10      25      LachesisXXgroup0        404525  10      36
BLOCK7  Chr01   308614  9       25      LachesisXXgroup0        308614  9       8
BLOCK8  Chr01   929669  10      71      LachesisXXgroup0        1682958 10      26
BLOCK9  Chr01   1606860 10      104     LachesisXXgroup0        2277577 10      123
BLOCK10 Chr01   38729   6       2       LachesisXXgroup0        143832  6       0
BLOCK11 Chr01   385605  7       44      LachesisXXgroup0        6302890 7       55
BLOCK12 Chr01   80612   6       6       LachesisXXgroup0        42799884        6       30
...

下游畫圖腳本

這里盡可能多的吧腳本都跑一遍吧, 實(shí)在沒法跑的以后補(bǔ)充
每個(gè)腳本運(yùn)行前需要.java 及.class 同時(shí)存在才可以

bar_plotter.java

• 配置文件

2000    //dimension (in pixels) of x axis
2000    //dimension (in pixels) of y axis
Chr01,Chr02,Chr03,Chr04,Chr05,Chr06,Chr07,Chr08,Chr09,Chr10,Chr11,Chr12,Chr13,Chr14,Chr15,Chr16,Chr17,Chr18,Chr19       //reference chromosomes
Lachesis_group0,Lachesis_group1,Lachesis_group2,Lachesis_group3,Lachesis_group4,Lachesis_group5,Lachesis_group6,Lachesis_group7,Lachesis_group8,Lachesis_group9,Lachesis_group10,Lachesis_group11,Lachesis_group12,Lachesis_group13,Lachesis_group14,Lachesis_group15,Lachesis_group16,Lachesis_group17,Lachesis_group18        //target chromosomes

• 運(yùn)行

$java dot_plotter -g gff_file -s collinearity_file -c control_file -o output_PNG_file

• 輸出

  
  
  圖不對(duì)題1

circle_plotter.java

• 配置文件

800     //plot width and height (in pixels)
Chr01,Chr02,Chr03,Chr04,Chr05,Chr06,Chr07,Chr08,Chr09,Chr10,Chr11,Chr12,Chr13,Chr14,Chr15,Chr16,Chr17,Chr18,Chr19,Lachesis_group0,Lachesis_group1,Lachesis_group2,Lachesis_group3,Lachesis_group4,Lachesis_group5,Lachesis_group6,Lachesis_group7,Lachesis_group8,Lachesis_group9,Lachesis_group10,Lachesis_group11,Lachesis_group12,Lachesis_group13,Lachesis_group14,Lachesis_group15,Lachesis_group16,Lachesis_group17,Lachesis_group18  //chromosomes in the circle

• 運(yùn)行

$java circle_plotter -g gff_file -s collinearity_file -c control_file -o output_PNG_file

• 輸出

  
  
  圖不對(duì)題2

dot_plotter.java

• 配置文件

800     //dimension (in pixels) of x axis
800     //dimension (in pixels) of y axis
sb1,sb2,sb3,sb4,sb5,sb6,sb7,sb8,sb9,sb10        //chromosomes in x axis
os1,os2,os3,os4,os5,os6,os7,os8,os9,os10,os11,os12      //chromosomes in y axis

• 運(yùn)行

$java dot_plotter -g gff_file -s collinearity_file -c control_file -o output_PNG_file

• 輸出

  
  
  圖不對(duì)題3

dual_synteny_plotter.java

• 配置文件

6000    //plot width (in pixels)
20000   //plot height (in pixels)
Chr01,Chr02,Chr03,Chr04,Chr05,Chr06,Chr07,Chr08,Chr09,Chr10,Chr11,Chr12,Chr13,Chr14,Chr15,Chr16,Chr17,Chr18,Chr19,      //chromosomes in the left column
Lachesis_group0,Lachesis_group1,Lachesis_group2,Lachesis_group3,Lachesis_group4,Lachesis_group5,Lachesis_group6,Lachesis_group7,Lachesis_group8,Lachesis_group9,Lachesis_group10,Lachesis_group11,Lachesis_group12,Lachesis_group13,Lachesis_group14,Lachesis_group15,Lachesis_group16,Lachesis_group17,Lachesis_group18,       //chromosomes in the right column

• 運(yùn)行

$java dual_synteny_plotter -g gff_file -s collinearity_file -c control_file -o output_PNG_file

• 輸出

  
  
  圖不對(duì)題4

最后

• 還有幾篇關(guān)于進(jìn)化的筆記隱藏了, 歡迎小伙伴們一起來交流學(xué)習(xí)

• 終于看到頭了

  
  
  For What?