从零开始学CIRCOS绘制圈图(四)

通常circos的中间部分不是空白区域,会用一条条线进行连接,表示两个染色体部分区域有关系。

数据格式

对于link,circos要求输入数据至少有6列,分别是chr1 start1 end1 chr2 start2 end2 [options]

举个例子

chr1	1000000	2000000	chr5	3000000	4000000

构建输入

这次会以A. lyrata 和 A.thalina的基因组为例,利用JCVI来构建Circos的输入.

新建一个项目文件夹

mkdir -p ath_aly && cd ath_aly

如果没有安装jcvi,用conda进行安装

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

数据下载

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

数据预处理

将GFF3转成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去重复

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文件中,已经包含了每个基因组的大小, 当然也可以各种工具从基因组序列序列中获取大小。

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

因为ENSMEBLE上Athalina和Alyrata的染色体命名都是1,2,3…,就会导致CIRCOS无法正确的区分来源,因此在原本的命名前加上了物种名缩写做为标签。

同时。我们要修改之前的bed文件

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

为了构建links文件,需要利用JCVI进行共线性分析.

确保有4个文件

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

用jcvi进行分析

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是我们后续要用到的文件。

$ 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里的基因名替换成实际的位置信息,将其变成符合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>>

运行circos -conf circos.conf的效果如下

circos-links-1

不难发现一个问题,里面线的颜色一模一样,不容易进行区分。虽然可以在输入ath.aly.anchors.simple_link.txt里增加颜色,但是circos提供了一个动态规则,可以更加方便的直接在配置文件里修改。

建立规则(rules)

circos的配置格式为

<rules>

<rule>
...
</rule>

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

</rules>

circos的规则可以很复杂,但是最简单的情况就是下面这种

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

condition = var(chr1) eq "ath1"表示,判断link文件中左侧染色体的名字(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>

运行circos -conf circos.conf的效果如下

circos-link-2

这张图还有一个问题是,通常别人的Circos染色体都是对称排列,右边第一个是ath1,那么左边第一个也最好是aly1,那么应该如何配置呢?以及两套基因组会分别占据两侧,中间有一个明显的空隙,如何做到的呢?

评论

Your browser is out-of-date!

Update your browser to view this website correctly. Update my browser now

×