共线性分析

安装jcvi

conda install jcvi 准备文件gene.gff3 cds.fa protein.fa

coline analysis

参考： https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version) https://www.cnblogs.com/zhanmaomao/p/12525411.html

conda activate mmdetection
conda install jcvi
conda install seqkit

安装lastal

网址：http://last.cbrc.jp
unzip last-1060.zip
cd last-1060
make

从Phytozome数据库获取物种数据（如果远端server无法联网，请直接下载基因组的相关数据，跳过此步）

python -m jcvi.apps.fetch phytozome Phytozome Login: xxxxxxxx Phytozome Password: 输入你的Phytozome用户名和密码。 python -m jcvi.apps.fetch phytozome Athaliana Phytozome的速度实在太慢，还是从本地用工具下载吧。

从gff3获取bed

以Spo，和Fve为例子

python -m jcvi.formats.gff bed --type=mRNA --key=ID Spo.genome.gff3 -o Spo.bed
python -m jcvi.formats.gff bed --type=mRNA --key=ID Fragaria_vesca_v4.0.a1.transcripts.gff3.gz -o Fve.bed

参数

--type：gff文件中第三列 --key：type对应的第九列信息前缀

我们分析只需要用到每个基因最长的转录本就行，在Fve中是以多个转录本进行存储，因为需要获取最长转录本

将Fve中bed进行去重复

python -m jcvi.formats.bed uniq Spo.bed
python -m jcvi.formats.bed uniq Fve.bed

cds和pep序列均可以进行共线性分析

conda install seqkit

根据上述得到的.bed文件调取对应cds和蛋白序列

# Spo
seqkit grep -f <(cut -f4 Spo.bed) Chr_genome_final.cds.fa  | seqkit seq -i >Spo.cds
seqkit grep -f <(cut -f4 Spo.bed) Spo.genome.protein.fa  | seqkit seq -i >Spo.pep

# Fve
seqkit grep -f <(cut -f4 Fve.bed) Fragaria_vesca_v4.0.a1_makerStandard_CDS.fasta.gz| seqkit seq -i >Fve.cds
seqkit grep -f <(cut -f4 Fve.bed) Fragaria_vesca_v4.0.a1_makerStandard_proteins.fasta.gz  | seqkit seq -i >Fve.pep

共线性分析

新建一个文件夹，方便在报错的时候，把全部都给删了

mv Fve.uniq.bed Fve.bed
mv Spo.uniq.bed Spo.bed
mkdir Spo_Fve && cd Spo_Fve
ln -s ../Fve.cds ./
ln -s ../Fve.bed ./
ln -s ../Fve.pep ./
ln -s ../Spo.cds ./
ln -s ../Spo.bed ./
ln -s ../Spo.pep ./
## 运行代码
python -m jcvi.compara.catalog ortholog (--dbtype prot[蛋白分析]) --no_strip_names Spo Fve
python -m jcvi.compara.catalog ortholog --dbtype prot --no_strip_names Spo Fve

结果： Spo.Fve.anchors：共线性block区块（高质量） Spo.Fve.last：原始的last输出 Spo.Fve.last.filtered：过滤后的last输出 Spo.Fve.pdf：点阵图

如果遇到报错，多半是要安装python包，更新Latex. 使用conda安装的latex果然有问题，自己动手安装latex参考。文件3G,安装比较慢。

#安装latex的目录是/share/softwares/Latex
export MANPATH=$MANPATH:/share/softwares/Latex/texmf-dist/doc/man
export INFOPATH=$INFOPATH:/share/softwares/Latex/texmf-dist/doc/info
export PATH=$PATH:/share/softwares/Latex/bin/x86_64-linux

~特殊情况的处理，当组装级别是scaffold时候，处理办法。~ 在设置seqids时候，面临Atr组装的是scaffold，是否需要展示全部的scaffold。文献中没有提及，但是都展示了部分的scaffold。此处是设置展示前120个scaffold。

cat Atr.bed|cut -f1|sort |uniq |head -120 |rev|cut -d " " -f1|rev >Atr.id
cat Spo.bed|cut -f1|sort|uniq|grep -v utg|rev|cut -d " " -f1|rev >Spo.id
##转置数据从列到行
cat Spo.id|awk 'BEGIN{c=0;} {for(i=1;i<=NF;i++) {num[c,i] = $i;} c++;} END{ for(i=1;i<=NF;i++){str=""; for(j=0;j<NR;j++){ if(j>0){str = str","} str= str""num[j,i]}printf("%s\n", str)} }' >Spo.ids
cat Atr.id|awk 'BEGIN{c=0;} {for(i=1;i<=NF;i++) {num[c,i] = $i;} c++;} END{ for(i=1;i<=NF;i++){str=""; for(j=0;j<NR;j++){ if(j>0){str = str","} str= str""num[j,i]}printf("%s\n", str)} }' >Atr.ids
cat Spo.ids Atr.ids >seqids

可视化

首先生成.sinple文件

python -m jcvi.compara.synteny screen --minspan=30 --simple Spo.Fve.anchors Spo.Fve.anchors.new

编辑配置文件seqids 和layout

注意：seqids上下的顺序和layout的顺序是一一对应的。而且seqids后面不要加注释，不然不会显示最后一条染色体。Chr1和Chr01格式一定要和bed里一样，不然也不会显示。 layout文件一定要按照格式写。不要多空格或tab，不然会报错。

设置需要展示染色体号

vi seqids
Chr01,Chr02,Chr03,Chr04,Chr05,Chr06,Chr07,Chr08,Chr09,Chr10,Chr11,Chr12,Chr13,Chr14,Chr15,Chr16,Chr17
Chr01,Chr02,Chr03,Chr04,Chr05,Chr06,Chr07,Chr08,Chr09,Chr10,Chr11,Chr12,Chr13,Chr14,Chr15,Chr16,Chr17

# 设置颜色，长宽等
vi layout
# y, xstart, xend, rotation, color, label, va, bed
 .6,    .1,    .8,    0,    red,    Spo,    top,     Spo.bed
 .4,    .1,    .8,    0,    blue,    Fve,    top,    Fve.bed
# edges
e, 0, 1, Spo.Fve.anchors.simple

注意, #edges下的每一行开头都不能有空格

运行代码

python -m jcvi.graphics.karyotype seqids layout

突出显示某一条，前面加上g*

vi Spo.Fve.anchors.simple
g*Spo03717      Spo03716        Bv3_048620_odzi.t1      Bv3_049090_cxmm.t1      46      +
Spo17356        Spo17350        Bv1_001140_tedw.t1      Bv1_001540_xzdn.t1      41      -
Spo13685        Spo13730        Bv5_123480_yfcy.t1      Bv5_123900_rucq.t1      46      -
Spo26250        Spo26280        Bv5_117270_qhwj.t1      Bv5_117680_iykf.t1      36      +
Spo19005        Spo06982        Bv7_173830_frmo.t1      Bv7_174150_pckw.t1      37      +
Spo19374        Spo20559        Bv4_081440_riqq.t1      Bv4_081750_yeuy.t1      32      -

运行

python -m jcvi.graphics.karyotype seqids layout

三个物种之间比较

常规布局

 vi layout
# y, xstart, xend, rotation, color, label, va,  bed
 .7,     .1,    .8,      15,      blue, F.vesca, top, Fve.bed
 .5,     .1,    .8,       0,      red, S.pohuashanensis, top, Spo.bed
 .3,     .1,    .8,     -15,      green, M.domestica, bottom, Mdo.bed
# edges
e, 0, 1, Spo.Fve.anchors.simple
e, 1, 2, Spo.Mdo.anchors.simple

vi seqids
Fvb1,Fvb2,Fvb3,Fvb4,Fvb5,Fvb6,Fvb7
Chr01,Chr02,Chr03,Chr04,Chr05,Chr06,Chr07,Chr08,Chr09,Chr10,Chr11,Chr12,Chr13,Chr14,Chr15,Chr16,Chr17
Chr01,Chr02,Chr03,Chr04,Chr05,Chr06,Chr07,Chr08,Chr09,Chr10,Chr11,Chr12,Chr13,Chr14,Chr15,Chr16,Chr17

python -m jcvi.graphics.karyotype seqids layout
##--nocircles 不画出圆圈的染色体底纹
python -m jcvi.graphics.karyotype seqids normal_layout --nocircles

三角形布局

# y, xstart, xend, rotation, color, label, va,  bed
 .5,      0.025,    0.625,      60,      , Grape, top, grape.bed
 .2,      0.2,    .8,       0,      , Peach, top, peach.bed
 .5,     0.375,    0.975,     -60,      , Cacao, top, cacao.bed
# edges
e, 0, 1, grape.peach.anchors.simple
e, 1, 2, peach.cacao.anchors.simple

运行

python -m jcvi.graphics.karyotype seqids layout

注意：anchors.simple是同义基因块（synteny blocks）.

每行是一个块（synteny blocks），每行前两列是物种1的起始基因和中止基因，3-4列是物种2的起始基因和中止基因,5列是得分，6列是正反链。

head grape.peach.anchors.simple
gra01G000550.1   gra01G001220.1   pea01G1001000    pea01G1008900    57      -
gra01G001300.1   gra01G011750.1   pea01G1009700    pea01G1130400    1025    +

所以下面这样取首尾共有基因作为代表，显示绿色并没有生物学意义。 ~~如果想让Peach的基因既在Grape又在Cacao上的基因块线条显示绿色（green）.~~

cat grape.peach.anchors.simple|awk '{print $1"\n"$2}' |sort -u >peach.txt
cat grape.cacao.anchors.simple|awk '{print $1"\n"$2}' |sort -u >cacao.txt
join peach.txt cacao.txt >gene.list

add_g.sh内容如下

simple=$1
cat gene.list| while read line
do
#备份原始simple文件
cp $simple ${simple}.origin
#先获取匹配行，再在匹配行首添加`g*`,因为可能存在一行添加两次的情况。
eval sed -i \'"/${line}/ s/^/g\*/"\' $simple
#替换行首添加两次的`g*g`为`g`
sed -i "s/^g\*g/g/" $simple
done

替换simple文件中对应基因的行首

bash add_g.sh grape.peach.anchors.simple
bash add_g.sh grape.peach.anchors.simple
#重新绘图
python -m jcvi.graphics.karyotype seqids layout

~~无意义代码结束~~ 经过咨询jcvi作者tanghaibao,可以使用.anchors构建.simple文件，把每个基因都作为一个block这样可以实现针对基因水平表现共线性。.simple需要6列，使用.anchors的第1，1，2，2，3，+组成新的.simple文件，后续和原来一样出图即可。

方法2：使用MCScanX分析和VGSC(可视化)

参考：https://www.jianshu.com/p/50f79e4d7460 VGSC地址直接下载文件，jar可以直接运行。

下载MCScanX

axel -n 16 http://chibba.pgml.uga.edu/mcscan2/MCScanX.zip
unzip MCScanX.zip

因为MCScanX不支持64位系统，需要运行下面三行命令，分别在三个文件的第一行前添加`#include <unistd.h>`

sed -i '1i\#include <unistd.h>' msa.h
sed -i '1i\#include <unistd.h>' dissect_multiple_alignment.h
sed -i '1i\#include <unistd.h>' detect_collinear_tandem_arrays.h
make -j4
./MCScanX -h

MCScanX也可以计算KaKs,但是需要安装clustalw conda install -c bioconda clustalw

构建gff文件和blast文件

gff文件格式

#head
 at3  AT3G19630  6818676  6820674
 at5  AT5G11220  3577057  3577854
 at2  AT2G29110  12506880  12510552
#tail
 vv1  GSVIVT01012140001  1098524  1102524
 vv18  GSVIVT01034900001  16234068  16271170
 vv8  GSVIVT01025554001  14085153  14089377

直接使用上面生成的bed文件，grep的目的是只提取chr的数据，scaffold的舍弃。参考链接：https://www.jianshu.com/p/50f79e4d7460

#根据你的基因组是否有chr来决定，没有chr用第二行，有chr用第一行。
sed 's/Chr/Spo/g' Spo.bed|grep Spo >Spo.gff
sed 's/^/Ath/g' Ath.bed |grep -v Mt|grep -v Pt >Ath.gff 
cat Ath.gff Spo.gff  >Ath_Spo.gff3
cat Ath_Spo.gff3|awk '{print $1"\t"$4"\t"$2"\t"$3}' >Ath_Spo.gff

准备blast文件

#蛋白文件来自于jcvi生成的
cat Ath.pep Spo.pep >Ath_Spo.pep
##使用blast+版本
makeblastdb -in Ath_Spo.pep -dbtype prot -out db_Ath_Spo -parse_seqids
blastp -query Ath_Spo.pep -db db_Ath_Spo -out Ath_Spo.blast -evalue 1e-5 -outfmt 6 -num_alignments 5  -num_threads 32

参数说明： -dbtype 后接序列类型，nucl为核酸，prot为蛋白 -parse_seqids 推荐加上 -out 后接数据库名 -logfile 日志文件，如果没有默认输出到屏幕 (-logfile xxx.log)

MCScanX运行

mkdir Ath_Spo
cp Ath_Spo.gff ./Ath_Spo
cp Ath_Spo.blast ./Ath_Spo
MCScanX ./Ath_Spo -m 5
#-m 5指定最大gap是5，默认是25.一般不用指定。

生成Ath_Spo.collinearity

运行VGSC可视化

如果java版不能使用，可以使用在线版

java -jar VGSC.jar -tp Bar -ig Ath_Spo.gff -is Ath_Spo.collinearity -ic Ath_Spo.ctl
-os Ath_Spo.svg

参数说明： -tp Bar|Circle| Dot |DS -ig gff -is os_sb.collinearity -ic bar.ctl |circle.ctl -os outputfile bar.ctl是出图配置文件，安装包里有对应的模板 Ath_Spo.ctl文件内容如下：

800     //dimension (in pixels) of x axis
800     //dimension (in pixels) of y axis
Ath1,Ath2,Ath3,Ath4,Ath5       //reference chromAthomes
Spo01,Spo02,Spo03,Spo04,Spo05,Spo06,Spo07,Spo08,Spo09,Spo10,Spo11,Spo12,Spo13,Spo14,Spo15,Spo16,Spo17    //target chromAthomes

共线性可视化可以使用基于ggplot2的`ggalluvial`

https://github.com/corybrunson/ggalluvial

共线性分析

安装jcvi

coline analysis

安装lastal

从Phytozome数据库获取物种数据（如果远端server无法联网，请直接下载基因组的相关数据，跳过此步）

从gff3获取bed

以Spo，和Fve为例子

参数

将Fve中bed进行去重复

cds和pep序列均可以进行共线性分析

根据上述得到的.bed文件调取对应cds和蛋白序列

共线性分析

新建一个文件夹，方便在报错的时候，把全部都给删了

可视化

首先生成.sinple文件

编辑配置文件seqids 和layout

设置需要展示染色体号

运行代码

突出显示某一条，前面加上g*

运行

三个物种之间比较

常规布局

三角形布局

运行

注意：anchors.simple是同义基因块（synteny blocks）.

方法2：使用MCScanX分析和VGSC(可视化)

下载MCScanX

因为MCScanX不支持64位系统，需要运行下面三行命令，分别在三个文件的第一行前添加#include <unistd.h>

构建gff文件和blast文件

准备blast文件

MCScanX运行

运行VGSC可视化

共线性可视化可以使用基于ggplot2的ggalluvial

因为MCScanX不支持64位系统，需要运行下面三行命令，分别在三个文件的第一行前添加`#include <unistd.h>`

共线性可视化可以使用基于ggplot2的`ggalluvial`