1.安装
1.1 依赖软件安装
依照github上的说明需要以下软件:
需要的python包:
1.2本体安装
MCScan其实就是python的一个包(jcvi),安装只需要用pip就可以。
pip install jcvi
MCScan是支持python3的,且马上要不支持python2了,所以最好使用python3。
如果安装失败可以尝试使用conda安装python,再使用pip安装。
或下载源代码安装:
cd ~/code # or any directory of your choice
git clone git://github.com/tanghaibao/jcvi.git
pip install -e .
示例文件运行
其实wiki上这部分写的很详细了,可以直接去上面看。
https://github.com/tanghaibao/jcvi/wiki/MCscan-(Python-version)
2.1 准备数据
MCScan需要的源文件是 cds序列文件(fasta格式)和gff3文件。且最好在phytozome上下载,用这个包自带的命令只能提phytozome数据库的数据。
示例使用的是葡萄,桃子吧(Vvinifera,Ppersica)。在phytozome上下载。需要账户,注册一个就好了。
python -m jcvi.apps.fetch phytozome
Phytozome Login: xxxxxxxx
Phytozome Password:
python -m jcvi.apps.fetch phytozome Vvinifera,Ppersica
这样就有了两个cds.fa和gff3文件(压缩)。
利用gff3文件生成.bed文件:
python -m jcvi.formats.gff bed --type=mRNA --key=Name Vvinifera_145_Genoscope.12X.gene.gff3.gz -o grape.bed
python -m jcvi.formats.gff bed --type=mRNA --key=Name Ppersica_298_v2.1.gene.gff3.gz -o peach.bed
生成cds文件:
python -m jcvi.formats.fasta format Vvinifera_145_Genoscope.12X.cds.fa.gz grape.cds
python -m jcvi.formats.fasta format Ppersica_298_v2.1.cds.fa.gz peach.cds
ls *.???
grape.bed grape.cds peach.bed peach.cds
这样就有了.cds文件 和.bed文件,后续共线性分析主要用这四个文件。
2.2 共线性分析
python -m jcvi.compara.catalog ortholog grape peach
这样就生成了如下几个文件。
ls grape.peach.*
grape.peach.anchors grape.peach.last.filtered grape.peach.pdf
grape.peach.last grape.peach.lifted.anchors
grape.peach.pdf:
2.3 绘图
在得到grape.peach.anchors 后就可以进行之后的绘图。
需要准备3个文件
第1个是 .simple 文件,由.anchors 生成。
python -m jcvi.compara.synteny screen --minspan=30 --simple grape.peach.anchors grape.peach.anchors.new
第2个是 seqids 文件。里面是染色体的名字:
chr1,chr2,chr3,chr4,chr5,chr6,chr7,chr8,chr9,chr10,chr11,chr12,chr13,chr14,chr15,chr16,chr17,chr18,chr19
Pp01,Pp02,Pp03,Pp04,Pp05,Pp06,Pp07,Pp08
注意染色体名字必须在之前的bed文件中存在。
第3个是 layout 文件。里面是绘图的相关设定。
# y, xstart, xend, rotation, color, label, va, bed
.6, .1, .8, 0, , Grape, top, grape.bed
.4, .1, .8, 0, , Peach, top, peach.bed
# edges
e, 0, 1, grape.peach.anchors.simple
准备好这3个文件就可以绘图了:
python -m jcvi.graphics.karyotype seqids layout
生成了一个新的pdf文件:
这样简单地使用示例文件进行共线性分析及绘图就完成了。
scripts
get_gene_bed.pl
use Getopt::Long;
my %opts;
use Data::Dumper;
GetOptions( \%opts, "in1=s","out=s","h");
if(!defined($opts{in1}) || !defined($opts{out})||defined($opts{h})){
$USAGE;
}
open (IN1,"$opts{in1}")or die "cannot open the file\n";
open (OUT,">$opts{out}")or die "cannot open the file\n";
while(<IN1>){
chomp;
my @a = split /\t/,$_;
if($a[2] eq "gene"){
#if($a[2] eq "mRNA"){
$a = ~m/ID=([^;]*)/;
$id =$1;
#$id=~s/\.TAIR10//;
print OUT "$a[0]\t$a[3]\t$a[4]\t$id\t$a[7]\t$a[6]\n";
}
}
get_fa_bed.pl
use Math::BigFloat;
use Bio::SeqIO;
use Bio::Seq;
#读入蛋白质序列
$in=Bio::SeqIO->new(
-file => "$ARGV[1]",
-format => 'Fasta'
);
#输出蛋白序列
$out = Bio::SeqIO->new(
-file => ">$ARGV[2]",
-format => 'Fasta'
);
##提取基因ID
my %keep=();
open IN,"$ARGV[0]"or die "cannot open the file";
while (<IN>){
chomp;
next if /^#/;
my @a=split /\t/;
###下一行适当更改
$keep{"$a[3].1"}=1;
}
close IN;
while(my $seq=$in->next_seq()){
my ($id,$sequence,$desc)=($seq->id, $seq->seq, $seq->desc);
if(exists $keep{$id}){
$out->write_seq($seq);
}
}
$in->close();
$out->close();