在linux系统上使用circos绘制圈图。
- circos的安装
安装是需要管理员权限的,需要添加到sudo
组中,我这里我的权限就是sudo
:
进入软件路径cd 01.software
创建circos的存放目录mkdir circos_soft && cd circos_soft
wget http://circos.ca/distribution/circos-0.69-6.tgz
tar xzf circos-0.69-6.tgz
创建一个软连接:
ln -s circos-0.69-3 current
检查缺少的安装模块:
cd current/bin
test.modules |grep missing> ../../missing_modules
cd ../../
less missing_modules
可以看到缺失的模块名字:
对缺失模块进行安装:
for mod in $(cat missing_modules);
do
cpan install $mod;
done
当提示需要高权限的时候,切换为sudo账号即可。
模块安装完就可以使用了,将软件所在路径添加到.bashrc,以便在哪个路径下都能直接使用,要不然必须使用绝对路径:
绝对路径就是这样:
而在~/.bashrc中加入:
PATH="/mnt/home/yanyt/01.software/circos_soft/current/bin:$PATH"
系统在哪里都会认识circos
,原理是将circos
与/mnt/home/yanyt/01.software/circos_soft/current/bin/
拼接在了一起,于是:
- circos的使用
这里我是基于宏基因组学的go功能注释文件的,先看看我的文件结构less go.tsv
:
其中列名是通路和样本名字,第一列是go通路,第2列到4列是样本的名字。
进行筛选,将未注释的通路所在行以及子通路所在的行删除:
less -S go.tsv|grep -v UNMAPPED|grep -v UNGROUPED|grep -v "|" >go_1.tsv
sed -i 's/# Gene Family/pathway/g' go_1.tsv
查看go_1.tsv:
而circos是有特定的文件输入格式,我的目的需要根据go_1.tsv的信息,构建circos的输入文件
nodes.tsv
和lines.tsv
。先看看构建好的文件结构:
less nodes.tsv
其中第一列和第二列是固定的,万年不用更改,第三列是每个节点的特有ID,我是直接字符串拼接的,第四列是节点的名字,样本名和丰度前十的通路名字,第五列和第六列是起始位置和终止位置,这里第五列都是0,而第六列取决于这个节点与其他多少个节点之间有联系,最后一列是颜色的信息。
less lines.tsv
lines.tsv是存储节点之间连线的信息,第一列是节点的ID名字,第2列和3列是起始节点的起始位置和终止位置,对应于nodes.tsv中的第五列和第六列,第4列是终止节点的ID,第5和6列则是终止节点的起始和终止位置。
那么如何通过go_1.tsv那样的表格,得到nodes.tsv和lines.tsv呢?在这里,我写了一个r脚本:
# 默认的quote会跳过2/3的数据,导致行减少产生NA,改默认值为空
taxonomy = read.table("go_1.tsv", header=T, sep="\t", quote = "", row.names=NULL, comment.char="")
print(paste0("All taxonomy annotations are ", dim(taxonomy)[1], " lines!"))
# 去除NA,否则无法计算
taxonomy = na.omit(taxonomy)
# 显示样本总数据,有冗余
# colSums(taxonomy)
grp = taxonomy[, "pathway", drop=F]
abu = taxonomy[,2:dim(taxonomy)[2]]
merge = cbind(abu, grp)
# group_by传变量,前面加".dots="
mergeTax = merge %>% group_by(.dots="pathway") %>% summarise_all(sum)
# 合并后表格转换为数据框
mergeTax = as.data.frame(mergeTax)
# 按丰度排序
idx = order(rowMeans(mergeTax[,2:dim(mergeTax)[2]]), decreasing = T)
mergeTax = mergeTax[idx,]
# 添加行名
rownames(mergeTax)=mergeTax[,1]
# remove rownames line
Top = mergeTax[,-1]
# normalization to percentage 100
Top = as.data.frame(t(t(Top)/colSums(Top,na=T) * 100))
# Select TopN line for plotting
Top = head(Top, n=10)
my_source<-c()
my_target<-c()
for(i in colnames(Top)){
sub<-Top[,i,drop=F]
sub[sub==0]<-NA
sub=na.omit(sub)
my_target<-c(my_target,rownames(sub))
my_source<-c(my_source,rep(colnames(sub),length(rownames(sub))))
}
my_table<-data.frame(source=my_source,target=my_target)
aa<-c()
bb<-c()
cc<-c()
dd<-c()
ee<-c()
ff<-c()
gg<-c()
j=14
for(i in unique(my_table$source)){
sub<-my_table[my_table$source==i,]
ff<-c(ff,length(rownames(sub)))
ee<-c(ee,0)
dd<-c(dd,i)
cc<-c(cc,paste0("ID",j))
bb<-c(bb,"-")
aa<-c(aa,"chr")
gg<-c(gg,paste0("chr",j))
j=j+1
}
j=4
for(i in unique(my_table$target)){
sub<-my_table[my_table$target==i,]
ff<-c(ff,length(rownames(sub)))
ee<-c(ee,0)
dd<-c(dd,i)
cc<-c(cc,paste0("ID",j))
bb<-c(bb,"-")
aa<-c(aa,"chr")
gg<-c(gg,paste0("chr",j))
j=j+1
}
colors<-data.frame(cc=paste0("ID",4:20),gg=c("255,41,41","255,63,63","0,255,0","0,153,237","151,143,255","255,182,106","255,244,80","0,204,0","237,0,0","0,136,220","255,41,41","255,63,63","0,255,0","0,153,237","151,143,255","255,182,106","0,204,0"))
table_1<-data.frame(aa=aa,bb=bb,cc=cc,dd=dd,ee=ee,ff=ff,gg=gg)
table_1$gg<-colors[match(table_1$cc,colors$cc),]$gg
write.table(table_1,file=paste( "nodes.tsv", sep = ""), append = FALSE, sep="\t", quote=F, row.names=F, col.names=F)
my_table$source_ID<-table_1[match(my_table$source,table_1$dd),"cc"]
my_table$target_ID<-table_1[match(my_table$target,table_1$dd),"cc"]
my_table$source_aa<-table_1[match(my_table$source,table_1$dd),"ee"]
my_table$source_bb<-table_1[match(my_table$source,table_1$dd),"ff"]
my_table$target_aa<-table_1[match(my_table$target,table_1$dd),"ee"]
my_table$target_bb<-table_1[match(my_table$target,table_1$dd),"ff"]
my_table<-my_table[,c("source_ID","source_aa","source_bb","target_ID","target_aa","target_bb")]
write.table(my_table,file=paste("lines.tsv", sep = ""), append = FALSE, sep="\t", quote=F, row.names=F, col.names=F)
其中有一些小细节,是基于我自己的需求,比如通路的ID是从ID3到ID12,这是为了便于我在后续控制连向通路的线条的颜色与通路的块的颜色始终保持一致。
接下来就是配置circos的文件,在lines.tsv与nodes.tsv所在路径下,vim circos.conf
输入以下配置信息:
#指定染色体文件(绝对/相对路径+文件名)
karyotype = nodes.tsv
chromosomes_units = 1 #定义最小单位,即100万bp为一个units,即1u = 100w,后面刻度线都是基于此的操作,如果染色体长度都在5kw以上,推荐用100w,>否则推荐用10w
#-----------------------------------------------------------------------------------
<ideogram> #这是定义染色体相关参数的标签,以</ideogram>结尾,其中包括要设置的参数
<spacing> #定义染色体间隙宽度的标签,以</spacing>,其中包括要设置的参数
default = 0.005r #r指的是圆的周长,设置0.5%圆的周长为间隙
#<pairwise hsY;hs1> #可以用<pairwise>标签特别指定某些染色体的间隙(用的是ID),因为在大多数文章中,都会留一个大间隙,来放label
#spacing = 20r #这里20r表示是相对default = 0.005r的20倍,也就是10%的圆的周长
#</pairwise> #标签都要以</>结尾,
</spacing> #间隙定义结束,下面是对染色体样式的调整
radius = 0.80r #轮廓的位置,这里的r指的是半径,由圆心到圆周上范围依次是0-1r,,超出部分将不再显示。
thickness = 40p #染色体整体的宽度,这里p指的是像素大小,也可以用r表示,1r=1500p
fill = yes #是否为染色体填充颜色,如果为yes,自动用第七列定义的颜色着色
stroke_color = dgrey #染色体边框的颜色,支持多种格式的输入,如:red或255,182,106
stroke_thickness = 2p #染色体边框的粗细
show_label = yes #选择yes表示要显示label
label_font = default # 字体可以再 etc/fonts.conf 查看所有,默认为CMUBright-Roman
label_radius = dims(image,radius)-60p #使用dims()获取图像半径的大小,从而定位染色体标签的位置,也可以直接定义1.075r
label_size = 30 #字体的大小
label_parallel = yes #将Label的方向设置为与染色体平行
</ideogram> #定义染色体属性的标签结束
show_ticks = yes #选择yes表示要显示刻度线
show_tick_labels = yes #选择yes表示要显示刻度线的数值
#定义刻度线的整体位置与形状
<ticks> #刻度线的转用标签,但凡是复数出现的,其下面的参数都表示全局参数,像下面的<tick>单数形式,都表示局部参数
radius = 1r #刻度线的位置,1r为最远距离,超过1r不再显示
color = black
thickness = 2p
multiplier = 1 #把刻度线标签(bp)缩小10万倍显示
format = %d #然后以整数的形式标记在刻度线上
#定义小的刻度线,且不显示数值
<tick>
spacing = 0.5u #最开始我们定义1u = 1000000,表示每500w bp显示一个小刻度线
size = 10p
show_label = no #由于小的刻度线展示出来太密集,因此我们no不展示,默认不展示
</tick>
#定义大的刻度线,显示数值
<tick>
spacing = 1u
size = 15p
show_label = no
#label_size = 20p
#label_offset = 10p #设置数值和刻度线之间的间隔
#format = %d
</tick>
</ticks>
<links>
<link>
file=lines.tsv
radius = 0.7r
bezier_radius = 0r
crest = 1
thickness = 5p
#color = eval(var(chr2))
#color = eval(sprintf("chr%s", var(chr)))
#color=blue
<rules>
<rule>
condition = var(chr2) eq "ID4"
color=255,41,41
</rule>
<rule>
condition = var(chr2) eq "ID5"
color=255,63,63
</rule>
<rule>
condition = var(chr2) eq "ID6"
color=0,255,0
</rule>
<rule>
condition = var(chr2) eq "ID7"
color=0,153,237
</rule>
<rule>
condition = var(chr2) eq "ID8"
color=151,143,255
</rule>
<rule>
condition = var(chr2) eq "ID9"
color=255,182,106
</rule>
<rule>
condition = var(chr2) eq "ID10"
color=255,244,80
</rule>
<rule>
condition = var(chr2) eq "ID11"
color=0,204,0
</rule>
<rule>
condition = var(chr2) eq "ID12"
color=237,0,0
</rule>
<rule>
condition = var(chr2) eq "ID13"
color=0,136,220
</rule>
</rules>
</link>
</links>
#-----------------------------------------------------------------------------------
#下面是每次都要复制粘贴上去的,他们属于circos自带的配置文件,用于调用颜色,距离,报错等信息
<image> #注意路径
<<include /mnt/home/yanyt/01.software/circos_soft/current/etc/image.conf>> #注意引用外部配置文件需要使用<<#>>
</image>
<<include /mnt/home/yanyt/01.software/circos_soft/current/etc/colors_fonts_patterns.conf>>
#官方没有提到下面的文件,但是没有这个文件会报错,所以还是加上去
<<include /mnt/home/yanyt/01.software/circos_soft/current/etc/housekeeping.conf>>
最后运行命令:circos -conf circos.conf
即可。
总结:circos也有windows的图形化界面,到官网下载windows版本即可。