從零開始學(xué)CIRCOS繪制圈圖(四)

通常circos的中間部分不是空白區(qū)域亲善，會用一條條線進(jìn)行連接，表示兩個染色體部分區(qū)域有關(guān)系。

數(shù)據(jù)格式

對于link犯建，circos要求輸入數(shù)據(jù)至少有6列，分別是chr1 start1 end1 chr2 start2 end2 [options]

舉個例子

chr1    1000000 2000000 chr5    3000000 4000000

構(gòu)建輸入

這次會以A. lyrata 和 A.thalina的基因組為例瓜客，利用JCVI來構(gòu)建Circos的輸入.

新建一個項(xiàng)目文件夾

mkdir -p ath_aly && cd ath_aly

如果沒有安裝jcvi适瓦，用conda進(jìn)行安裝

conda create -n jcvi jcvi
conda activate jcvi
conda install last scipy
wget https://raw.githubusercontent.com/tanghaibao/jcvi-bin/master/bin/scip
chmod 755 scip
mv scip ~/miniconda3/bin

數(shù)據(jù)下載

從http://plants.ensembl.org/index.html分別下載這兩個物種的GFF文件和CDS序列，

# Athaliana
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.44.gff3.gz
# Alyrata
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_lyrata/cds/Arabidopsis_lyrata.v.1.0.cds.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_lyrata/Arabidopsis_lyrata.v.1.0.44.gff3.gz

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

將GFF3轉(zhuǎn)成BED格式

python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Arabidopsis_thaliana.TAIR10.44.gff3.gz > ath.bed
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Arabidopsis_lyrata.v.1.0.44.gff3.gz  > aly.bed

將BED去重復(fù)

python -m jcvi.formats.bed uniq ath.bed
python -m jcvi.formats.bed uniq aly.bed

提取CDS序列

seqkit grep -f <(cut -f 4 ath.bed ) Athaliana_167_TAIR10.cds.fa.gz | seqkit seq -i > ath.cds
seqkit grep -f <(cut -f 4 aly.bed )  Arabidopsis_lyrata.v.1.0.cds.all.fa.gz | seqkit seq -i > aly.cds

karyotype

從ENSEMBLE下載的GFF文件中谱仪，已經(jīng)包含了每個基因組的大小, 當(dāng)然也可以各種工具從基因組序列序列中獲取大小玻熙。

karyotype.aly.txt

chr -   aly1    aly1    0   33132539    rdylbu-11-div-1
chr -   aly2    aly2    0   19320864    rdylbu-11-div-2
chr -   aly3    aly3    0   24464547    rdylbu-11-div-3
chr -   aly4    aly4    0   23328337    rdylbu-11-div-4
chr -   aly5    aly5    0   21221946    rdylbu-11-div-5
chr -   aly6    aly6    0   25113588    rdylbu-11-div-6
chr -   aly7    aly7    0   24649197    rdylbu-11-div-7
chr -   aly8    aly8    0   22951293    rdylbu-11-div-8

karyotype.tair10.txt

chr -   ath1    ath1    0   30427671    brbg-10-div-1
chr -   ath2    ath2    0   19698289    brbg-10-div-3
chr -   ath3    ath3    0   23459830    brbg-10-div-5
chr -   ath4    ath4    0   18585056    brbg-10-div-7
chr -   ath5    ath5    0   26975502    brbg-10-div-9

因?yàn)镋NSMEBLE上Athalina和Alyrata的染色體命名都是1,2,3…，就會導(dǎo)致CIRCOS無法正確的區(qū)分來源疯攒，因此在原本的命名前加上了物種名縮寫做為標(biāo)簽嗦随。

同時。我們要修改之前的bed文件

sed -e 's/^/ath/' ath.uniq.bed  > ath.bed
sed -e 's/^/aly/' aly.uniq.bed  > aly.bed

links

為了構(gòu)建links文件敬尺，需要利用JCVI進(jìn)行共線性分析.

確保有4個文件

$ ls ???.???
aly.bed aly.cds ath.bed ath.cds

用jcvi進(jìn)行分析

python -m jcvi.compara.catalog ortholog --no_strip_names ath aly
python -m jcvi.compara.synteny screen --minspan=30 --simple ath.aly.anchors ath.aly.anchors.new

其中ath.aly.anchors.simple是我們后續(xù)要用到的文件枚尼。

$ head -n 1 ath.aly.anchors.simple
AT1G24260.1 AT1G27280.1 fgenesh1_pm.C_scaffold_1002045  fgenesh2_kg.1__2877__AT1G24260.1    225 -

我們需要將ath.aly.anchors.simple里的基因名替換成實(shí)際的位置信息，將其變成符合Circos的輸入信息砂吞。

我寫了一個simple2links.py腳本署恍，代碼在我的GitHub上，https://github.com/xuzhougeng/myscripts蜻直。

python ~/myscripts/simple2links.py ath.aly.anchors.simple

最終會輸出ath.aly.anchors.simple_link.txt

配置circos

接下來做如下配置盯质。新建一個etc文件夾，在里面建立一個links.conf用來配置links, 建立一個ticks.conf配置ticks

mkdir etc
# vim etc/links.conf
<links>

<link>
file          = ath.aly.anchors.simple_link.txt
radius        = 0.61r
color         = blue_a4
ribbon = yes
</link>

</links>

# vim etc/ticks.conf

show_ticks          = yes
show_tick_labels    = yes

<ticks>

radius           = 1r
color            = black
thickness        = 2p
multiplier       = 1e-6
format           = %d

<tick>
spacing        = 1u
size           = 5p
</tick>

<tick>
thickness      = 4p
spacing        = 5u
size           = 10p
show_label     = yes
label_size     = 10p
label_offset   = 10p
format         = %d
</tick>

</ticks>

編輯circos.conf

karyotype = karyotype.tair10.txt,karyotype.aly.txt

chromosomes_color = chr1=rdylbu-11-div-1,chr2=rdylbu-11-div-3,chr3=rdylbu-11-div-5,chr4=rdylbu-11-div-7,chr5=rdylbu-11-div-9

chromosomes_units = 1000000
<<include ./etc/ticks.conf>>

<ideogram>
<spacing>
default = 0.005r
</spacing>
radius           = 0.90r
thickness        = 20p
fill             = yes
stroke_color     = dgrey
stroke_thickness = 2p

show_label     = yes #展示label
label_font     = default # 字體
label_radius   = dims(ideogram,radius) + 0.08r #位置
label_size     = 16 # 字體大小
label_parallel = yes # 是否平行

label_format   = eval(sprintf("%s",var(chr))) # 格式
</ideogram>

<<include ./etc/links.conf>>

<image>
dir*    = .    # 輸出文件夾
radius* = 500p # 圖片半徑
svg*    = no   # 是否輸出svg
<<include etc/image.conf>>
</image>

<<include etc/colors_fonts_patterns.conf>>
<<include etc/housekeeping.conf>>

運(yùn)行circos -conf circos.conf的效果如下

circos-links-1

不難發(fā)現(xiàn)一個問題概而，里面線的顏色一模一樣呼巷，不容易進(jìn)行區(qū)分。雖然可以在輸入ath.aly.anchors.simple_link.txt里增加顏色赎瑰，但是circos提供了一個動態(tài)規(guī)則王悍，可以更加方便的直接在配置文件里修改。

建立規(guī)則(rules)

circos的配置格式為

<rules>

<rule>
...
</rule>

<rule>
...
</rule>
...

</rules>

circos的規(guī)則可以很復(fù)雜乡范，但是最簡單的情況就是下面這種

<rule>
condition = var(chr1) eq "ath1"
color=rdylgn-5-div-1
</rule>

condition = var(chr1) eq "ath1"表示配名，判斷l(xiāng)ink文件中左側(cè)染色體的名字(var(chr1))是不是(eq)"ath1"啤咽，如果是的話，那么顏色就是rdylgn-5-div-1

我們可以在etc/links.conf中增加五個條件渠脉，修改后的links.conf如下

<link>
file          = ath.aly.anchors.simple_link.txt
radius        = 0.61r
color         = blue_a4
ribbon = yes

<rules>
<rule>
condition = var(chr1) eq "ath1"
color=rdylgn-5-div-1
</rule>
<rule>
condition = var(chr1) eq "ath2"
color=rdylgn-5-div-2
</rule>
<rule>
condition = var(chr1) eq "ath3"
color=rdylgn-5-div-3
</rule>
<rule>
condition = var(chr1) eq "ath4"
color=rdylgn-5-div-4
</rule>
<rule>
condition = var(chr1) eq "ath5"
color=rdylgn-5-div-5
</rule>
</rules>

</link>
</links>