一：数据下载
从http://plants.ensembl.org/index.html分别下载这两个物种的GFF文件和CDS序列.
# Athaliana
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_thaliana/cds/Arabidopsis_thaliana.TAIR10.cds.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_thaliana/Arabidopsis_thaliana.TAIR10.44.gff3.gz
# Alyrata
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/fasta/arabidopsis_lyrata/cds/Arabidopsis_lyrata.v.1.0.cds.all.fa.gz
wget ftp://ftp.ensemblgenomes.org/pub/plants/release-44/gff3/arabidopsis_lyrata/Arabidopsis_lyrata.v.1.0.44.gff3.gz
二：数据预处理
将GFF3转成BED格式
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Arabidopsis_thaliana.TAIR10.44.gff3.gz > ath.bed
python -m jcvi.formats.gff bed --type=mRNA --key=transcript_id Arabidopsis_lyrata.v.1.0.44.gff3.gz > aly.bed
然后将bed进行去重复
python -m jcvi.formats.bed uniq ath.bed
python -m jcvi.formats.bed uniq aly.bed
最后我们得到了ath.uniq.bed和osa.uniq.bed, 根据bed文件第4列就可以用于提取cds序列和蛋白序列。
seqkit grep -f <(cut -f 4 aly.bed ) Arabidopsis_lyrata.v.1.0.cds.all.fa.gz | seqkit seq -i > aly.cds
seqkit grep -f <(cut -f 4 ath.bed ) Arabidopsis_thaliana.TAIR10.cds.all.fa.gz | seqkit seq -i > ath.cds
blast比对
makeblastdb -in aly.cds -out db/aly -dbtype nucl
blastn -num_threads 20 -query ath.cds -db db/aly -outfmt 6 -evalue 1e-5 -num_alignments 5 > ath_aly.blast
用jcvi.compara.blastfilter对结果进行过滤
python -m jcvi.compara.blastfilter --no_strip_names ath_aly.blast --sbed aly.bed --qbed ath.bed
最后输出aly_ath.blast.filtered用于做图
python -m jcvi.graphics.blastplot ath_aly.blast.filtered --sbed aly.bed --qbed ath.bed
三：结果：
ath_aly
