基因组共线性分析-JCVI

由于忙课题的原因，断更了好久，虽然断更的时间内学习到了很多东西，但细想一下还是记录在自己博客里比较好，方便自己查阅也方便其他人学习

JCVI这个软件网上有太多教程，再记录一次是想按照自己平时的使用习惯把整个流程从头到尾捋一遍

安装

能conda就conda，就一个字，省心！

conda create -n jcvi
conda install jcvi

准备文件

分析需要用到 最长的转录本，很多软件分析时其实已经内置了相应的程序，可以使用以下的方法来提取

less your_gff.gff |grep "mRNA"|cut -d ";" -f1|cut -f1,4,5,7,9|awk -vOFS="\t" '{print$1,($2-1),$3,$5,"0",$4}'|sed 's/ID=//g' > your.bed
conda activate jcvi
python -m jcvi.formats.bed uniq your.bed
seqkit grep -f <(cut -f4 your.uniq.bed) your_pep.fa > your_pep.long.fa

利用上面的方法把需要分析的物种的bed文件以及蛋白或cds的fasta文件准备好
姑且称之为A.bed, A.pep, B.bed, B.pep 吧

分析

python -m jcvi.compara.catalog ortholog A B  --dbtype prot --cscore=.99

我用的是蛋白序列文件，所以加个 --dbtype prot，cds序列则不需要

1，先分析整体的共线性（可视化）

1.1
对上面分析的结果进行一个简单过滤

python -m jcvi.compara.synteny screen --minspan=30 --simple A.B.anchors A.B.anchors.new

--minspan: remove blocks with less span than this.跨度小于30个基因的block将会移除

1.2
创建seqid文件，想展示哪些染色体就将其写入到这个文件中

Chr1A,Chr2A,Chr3A,Chr4A,Chr5A
Chr1B,Chr2B,Chr3B,Chr4B,Chr5B

1.3
创建layout文件，用于设置绘图参数，最终出图都是通过调整这个文件来达到自己想要的效果。

# y, xstart, xend, rotation, color, label, va,  bed
 .6,     .2,    .8,       0,      , A, top, A.bed
 .4,     .2,    .8,       0,      , B, top, B.bed
# edges
e, 0, 1, aly.ath.anchors.simple

y为纵轴位置，x为横轴位置，rotation为旋转角度，比如45或者-45啥的，看自己需求，label就是基因组名字了，后面的top表示名字在染色体上方表示，也可以设置left，right，bottom。
edges那里就是连线，0是第一个基因组，1是第二个，然后就是3，4...,想让哪两个连线就设置好对应的edges和simple文件

最后运行命令

python -m jcvi.graphics.karyotype seqids layout

2，局部共线性分析（可视化）

获得block文件

python -m jcvi.compara.synteny mcscan A.bed A.B.lifted.anchors --iter=1 -o A.B.i1.blocks

解释一下的话就是，以A基因组作为标准，从B中得到对应的共线性基因，假如B中有的基因而A中没有的话是无法知道的，--iter=1表示基因是一对一，参数调整范围为1-100
如果有更多基因组的话，可以继续按照上面的流程得到blocks

python -m jcvi.compara.synteny mcscan A.bed A.C.lifted.anchors --iter=1 -o A.C.i1.blocks
#然后将这些block合并到一起
paste *.i1.blocks|cut -f1,2,4 > all.i1.blocks
#合并bed文件
cat A.bed B.bed C.bed > all.bed

准备layout文件

# x,   y, rotation,     ha,     va, color, ratio,            label
0.2, 0.3,        0,  center,    bottom,   ,    .8,         1A
0.8, 0.3,        0,  center,    bottom,   ,    .8,         1B
0.47, 0.3,        0,  center,    bottom,   ,    .8,         1C
# edges
e, 0, 1
e, 1, 2
e, 2, 3

三个基因组的话差不多就这样，跟上面的差不多，不做过多解释了，这些数字是我乱写的，画的时候多试试就知道了

最后！画图！

python -m jcvi.graphics.synteny all.i1.blocks all.bed layout

写在最后

以上就是全部内容了，除了没放结果展示，基本上分析下来是没有问题的，当然主要是写给我看的，所以不会像其他博主一样面面俱到

更新一下内容

之前一直用的旧版本的JCVI，今天才发现有更新
如果不想在共线性图上展示染色体编号的话可以加一个--nocircles

python -m jcvi.graphics.karyotype seqids layout.2 --nocircle

如果不想用曲线，可以加--shadestyle=line
最后展示下多基因组的共线性图

测试下

再次更新

展示整个基因组内感兴趣基因的共线性
近期想要在两个小麦基因组间看一下感兴趣基因的共线性，并需要把相应的基因标识出来，有需求就有更新，这个操作跟全基因组展示基因家族的共线性是一个道理，做基因家族的朋友也可以参考一下。要画这个图其实很简单，跟 1.3 的流程是一样的，关键是获得xxx.simple文件，然后才能开始下面的绘图，现在需要获得两个基因组感兴趣基因的1对1文件，那就需要走一遍 2，局部共线性分析（可视化） 这个流程，得到i1.blocks后，根据感兴趣基因的列表来匹配1对1blocks中的基因，很简单，写个程序就可以：

import sys
gene_list = open(sys.argv[1],'r')
blocks = open(sys.argv[2],'r')
list = []
for i in gene_list:
    i = i.strip()
    list.append(i)
for line in blocks:
    line = line.strip()
    item1 = line.split("\t")[0]
    item2 = line.split("\t")[1]
    if item1 in list and item2 != ".":
        print(item1+"\t"+item2)
    else:
        continue

python get_interestedGene_from_i1_blocks.py gene_list all.i1.blocks > interestedGene.i1.blocks

然后根据这个新blcoks弄成一个simple文件，可以使用以下代码：

cat interestedGene.i1.blocks|awk -vOFS="\t" '{print"#FF0000*"$1,$1,$2,$2,"10","+"}' > interestedGene.simple

如果想在circos上展示共线性，那就需要弄一个link文件，需要基因的bed文件，以及上面的interestedGene.i1.blocks，这里也把代码贴上：

from sys import argv
bed = argv[1]
syn = argv[2]
ref_dict = {line.split("\t")[3]:line.split("\t")[0:3] for line in open(bed)}
for line in open(syn):
    items = line.strip().split("\t")
    s_gene = items[0]
    e_gene = items[1]
    s_chr,s_gene_s,s_gene_e = ref_dict[s_gene][0:3]
    e_chr,e_gene_s,e_gene_e = ref_dict[e_gene][0:3]
    circos_input = [s_chr,s_gene_s,s_gene_e,e_chr,e_gene_s,e_gene_e]
    out = '\t'.join(circos_input)
    print(out)

python get_link_from_i1.blocks.py all.gene.bed interestedGene.i1.blocks > interestedGene.links

微观共线性分析就简单了，直接在i1.blocks文件前面加上颜色就可以了，想看哪一段blocks，就去bed文件里找基因位置，然后把这一段的截出来就行了

更新这一段内容是因为前段时间看了发表在MP上的Genetribe软件文章，这篇文章使用优化的算法鉴定不同基因组之前的共线性基因，并将这些基因分为1对1最佳（RBH）基因，单向最佳匹配基因（SBH），以及1对多基因等，最近在用这个软件分析自己的课题，最后的可视化用jcvi的程序就可以做到，感兴趣的可以看看https://doi.org/10.1016/j.molp.2020.09.019