生信 | circos文件准备大全（以拟南芥为例）

一、写在前面

$\color{red}{纯手撸linux代码，首次公开，不依赖SPDE和TBtools或其他语言脚本}$
$\color{red}{学完本教程，你将了解所有能用于circos的数据是如何产生的}$
$\color{red}{如果依然没有你想要的数据格式，请加群说出你的诉求，我将1天内更新完毕}$
$\color{red}{当然，需要你有一点点Linux知识，不过小白也没事，掉头发的事情我已经做完了}$
扣扣群：559758504（目前无需验证，一键即可加入）
如果你不想学，老板又催着要，你可以加群私聊我寻求帮助( $ _ $ )，但还是自己学，知识都是自己的。
实战代码，请看：【生信 | Circos实战】

目前已有格式总结（欢迎补充）：
1.染色体文件
2.统计分隔区间
3.统计基因密度(序列重复密度/TE转座子密度/RNA分布密度均同理）
4.统计GC含量
5.统计区间内基因表达/SNP/peak丰度等（只要能测序的都可以）
6.共线性文件（种内，种间都可）
7.文字文件（常见基因ID标注）

二、准备工作（可跳过）

1. 安装必须的软件（均在conda环境中安装）：

安装Circos、bedtools、BLAST、pyfaidx、samtools（装过的conda会自动跳过）

conda install -c bioconda -y circos 
conda install -c bioconda -y bedtools 
conda install -c bioconda -y blast 
conda install -c bioconda -y pyfaidx
conda install -c bioconda -y samtools

安装MCScanX（要做共线性分析还需要安装下面的软件，不需要则略过）

wget https://archive.fastgit.org/wyp1125/MCScanX/archive/refs/heads/master.zip
unzip master.zip
cd MCScanX-master
make
# 测试安装是否成功（能打印出来help信息则表示成功）
./MCScanX -h
# 成功则添加到环境变量
echo 'export PATH=$PATH:'"$(pwd)" >> ~/.bashrc
source ~/.bashrc

2. 下载拟南芥genome，gff3，pep文件

wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_chromosome_files/TAIR10_chr_all.fas
wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_gff3/TAIR10_GFF3_genes.gff
wget https://www.arabidopsis.org/download_files/Genes/TAIR10_genome_release/TAIR10_blastsets/TAIR10_pep_20101214_updated
# 微调一些官网上的文件格式
cut -d " " -f1 TAIR10_chr_all.fas | \
awk '{if($0~">"){print $0}else{printf "%s",$1}}' | \
sed 's/>/\n>/g' | grep -A1 -E ">[0-9]" | \
sed 's/>/>Chr/g' > tair10.genome.fa
cut -d " " -f1 TAIR10_pep_20101214_updated.txt > tair10.pep.fa
# 改名并删除没用的文件
mv TAIR10_GFF3_genes.gff tair10.gff3
rm TAIR10_chr_all.fas TAIR10_pep_20101214_updated.txt

三、文件准备

先说明一点：circos中的数据都以\t分割，数据表无表头

1. 染色体文件，格式如下：

前2列：固定为`chr -`，表明我们要画染色体轮廓
第3列：染色体ID，像身份证号，`注意：如果ID中除了“_”有其他符号一律会报错`
第4列：染色体Label，像人名，会显示在图上而ID不会
第5列：染色体起始位置
第6列：染色体终止位置
第7列：染色体的颜色，chr1-4,chry代表的都是circos内置颜色的名称
chr -   Chr1    C1  0   30427671    chr1
chr -   Chr2    C2  0   19698289    chr2
chr -   Chr3    C3  0   23459830    chr3
chr -   Chr4    C4  0   18585056    chr4
chr -   Chr5    C5  0   26975502    chr5

过程

genome=tair10.genome.fa
faidx ${genome} -i chromsizes > tair10.genome.size
awk -vFS="\t" -vOFS="\t" '{print "chr","-",$1,"C"NR,"0",$2,"chr"NR}' tair10.genome.size > circos.chromosome.txt
cat circos.chromosome.txt
# 结果如上表，由于chr1-5的颜色过于丑陋，因此我们'手动'修改为好看的颜色，
# 后面就可以根据染色体的颜色对其他物体着色啦，如下：
chr -   Chr1    C1  0   30427671    142,212,202
chr -   Chr2    C2  0   19698289    180,224,101
chr -   Chr3    C3  0   23459830    247,160,151
chr -   Chr4    C4  0   18585056    131,178,210
chr -   Chr5    C5  0   26975502    253,130,115

-vFS="\t"：输入文件的分隔符\t
-vOFS="\t"：输出文件的分隔符\t
NR：表示行号，正好可以和染色体的123等对应上，因此可以拿来用，但是有x,y染色体的时候需要自己再手动调整一下

2. 统计分隔区间，格式如下：

第1列：染色体ID
第2列：起始分隔点
第3列：终止分割点
Chr1    0   100000
Chr1    100000  200000
Chr1    200000  300000

过程

genomeSize=tair10.genome.size
windowSize=100000
bedtools makewindows -g ${genomeSize} -w ${windowSize} | \
awk -vFS="\t" -vOFS="\t" '{print $1,$2,$3}' | \
bedtools sort -i - > ${genomeSize}.100kb

windowSize：要根据你基因组的大小，动态调整，因为拟南芥的每条染色体的序列长度都小于5kw，因此我这里设置10万为一个小区间，后面所有的数据统计都是基于此区间统计的，因此合适的区间选择十分重要。如果你研究的基因组每个都大于5kw，则可以设置为100万为一个区间。
bedtools sort：由于后面要使用bedtools做一些统计，因此需要排序，不过默认也是排好的，以防万一。

3. 统计基因密度，格式如下：

第1列：染色体ID
第2列：在染色体上的起始位置
第3列：在染色体上的终止位置
第4列：区间内基因的数目
Chr1    0   100000  29
Chr1    100000  200000  36
Chr1    200000  300000  31

过程

# 定义一些变量，后面修改方便
gff3=tair10.gff3
genomeSize=tair10.genome.size.100kb
grep '[[:blank:]]gene[[:blank:]]' ${gff3} | \
awk '{print $1"\t"$4"\t"$5}'| \
bedtools coverage -a ${genomeSize} -b - | \
cut -f 1-4 > circos.gene.density.txt
# 按理来说，上面的circos.gene.density.txt已经足够了，不过我们可以修改文件批量美化
# 增加第5列，颜色信息（非必须，可选）
awk -vFS="\t" -vOFS="\t" 'NR==FNR{arr[$3]=$3;brr[$3]=$7}NR!=FNR{if(arr[$1]=$1){print $0,"fill_color="brr[$1]}}' circos.chromosome.txt circos.gene.density.txt > circos.gene.density.color.txt
# 结果文件形式如下
Chr1    0   100000  29  fill_color=142,212,202
Chr1    100000  200000  36  fill_color=142,212,202
Chr1    200000  300000  31  fill_color=142,212,202
Chr1    300000  400000  31  fill_color=142,212,202

FNR：每打开一个文件便重新计数，而NR会一直累计。推荐阅读：【awk之NR==FNR问题】

\color{red}{由于统一设置颜色很简单，而用perl控制颜色又复杂，所以后面就着重讲解如何在文件中配置颜色}

再啰嗦一句，使用这种方法，还可以刻意的突出你想展示的区域，对于不会perl的同学灵活度很大，而且还算比较方便

4. 统计GC含量，格式如下：

第1列：染色体ID
第2列：在染色体上的起始位置
第3列：在染色体上的终止位置
第4列：区间内GC碱基的数目
Chr1    0   100000  35627
Chr1    100000  200000  37655
Chr1    200000  300000  37903

代码

genome=tair10.genome.fa
genomeSize=tair10.genome.size.100kb
bedtools nuc -fi ${genome} -bed ${genomeSize} | \
cut -f 1-3,5 | grep -v "#" | \
awk -vFS="\t" -vOFS="\t" '{print $1,$2,$3,$4*($3-$2)}' > circos.gc.density.txt
# 同样在第5列添加和对应染色体相同的颜色（非必须，可选）
awk -vFS="\t" -vOFS="\t"  'NR==FNR{arr[$3]=$3;brr[$3]=$7}NR!=FNR{if(arr[$1]=$1){print $0,"fill_color="brr[$1]}}' circos.chromosome.txt circos.gc.density.txt > circos.gc.density.color.txt
# 结果文件形式如下
Chr1    0   100000  35627   fill_color=142,212,202
Chr1    100000  200000  37655   fill_color=142,212,202
Chr1    200000  300000  37903   fill_color=142,212,202

awk：计算GC含量的时候，使用$4*($3-$2)可以有效避免染色体末端含量异常，当然有时候末端比较短依然会有异常值，那仅保留$4再次尝试一下

$\color{red}{其实，circos的柱形图，热图，点图，折线图等所需的格式都是相同的，因此这里不再展开其他数据}$
注：其他格式如，区域内基因的表达丰度，SNP丰度，peak的丰度，等等等，你能转换为上面的格式就行，算了，还是再简单的举一个基因表达的例子吧

5. 统计区间内基因表达/SNP/peak丰度等，以基因表达为例，格式如下：

第1列：染色体ID
第2列：在染色体上的起始位置
第3列：在染色体上的终止位置
第4列：区间内基因表达丰度（也就是reads数目的多少）
Chr1    0   100000  3562
Chr1    100000  200000  3765
Chr1    200000  300000  3790

代码

# 主要就是对比对后的bam文件，进行的区间统计
genomeSize=tair10.genome.size.100kb
bam=SRR3081110.bam
samtools index ${bam}
bedtools multicov -bams ${bam} -bed ${genomeSize}  > circos.RNA.level.txt
# 结果文件格式和上面一样
# 如果还想添加第5列颜色信息，也可以加上（非必须）
awk -vFS="\t" -vOFS="\t"  'NR==FNR{arr[$3]=$3;brr[$3]=$7}NR!=FNR{if(arr[$1]=$1){print $0,"fill_color="brr[$1]}}' circos.chromosome.txt circos.RNA.level.txt > circos.RNA.level.color.txt
# 结果文件形式如下，一般不用 fill_color，因为丰度一般用热图表示最好，所以这一步不用做
Chr1    0   100000  3562   fill_color=142,212,202
Chr1    100000  200000  3765   fill_color=142,212,202
Chr1    200000  300000  3790   fill_color=142,212,202

$\color{red}{其实，circos的最后一列可以是可用的任意属性}$
$\color{red}{不仅可以调fill\ color还可以调color等等，凡是属性的值都能加到最后一列，多个属性要用【，】隔开}$

6. 共线性文件，格式如下：

第1列：染色体ID，指连线起始的染色体
第2列：连线起始的染色体的起始位置
第3列：连线起始的染色体的终止位置
第4列：染色体ID，指连线的终止染色体
第5列：连线终止的染色体的起始位置
第6列：连线终止的染色体的终止位置
Chr1    6265416 6266937 Chr1    27686992    27688127
Chr1    6331398 6333743 Chr1    27759973    27761588

物种内做共线性
这里不解释细节，看原理戳这里【物种内共线性分析】，或有疑问，加群提出
不过，有一点我再提一下，因为大多数教程都没有说，就是不要一味的使用MCScanX+前缀而不自己设定参数，-k,-g,-s,-m参数很重要，会极大的影响结果，它默认的并不一定是最好的
2021年12月2日更新：for循环可以并行了，速度提升约100倍，几乎秒出结果。

mkdir blastdb
mkdir blastresult
pep=tair10.pep.fa
gff3=tair10.gff3
makeblastdb -in ${pep} -dbtype prot -out ./blastdb/tair10
blastp -query ${pep} -db ./blastdb/tair10 -out ./blastresult/tair10.blast -evalue 1e-10 -num_threads 10 -outfmt 6 -num_alignments 5
# 格式化gff3文件，后面用的gff3都是格式化好的，千万注意
awk -vFS="\t" -vOFS="\t" '{if($3=="mRNA"){match($9,/ID=([^;]+)/,a);sub(/ID=/,"",a[0]);print $1,a[0],$4,$5}}' ${gff3} > ./blastresult/tair10.gff
cd blastresult
MCScanX ./blastresult/tair10
# 为了能从collinearity中提取用于circos的数据格式，这里讲解两种方法
# 第一种：以单个基因对单个基因来做（比较简单，但最后circos图很乱）
cd ./blastresult
grep -v "#" tair10.collinearity | cut -f2,3 | \
awk -vFS="\t" -vOFS="\t" 'NR==FNR{arr[$2]=$2;brr[$2]=$1"\t"$3"\t"$4}NR!=FNR{if(arr[$1]==$1)print brr[$1],$2}' tair10.gff - | \
awk -vFS="\t" -vOFS="\t" 'NR==FNR{arr[$2]=$2;brr[$2]=$1"\t"$3"\t"$4}NR!=FNR{if(arr[$4]==$4)print $1,$2,$3,brr[$4]}' tair10.gff - > ../circos.link.txt
# 增加颜色信息（非必须，但推荐做）
awk -vFS="\t" -vOFS="\t"  'NR==FNR{arr[$3]=$3;brr[$3]=$7}NR!=FNR{if(arr[$1]=$1){print $0,"color="brr[$1]}}' ../circos.chromosome.txt ../circos.link.txt > ../circos.link.color.txt
# 最终文件格式如下
Chr1    6265416 6266937 Chr1    27686992    27688127    color=142,212,202
Chr1    6331398 6333743 Chr1    27759973    27761588    color=142,212,202
Chr1    6336528 6342460 Chr1    27770986    27776857    color=142,212,202
# 第二种：以block来做（使用circos的丝带效果更好看，更推荐，但相对来说难理解一点，不过你能跑通就行，掉头发的事情我来做）
# 你只要是按照上面的流程跑下来的，定义好下面的这两个变量，直接执行for循环即可
collinearity=tair10.collinearity
gff=tair10.gff
# ------------------------------------------------------------
for i in $(tail -1 ${collinearity} | cut -d "-" -f1);do
for ((j=0;j<=$i;j++));do
{ grep -w "${j}-" ${collinearity} | cut -f2,3 | \
awk -vFS="\t" -vOFS="\t" 'NR==FNR{arr[$2]=$2;brr[$2]=$1"\t"$3"\t"$4}NR!=FNR{if(arr[$1]==$1)print brr[$1],$2}' ${gff} - | \
awk -vFS="\t" -vOFS="\t" 'NR==FNR{arr[$2]=$2;brr[$2]=$1"\t"$3"\t"$4}NR!=FNR{if(arr[$4]==$4)print $1,$2,$3,brr[$4]}' ${gff} - > block.${j}
chr_l=$(head -1 block.${j}|cut -f1)
min_l=$(sort -k2n block.${j}|head -1|cut -f2)
max_l=$(sort -k2nr block.${j}|head -1|cut -f3)
chr_r=$(head -1 block.${j}|cut -f4)
min_r=$(sort -k5n block.${j}|head -1|cut -f5)
max_r=$(sort -k6nr block.${j}|head -1|cut -f6)
echo -e ${chr_l}"\t"${min_l}"\t"${max_l}"\t"${chr_r}"\t"${min_r}"\t"${max_r} >> block.merge; }&
done; wait; clear
cp block.merge ../circos.link.block.txt
rm block.*
echo -e "\033[1;32mSuccess: \033[36m you can find it from [../circos.link.block.txt]\033[0m"
done
# ------------------------------------------------------------
# 最终文件(../circos.link.block.txt)格式如下（可以看出，区间更大了）
Chr1    6265416 7292507 Chr1    27686992    28717416
Chr1    8845231 9469634 Chr1    25524718    26389916
Chr1    5896416 6179289 Chr1    27217477    27603317
# 当然依旧可以增加颜色信息（非必须，但推荐做）
awk -vFS="\t" -vOFS="\t"  'NR==FNR{arr[$3]=$3;brr[$3]=$7}NR!=FNR{if(arr[$1]=$1){print $0,"color="brr[$1]}}' ../circos.chromosome.txt ../circos.link.block.txt > ../circos.link.block.color.txt
Chr1    6265416 7292507 Chr1    27686992    28717416    color=142,212,202
Chr1    8845231 9469634 Chr1    25524718    26389916    color=142,212,202
Chr1    5896416 6179289 Chr1    27217477    27603317    color=142,212,202

物种间做共线性
原理与物种内类似，就是一个建库，用另一个比对，然后将二者的处理过的gff文件（不是gff3）合并操作即可，不过为了展示在circos图上，有一点不太一样的是，你需要为每一个物种准备上面的所有东西，相信你应该可以触类旁通，下面是我简单做了四个物种间的共线性展示。

7. 文字文件，格式如下：

第1列：染色体ID
第2列：在染色体上的起始位置
第3列：在染色体上的终止位置
第4列：基因ID或者其他任何文字均可
hs1 100425066 100487997 DBT
hs1 10381671 10402787 PGD

代码，这种应该非常简单，你应该能搞定吧

# 首先，这种情况就是你在有一堆基因列表时，想获取基因区间的时候用到
# 比如说我先在有以下的基因名(必须保证gff3文件中有你的基因名)
head tair10.gene.list.some.txt
AT1G01010
AT2G01008
AT3G01010
AT4G00005
AT5G01010
# 从gff3中提取基因列表
gff3=tair10.gff3
awk -vFS="\t" -vOFS="\t" '{if($3=="gene"){match($9,/ID=([^;]+)/,a);sub(/ID=/,"",a[0]);print $1,$4,$5,a[0]}}' ${gff3} > ${gff3%%.*}.gene.list.all.txt 
head ${gff3%%.*}.gene.list.all.txt 
Chr1    3631    5899    AT1G01010
Chr1    5928    8737    AT1G01020
Chr1    11649   13714   AT1G01030
Chr1    23146   31227   AT1G01040
# 使用grep提取即可（注意前面一个文件（*some.txt）最后不能有空行，否则会全部提出来）
grep -f tair10.gene.list.some.txt tair10.gene.list.all.txt > circos.text.txt
# 甚至你想标记颜色的话，还可以添加颜色信息（非必须）
awk -vFS="\t" -vOFS="\t"  'NR==FNR{arr[$3]=$3;brr[$3]=$7}NR!=FNR{if(arr[$1]=$1){print $0,"color="brr[$1]}}' circos.chromosome.txt circos.text.txt > circos.text.color.txt
Chr1    3631    5899    AT1G01010   color=142,212,202
Chr2    1025    2810    AT2G01008   color=180,224,101
Chr3    4342    4818    AT3G01010   color=247,160,151
Chr4    1180    1536    AT4G00005   color=131,178,210
Chr5    1223    5061    AT5G01010   color=253,130,115

四、写在最后

$\color{red}{遗漏之处，欢迎私信，评论，或加群提问。}$

最后编辑于：2022.06.02 18:52:05

人面猴
序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
沈念sama阅读 194,524评论 5赞 460
死咒
序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
沈念sama阅读 81,869评论 2赞 371
救了他两次的神仙让他今天三更去死
文/潘晓璐我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
开封第一讲书人阅读 141,813评论 0赞 320
道士缉凶录：失踪的卖姜人
文/不坏的土叔我叫张陵，是天一观的道长。经常有香客问我，道长，这世上最难降的妖魔是什么？我笑而不...
开封第一讲书人阅读 52,210评论 1赞 263
港岛之恋（遗憾婚礼）
正文为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
茶点故事阅读 61,085评论 4赞 355
恶毒庶女顶嫁案：这布局不是一般人想出来的
文/花漫我一把揭开白布。她就那样静静地躺着，像睡着了一般。火红的嫁衣衬着肌肤如雪。梳的纹丝不乱的头发上，一...
开封第一讲书人阅读 46,117评论 1赞 272
城市分裂传说
那天，我揣着相机与录音，去河边找鬼。笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。我是一名探鬼主播，决...
沈念sama阅读 36,533评论 3赞 381
双鸳鸯连环套：你想象不到人心有多黑
文/苍兰香墨我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
开封第一讲书人阅读 35,219评论 0赞 253
万荣杀人案实录
序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
沈念sama阅读 39,487评论 1赞 290
护林员之死
正文独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋以下内容为张勋视角年9月15日...
茶点故事阅读 34,582评论 2赞 309
白月光启示录
正文我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
茶点故事阅读 36,362评论 1赞 326
活死人
序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
沈念sama阅读 32,218评论 3赞 312
日本核电站爆炸内幕
正文年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
茶点故事阅读 37,589评论 3赞 299
男人毒药：我在死后第九天来索命
文/蒙蒙一、第九天我趴在偏房一处隐蔽的房顶上张望。院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
开封第一讲书人阅读 28,899评论 0赞 17
一桩弑父案，背后竟有这般阴谋
文/苍兰香墨我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。一阵脚步声响...
开封第一讲书人阅读 30,176评论 1赞 250
情欲美人皮
我被黑心中介骗来泰国打工，没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。一个月前我还...
沈念sama阅读 41,503评论 2赞 341
代替公主和亲
正文我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
茶点故事阅读 40,707评论 2赞 335