不多bb,直接上小脚本。
本分析流程的数据来源于《Allele-aware chromosome-scale assembly of the allopolyploid genome of hexaploid Ma bamboo (Dendrocalamus latiflorus Munro)》
(1)配置config
每一行内为一组需要进行共线性分析的基因组区域。
vim config
A B
B C
(2)写个循环
cat config | while read id
do
arr=($id)
name1=${arr[0]}
name2=${arr[1]}
dir=${name1}_${name2}
mkdir ${dir}
cp ${name1}.uniq.bed ${dir}/${name1}.bed
cp ${name2}.uniq.bed ${dir}/${name2}.bed
cp ${name1}.cds ${dir}/
cp ${name2}.cds ${dir}/
cd ${dir}
python -m jcvi.compara.catalog ortholog --no_strip_names ${name1} ${name2} --cscore=.99
python -m jcvi.compara.synteny screen --minspan=30 --simple ${name1}.${name2}.anchors ${name1}.${name2}.anchors.new
tmpstr1=`less -S ${name1}.bed | awk '{print $1}' | uniq | tr '\n' ','`
str1=${tmpstr1%?}
echo $str1 > seqids
tmpstr2=`less -S ${name2}.bed | awk '{print $1}' | uniq | tr '\n' ',' `
str2=${tmpstr2%?}
echo $str2 >> seqids
content="
# y, xstart, xend, rotation, color, label, va, bed\n .6, .1, .8, 0, , ${name1#*_}, top, ${name1}.bed\n .4, .1, .8, 0, , ${name2#*_}, top, ${name2}.bed\n# edges\ne, 0, 1, ${name1}.${name2}.anchors.simple
"
echo -e $content > layout
python -m jcvi.graphics.karyotype seqids layout
cd ../
done
上述代码已经自行绘制了基因组两两区域之间的宏观共线性图。
(3)设置更高级一点的画布
用下面给出的画布,就能近似的绘制出文章中的图。
对画布的配置进行一个简要的说明:
-
# y, xstart, xend, rotation, color, label, va, bed
这一列为track对应需要配置的注释,接这一行对每一个track进行配置 - track的索引和Python中的索引一样,从0开始
- 共线性关系(edges)所使用的索引与track的索引数一样
- layout中不能有空行
# y, xstart, xend, rotation, color, label, va, bed
.7, 0.05, 0.45, 15, , A1, top, dla_A1.bed
.4, 0.05, 0.45, 15, , B1, top, dla_B1.bed
.1, 0.05, 0.45, 15, , C1, top, dla_C1.bed
.7, 0.55, 0.95, -15, , A2, top, dla_A2.bed
.4, 0.55, 0.95, -15, , B2, top, dla_B2.bed
.1, 0.55, 0.95, -15, , C2, top, dla_C2.bed
.9, 0.4, 0.8, 0, , O.sativa, top, osa.bed
# edges
e, 0, 1, dla_A1.dla_B1.anchors.simple
e, 1, 2, dla_B1.dla_C1.anchors.simple
e, 3, 4, dla_A2.dla_B2.anchors.simple
e, 4, 5, dla_B2.dla_C2.anchors.simple
e, 0, 6, dla_A1.osa.anchors.simple
e, 3, 6, dla_A2.osa.anchors.simple
放在一起对比一下:
有一些细节还是需要接着调调,但是大体已经差不多了(明天还得考六级,没时间了,先不整。)
网友评论