由于忙課題的原因国章，斷更了好久，雖然斷更的時(shí)間內(nèi)學(xué)習(xí)到了很多東西豆村，但細(xì)想一下還是記錄在自己博客里比較好液兽，方便自己查閱也方便其他人學(xué)習(xí)
JCVI這個(gè)軟件網(wǎng)上有太多教程，再記錄一次是想按照自己平時(shí)的使用習(xí)慣把整個(gè)流程從頭到尾捋一遍
安裝
能conda就conda掌动，就一個(gè)字四啰，省心！
conda create -n jcvi
conda install jcvi
準(zhǔn)備文件
分析需要用到 最長的轉(zhuǎn)錄本粗恢，很多軟件分析時(shí)其實(shí)已經(jīng)內(nèi)置了相應(yīng)的程序柑晒，可以使用以下的方法來提取
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文件準(zhǔn)備好
姑且稱之為A.bed, A.pep, B.bed, B.pep 吧
分析
python -m jcvi.compara.catalog ortholog A B --dbtype prot --cscore=.99
我用的是蛋白序列文件，所以加個(gè) --dbtype prot眷射，cds序列則不需要
1匙赞，先分析整體的共線性（可視化）
1.1
對(duì)上面分析的結(jié)果進(jìn)行一個(gè)簡單過濾
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個(gè)基因的block將會(huì)移除
1.2
創(chuàng)建seqid文件佛掖，想展示哪些染色體就將其寫入到這個(gè)文件中
Chr1A,Chr2A,Chr3A,Chr4A,Chr5A
Chr1B,Chr2B,Chr3B,Chr4B,Chr5B
1.3
創(chuàng)建layout文件，用于設(shè)置繪圖參數(shù)涌庭，最終出圖都是通過調(diào)整這個(gè)文件來達(dá)到自己想要的效果芥被。
# 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為旋轉(zhuǎn)角度拴魄，比如45或者-45啥的，看自己需求席镀，label就是基因組名字了匹中，后面的top表示名字在染色體上方表示，也可以設(shè)置left愉昆，right职员，bottom麻蹋。
edges那里就是連線跛溉，0是第一個(gè)基因組，1是第二個(gè)扮授，然后就是3芳室，4...,想讓哪兩個(gè)連線就設(shè)置好對(duì)應(yīng)的edges和simple文件
最后運(yùn)行命令
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基因組作為標(biāo)準(zhǔn)堪侯，從B中得到對(duì)應(yīng)的共線性基因，假如B中有的基因而A中沒有的話是無法知道的荔仁，--iter=1表示基因是一對(duì)一伍宦，參數(shù)調(diào)整范圍為1-100
如果有更多基因組的話，可以繼續(xù)按照上面的流程得到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
準(zhǔn)備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
三個(gè)基因組的話差不多就這樣乏梁，跟上面的差不多次洼，不做過多解釋了，這些數(shù)字是我亂寫的遇骑，畫的時(shí)候多試試就知道了
最后卖毁！畫圖！
python -m jcvi.graphics.synteny all.i1.blocks all.bed layout
寫在最后
以上就是全部內(nèi)容了落萎，除了沒放結(jié)果展示亥啦，基本上分析下來是沒有問題的，當(dāng)然主要是寫給我看的练链，所以不會(huì)像其他博主一樣面面俱到
更新一下內(nèi)容
之前一直用的舊版本的JCVI翔脱，今天才發(fā)現(xiàn)有更新
如果不想在共線性圖上展示染色體編號(hào)的話可以加一個(gè)--nocircles
python -m jcvi.graphics.karyotype seqids layout.2 --nocircle
如果不想用曲線，可以加--shadestyle=line
最后展示下多基因組的共線性圖
再次更新
展示整個(gè)基因組內(nèi)感興趣基因的共線性
近期想要在兩個(gè)小麥基因組間看一下感興趣基因的共線性媒鼓，并需要把相應(yīng)的基因標(biāo)識(shí)出來碍侦，有需求就有更新粱坤，這個(gè)操作跟全基因組展示基因家族的共線性是一個(gè)道理，做基因家族的朋友也可以參考一下瓷产。要畫這個(gè)圖其實(shí)很簡單站玄，跟 1.3 的流程是一樣的，關(guān)鍵是獲得xxx.simple文件濒旦，然后才能開始下面的繪圖株旷，現(xiàn)在需要獲得兩個(gè)基因組感興趣基因的1對(duì)1文件，那就需要走一遍 2尔邓，局部共線性分析（可視化） 這個(gè)流程晾剖，得到i1.blocks后，根據(jù)感興趣基因的列表來匹配1對(duì)1blocks中的基因梯嗽，很簡單齿尽，寫個(gè)程序就可以：
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
然后根據(jù)這個(gè)新blcoks弄成一個(gè)simple文件，可以使用以下代碼：
cat interestedGene.i1.blocks|awk -vOFS="\t" '{print"#FF0000*"$1,$1,$2,$2,"10","+"}' > interestedGene.simple
如果想在circos上展示共線性灯节，那就需要弄一個(gè)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文件里找基因位置全跨，然后把這一段的截出來就行了
更新這一段內(nèi)容是因?yàn)榍岸螘r(shí)間看了發(fā)表在MP上的Genetribe軟件文章，這篇文章使用優(yōu)化的算法鑒定不同基因組之前的共線性基因亿遂，并將這些基因分為1對(duì)1最佳（RBH）基因浓若，單向最佳匹配基因（SBH），以及1對(duì)多基因等蛇数，最近在用這個(gè)軟件分析自己的課題挪钓，最后的可視化用jcvi的程序就可以做到，感興趣的可以看看https://doi.org/10.1016/j.molp.2020.09.019
