一,安装circos
conda install circos
二,数据准备
1,计算染色体长度
seqtk comp genome.fa |less #可看到每一行的信息
seqtk comp genome.fa |grep "^Chr" |less
seqtk comp genome.fa |grep "^Chr" |awk '{print $1"\t"$2}' |less
seqtk comp genome.fa |grep "^Chr" |awk '{print $1"\t"$2}' > genome.len
2,生成染色体文件 7列 karyotype.txt
awk '{print "chr\t-\t"$1"\tchr"NR"\t0\t"$2"\tchr"NR}' genome.len > karyotype.txt
3, 生成窗口文件 每50k一个窗口
bedtools makewindows -w 50000 -g genome.len > genome.window.bed
4,计算每个窗口平均GC含量
先把每个窗口里的序列提取出来
计算每个窗口的GC含量
awk -F ":|-" 意思是把:或—都看作分隔符。
genome.fa里没有挂载到染色体的序列,要去掉
seqtk comp genome.fa |less #查看每一行的情况
sed -i '23,$'d genome.fa #删除从23行开始到最后一行的内容
seqtk subseq genome.fa genome.window.bed > genome.window.fa
seqtk comp genome.window.fa |awk '{print $1 "\t" ($4+$5)/($3+$4+$5+$6)}' |awk -F ":|-" '{print $1"\t"$2"\t"$3"\t"$4}' > iyun_gc.txt
5,计算每个窗口基因条数
注释文件gff转bed文件之后,记得提取挂载到染色体上的序列。
grep "^Hic_asm" iyun.bed > iyun_chr.bed
bedtools coverage -a genome.window.bed -b iyun_chr.bed -c -F 0.5 > iyun_genecount.txt
bedtools coverage -a genome.window.bed -b repeat.gff -c -F 0.5 > iyun_genecount.txt
6,计算每个窗口重复序列含量
bedtools coverage -a genome.window.bed -b repeat.gff |awk '{print $1 "\t"$2"\t"$3"\t"$7}' > iyun_repeat.txt
7,生成共线性link文件
8,LTR转录因子文件,基因个数
显示匹配到转录因子Gypsy的行
sed -n '/Gypsy/'p All_Repeat_without_trf.gff > Gypsy.gff
bedtools coverage -a genome.window.bed -b Gypsy.gff -c -F 0.5 > Gypsy.txt
9. 取对数
awk '{print $1"\t"$2"\t"$3"\t"log($4)/log(10)}' iyun_genecount.txt | sed 's/-inf/0/g' > iyun_genecount_log10.txt
三,运行
circos -conf circos.conf