1.划分窗口
bedtools makewindows -g Chr.length -w 50000 > 50k.windows
Chr.length就是每条染色体的长度
2.计算每个滑窗内基因的数量 #同理可以换成任何其余东西比如SNP
grep -w "gene" input.gff | awk '{print 4"\t"$5}' > gene.pos
gene.pos长这样，每个基因的位置信息, 只要前三列的信息就行，其余无所谓
bedtools intersect -a 50k.windows -b gene.pos -c > out
最后的结果和TBtools输出的一致，光拿基因密度来说
如果不需要基因密度为0的窗口的信息，还是用TBtools方便一点，后续画什么Circos图啥的
TBtools Ref: https://www.jianshu.com/p/801807865864
- 滑窗统计基因组GC含量
seqkit sliding -s 100000 -W 100000 input.fa | seqkit fx2tab -n -g > out
用TBtools输出文件看着舒服很多，还顺带有N，GC skew两个参数
https://www.jianshu.com/p/de97067136a9
3.1. 滑窗统计基因组GC含量
但是上述两个方法会有一个细节问题，比如我以50kb滑窗计算GC含量
如果最后一个窗口没有50kb这么长，seqkit会跳过这个窗口，TBtools则是会把最后一个窗口并入到前一个窗口中
所以最后改用 https://www.cnblogs.com/liujiaxin2018/p/16567643.html提供的脚本
即使最后一个窗口没有50kb，也会照常计算, 唯一的缺点是运行的相对慢一点
网友评论