多倍化以及后續(xù)的基因丟失和二倍化現(xiàn)象存在于大部分的物種中, 是物種進化的重要動力。如果一個物種在演化過程中發(fā)生過多倍化表悬，那么在基因組上就會存在一些共線性區(qū)域（即兩個區(qū)域之間的基因是旁系同源基因校读，其基因的排布順序基本一致）沼侣。

例如擬南芥經(jīng)歷了3次古多倍化，包括2次二倍化歉秫，一次3倍化 (Tang. et.al 2008 Science).

..For example, Arabidopsis thaliana (thale cress) has undergone three paleo-polyploidies, including two doublings (5) and one tripling (12), resulting in ~12 copies of its ancestral chromosome set in a ~160-Mb genome...

當然之前只是記住了結(jié)論蛾洛，現(xiàn)在我想的是，如何復現(xiàn)出這個分析結(jié)果呢雁芙？

前期準備

首先轧膘，我們需要使用conda安裝wgdi。 我一般會新建一個環(huán)境避免新裝軟件和其他軟件產(chǎn)生沖突

# 安裝軟件
conda create -c bioconda -c conda-forge -n wgdi wgdi
# 啟動環(huán)境
conda activate wgdi 

之后兔甘，我們從ENSEMBL上下載基因組谎碍，CDS, 蛋白和GFF文件

#Athaliana
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/dna/Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
gunzip Arabidopsis_thaliana.TAIR10.dna.toplevel.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
gunzip Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/pep/Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
gunzip Arabidopsis_thaliana.TAIR10.pep.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.44.gff3.gz
gunzip Arabidopsis_thaliana.TAIR10.44.gff3.gz

注意: cDNA和CDS并不相同，CDS只包括起始密碼子到終止密碼子之間的序列洞焙，而cDNA還會包括UTR.

數(shù)據(jù)預處理

數(shù)據(jù)的預處理是用時最長的步驟蟆淀，因為軟件要求的輸入格式并非是你手頭所擁有的格式，你通常都需要進行一些轉(zhuǎn)換才能得到它所需的形式澡匪。

wgdi需要三種信息熔任，分別是BLAST, 基因的位置信息和染色體長度信息，要求格式如下

1. BLAST 的 -outfmt 6輸出的文件

2. 基因的位置信息: 以tab分隔唁情，分別為chr疑苔，id，start甸鸟，end惦费，strand赛惩，order，old_id趁餐。（并非真正意義上的GFF格式）

3. 染色體長度信息和染色體上的基因個數(shù)喷兼，格式為 chr, length, gene number

同時，對于每個基因我們只需要一個轉(zhuǎn)錄本后雷，我通常使用最長的轉(zhuǎn)錄本作為該基因的代表季惯。之前我寫過一個腳本 (get_the_longest_transcripts.py) 提取每個基因的最長轉(zhuǎn)錄本，我在此處上寫了一個新的腳本用來根據(jù)參考基因組和注釋的GFF文件生成wgdi的兩個輸入文件臀突，腳本地址為https://github.com/xuzhougeng/myscripts/blob/master/comparative/generate_conf.py

pytho generate_conf.py -p ath Arabidopsis_thaliana.TAIR10.dna.toplevel.fa Arabidopsis_thaliana.TAIR10.44.gff3

輸出文件是ath.gff 和 ath.len.

由于ENSEMBL上GFF的ID編號和pep.fa 和cds.fa 不太一致勉抓，簡單的說就是，編號之前中有一個 "gene:" 和 "transcript:" . 如下

$ head ath.gff 
1  gene:AT1G01010  3631  5899  +  1  transcript:AT1G01010.1
1  gene:AT1G01020  6788  9130  -  2  transcript:AT1G01020.1
1  gene:AT1G01030  11649  13714  -  3  transcript:AT1G01030.1
1  gene:AT1G01040  23416  31120  +  4  transcript:AT1G01040.2
1  gene:AT1G01050  31170  33171  -  5  transcript:AT1G01050.1
1  gene:AT1G01060  33379  37757  -  6  transcript:AT1G01060.3
1  gene:AT1G01070  38444  41017  -  7  transcript:AT1G01070.1
1  gene:AT1G01080  45296  47019  -  8  transcript:AT1G01080.2
1  gene:AT1G01090  47234  49304  -  9  transcript:AT1G01090.1
1  gene:AT1G01100  49909  51210  -  10  transcript:AT1G01100.2

于是我用 sed 刪除了這些信息

sed -i -e 's/gene://' -e 's/transcript://' ath.gff

另外候学，pep.fa, cds.fa里面包含所有的基因藕筋，且命名特別的長， 同時我們也不需要基因中的 ".1", ".2" 這部分信息

$ head -n 5 Arabidopsis_thaliana.TAIR10.cds.all.fa 
>AT3G05780.1 cds chromosome:TAIR10:3:1714941:1719608:-1 gene:AT3G05780 gene_biotype:protein_coding transcript_biotype:protein_coding gene_symbol:LON3 description:Lon protease homolog 3, mitochondrial [Source:UniProtKB/Swiss-Prot;Acc:Q9M9L8]
ATGATGCCTAAACGGTTTAACACGAGTGGCTTTGACACGACTCTTCGTCTACCTTCGTAC
TACGGTTTCTTGCATCTTACACAAAGCTTAACCCTAAATTCCCGCGTTTTCTACGGTGCC
CGCCATGTGACTCCTCCGGCTATTCGGATCGGATCCAATCCGGTTCAGAGTCTACTACTC
TTCAGGGCACCGACTCAGCTTACCGGATGGAACCGGAGTTCTCGCGATTTATTGGGTCGTb

借助于seqkit, 我對原來的數(shù)據(jù)進行了篩選和重命名

seqkit grep -f <(cut -f 7 ath.gff ) Arabidopsis_thaliana.TAIR10.cds.all.fa | seqkit seq --id-regexp "^(.*?)\\.\\d" -i > ath.cds.fa
seqkit grep -f <(cut -f 7 ath.gff ) Arabidopsis_thaliana.TAIR10.pep.all.fa | seqkit seq --id-regexp "^(.*?)\\.\\d" -i > ath.pep.fa  

通過NCBI的BLASTP 或者DIAMOND 進行蛋白之間相互比對梳码，輸出格式為 -outfmt 6

makeblastdb -in ath.pep.fa -dbtype prot
blastp -num_threads 50 -db ath.pep.fa -query ath.pep.fa -outfmt 6 -evalue 1e-5 -num_alignments 20  -out ath.blastp.txt & 

共線性分析

繪制點陣圖

在上面步驟結(jié)束后隐圾，其實繪制點陣圖就非常簡單了，只需要創(chuàng)建配置文件掰茶，然后修改配置文件暇藏，最后運行wgdi。

基本wgdi的其他分析也都是三部曲濒蒋，創(chuàng)建配置盐碱，修改配置，運行程序沪伙。

第一步瓮顽，創(chuàng)建配置文件

wgdi -d \? > ath.conf

配置文件的信息如下(下面信息只是為了講解參數(shù)，不需要復制）

[dotplot]
blast = blast file
gff1 =  gff1 file
gff2 =  gff2 file
lens1 = lens1 file
lens2 = lens2 file
genome1_name =  Genome1 name
genome2_name =  Genome2 name
multiple  = 1   # 最好的同源基因數(shù), 用輸出結(jié)果中會用紅點表示
score = 100     # blast輸出的score 過濾 
evalue = 1e-5   # blast輸出的evalue 過濾 
repeat_number = 20  # genome2相對于genome1的最多同源基因數(shù)
position = order
blast_reverse = false
ancestor_left = none
ancestor_top = none
markersize = 0.5  # 點的大小
figsize = 10,10   # 圖片大小
savefig = savefile(.png,.pdf)

第二步围橡，修改配置文件暖混。因為是自我比對，所以 gff1和gff2內(nèi)容一樣某饰， lens1和lens2內(nèi)容一樣

[dotplot]
blast = ath.blastp.txt
gff1 =  ath.gff
gff2 =  ath.gff
lens1 = ath.len
lens2 = ath.len
genome1_name =  A. thaliana
genome2_name =  A. thaliana
multiple  = 1
score = 100
evalue = 1e-5
repeat_number = 5 
position = order
blast_reverse = false
ancestor_left = none
ancestor_top = none
markersize = 0.5
figsize = 10,10
savefig = ath.dot.png

最后運行如下命令

wgdi -d ath.conf

輸出的png如下所示

點陣圖

從圖中儒恋，我們很容易地就能觀察到三種顏色，分別是紅色黔漂，藍色诫尽，和灰色。在WGDI中炬守，紅色表示genome2的基因在genome1中的最優(yōu)同源（相似度最高）的匹配牧嫉，次好的四個基因是藍色，而剩余部分是為灰色。圖中對角線出現(xiàn)的片段并非是自身比對酣藻，因為WGDI已經(jīng)過濾掉自身比自身的結(jié)果曹洽。那么問題來了，這些基因是什么呢辽剧？這個問題會在后續(xù)的分析中進行解答送淆。

如果基因組存在加倍事件，那么對于一個基因組的某個區(qū)域可能會存在另一個區(qū)域與其有著相似的基因（同源基因）怕轿，并且這些基因的排布順序較為一致偷崩。反應到點陣圖上，就是能夠觀察到有多個點所組成的"線"撞羽。這里的線需要有一個引號阐斜，因為當你放大之后，你會發(fā)現(xiàn)其實還是點而已诀紊。

部分區(qū)域放大

由于進一步鑒定共線性區(qū)域就建立在這些"點"的基礎(chǔ)上谒出，那么影響這些點是否為同源基因的參數(shù)就非常重要了，即配置文件中的score,evalue,repeat_number邻奠。

共線性分析和Ks計算

對于同線性(synteny)和共線性(collinearity)笤喳，一般認為同線性指的是兩個區(qū)域有一定數(shù)量的同源基因，對基因順序排布無要求惕澎。共線性是同線性的特殊形式莉测，要求同源基因的排布順序也相似。

wgdi開發(fā)了-icl模塊（Improved version of ColinearScan）用于進行共線性分析唧喉，使用起來非常簡單，也是先建立配置文件忍抽。

wgdi -icl \? >> ath.conf

然后對配置文件進行修改.

[collinearity]
gff1 = ath.gff
gff2 = ath.gff
lens1 = ath.len
lens2 = ath.len
blast = ath.blastp.txt
blast_reverse = false
multiple  = 1
process = 8
evalue = 1e-5
score = 100
grading = 50,40,25
mg = 40,40
pvalue = 0.2
repeat_number = 10
positon = order
savefile = ath.collinearity.txt

其中的evalue, score和BLAST文件的過濾有關(guān)八孝，用于篩選同源基因?qū)Αultiple確定最佳比對基因鸠项，repeat_number表示最多允許多少基因是潛在的同源基因, grading根據(jù)同源基因的匹配程度進行打分, mg 指的是共線性區(qū)域中所允許最大空缺基因數(shù)干跛。

運行結(jié)果后，會得到ath.collinearity.txt祟绊，記錄共線性區(qū)域楼入。以 # 起始的行記錄共線性區(qū)域的元信息，例如得分(score)牧抽，顯著性(pvalue)嘉熊，基因?qū)?shù)(N)等。

grep '^#' ath.collinearity.txt | tail -n 1
# Alignment 959: score=778.0 pvalue=0.0003 N=16 Pt&Pt minus 

緊接著扬舒，我們就可以根據(jù)共線性結(jié)果來計算Ks阐肤。

同樣的也是先生成配置文件

wgdi -ks \? >> ath.conf

然后修改配置文件。我們需要提供cds和pep文件，共線性分析輸出文件(支持MCScanX的共線性分析結(jié)果）孕惜。WGDI會用muscle根據(jù)蛋白序列進行聯(lián)配愧薛，然后使用pal2pal.pl 基于cds序列將蛋白聯(lián)配轉(zhuǎn)成密碼子聯(lián)配，最后用paml中的yn00計算ka和ks衫画。

[ks]
cds_file =   ath.cds.fa
pep_file =   ath.pep.fa
align_software = muscle
pairs_file = ath.collinearity.txt
ks_file = ath.ks

ks計算需要一段時間毫炉，那么不妨在計算時了解下Ka和Ks。Ka和Ks分別指的是非同義替換位點數(shù)削罩，和同義替換位點數(shù)瞄勾。根據(jù)中性演化假說，基因的變化大多是中性突變鲸郊，不會影響生物的生存丰榴，那么當一對同源基因分開的時間越早，不影響生存的堿基替換數(shù)就越多(Ks)秆撮，反之越多四濒。

運行結(jié)果后，輸出的Ks文件共有6列职辨，對應的是每個基因?qū)Φ腒a和Ks值盗蟆。

$ head ath.ks
id1  id2  ka_NG86  ks_NG86  ka_YN00  ks_YN00
AT1G72300  AT1G17240  0.1404  0.629  0.1467  0.5718
AT1G72330  AT1G17290  0.0803  0.5442  0.079  0.6386
AT1G72350  AT1G17310  0.3709  0.9228  0.3745  1.071
AT1G72410  AT1G17360  0.2636  0.875  0.2634  1.1732
AT1G72450  AT1G17380  0.2828  1.1068  0.2857  1.5231
AT1G72490  AT1G17400  0.255  1.2113  0.2597  2.0862
AT1G72520  AT1G17420  0.1006  0.9734  0.1025  1.0599
AT1G72620  AT1G17430  0.1284  0.7328  0.1375  0.603
AT1G72630  AT1G17455  0.0643  0.7608  0.0613  1.2295

鑒于共線性和Ks值輸出結(jié)果信息量太大，不方便使用舒裤，更好的方法是將兩者進行聚合喳资，得到關(guān)于各個共線性區(qū)的匯總信息。

WGDI提供 -bi 參數(shù)幫助我們進行數(shù)據(jù)整合腾供。同樣也是生成配置文件

wgdi -bi ? >> ath.conf

修改配置文件. 其中collinearity 和ks是我們前兩步的輸出文件仆邓，而 ks_col 則是聲明使用ks文件里的哪一列。

[blockinfo]
blast = ath.blastp.txt
gff1 =  ath.gff
gff2 =  ath.gff
lens1 = ath.len
lens2 = ath.len
collinearity = ath.collinearity.txt
score = 100
evalue = 1e-5
repeat_number = 20
position = order
ks = ath.ks
ks_col = ks_NG86
savefile = ath_block_information.csv

運行即可

wgdi -bi ath.conf

輸出文件以csv格式進行存放伴鳖，因此可以用EXCEL直接打開节值，共11列。

1. id 即共線性的結(jié)果的唯一標識

2. chr1,start1,end1 即參考基因組（點圖的左邊）的共線性范圍

3. chr2,start2,end2 即參考基因組（點圖的上邊）的共線性范圍

4. pvalue 即共線性結(jié)果評估榜聂，常常認為小于0.01的更合理些

5. length 即共線性片段長度

6. ks_median 即共線性片段上所有基因?qū)?code>ks的中位數(shù)（主要用來評判ks分布的）

7. ks_average 即共線性片段上所有基因?qū)?code>ks的平均值

8. homo1,homo2,homo3,homo4,homo5與multiple參數(shù)有關(guān)搞疗，即一共有homo+multiple列

主要規(guī)則是：基因?qū)θ绻邳c圖中為紅色，賦值為1须肆，藍色賦值為0匿乃，灰色賦值為-1（顏色參考前述 wgdi 點圖 相關(guān)推文）。共線性片段上所有基因?qū)x值后求平均值豌汇，這樣就賦予該共線性一個-1,1的值幢炸。如果共線性的點大部分為紅色，那么該值接近于1瘤礁⊙舳可以作為共線性片段的篩選。

1. block1,block2分別為共線性片段上基因order的位置。

2. ks共線性片段的上基因?qū)Φ?code>ks值

3. density1,density2 共線性片段的基因分布密集程度岩调。值越小表示稀疏

這個表格是后續(xù)探索分析的一個重點巷燥，比如說我們可以用它來分析我們點圖中出現(xiàn)在對角線上的基因。下面的代碼就是提取了第一個共線性區(qū)域的所有基因号枕，然后進行富集分析缰揪，繪制了一個氣泡圖。

df <- read.csv("ath_block_information.csv")

tandem_length <- 200

df$start <- df$start2 - df$start1
df$end <- df$end2 - df$end1

df_tandem <- df[abs(df$start) <= tandem_length |
     abs(df$end) <= tandem_length,]

gff <- read.table("ath.gff")

syn_block1 <-  df_tandem[1,c("block1", "block2")]

gene_order <- unlist(
  lapply(syn_block1, strsplit, split="_", fixed=TRUE, useBytes=TRUE)
)

gene_order <- unique(as.numeric(gene_order))


syn_block1_gene <- gff[ gff$V6 %in% gene_order, "V2"]

#BiocManager::install("org.At.tair.db")
library(org.At.tair.db)
library(clusterProfiler)

ego <- enrichGO(syn_block1_gene, 
                OrgDb = org.At.tair.db, 
                keyType = "TAIR",ont = "BP"
                )
dotplot(ego)

氣泡圖

進一步觀察這些基因葱淳，可以發(fā)現(xiàn)這些基因的編號都是前后相連钝腺，說明這個基因簇和花粉識別有一定關(guān)系。

> as.data.frame(ego)[1,]
                   ID           Description GeneRatio  BgRatio       pvalue
GO:0048544 GO:0048544 recognition of pollen    11/560 43/21845 7.794138e-09
               p.adjust      qvalue
GO:0048544 8.121073e-06 7.95418e-06
                                                                                                                  geneID
GO:0048544 AT1G61360/AT1G61370/AT1G61380/AT1G61390/AT1G61400/AT1G61420/AT1G61440/AT1G61480/AT1G61490/AT1G61500/AT1G61550
           Count
GO:0048544    11

A gene cluster is a group of two or more genes found within an organism's DNA that encode similar polypeptides, or proteins, which collectively share a generalized function and are **often located within a few thousand base pairs **of each other. --維基百科

接下來赞厕，我們還可以根據(jù)這個表格中的Ks對點陣圖進行上色艳狐， 以及正確的計算Ks峰值。

如果對WGDI的使用有疑問皿桑，可以加入qq群

參考資料：
https://wgdi.readthedocs.io/en/latest/index.html