其实滑窗的用法离你很近，比如说circle图、定位分析。滑窗按区域来统计，非nt。
先自己科普什么是CNAs和CNA？！
现在有一个需求，假如在重测序中，需要统计样本1 的bam(比对文件)里哪些区域内的测到的reads比较多？
其实原需求是统计拷贝数变异，但是拷贝数变异在二代测序中就是统计测序reads的深度信息。
测序得到的比较多的reads，就容易得到两种结果，一种是该区域随机测的比较多，另外一种也是我们想要的，那就是这个区域是结构变异中的拷贝数变异。
（如果在RNA-seq中，测序的是cDNA，那么测序得到的reads的区域越多，那就是该区域的基因表达的高，可是在重测序中，没有基因的位置这一统计单位，就算我们从gtf中找到了，这样的统计也是没有实际意义的，因为我们找的是拷贝数变异的区域，这个区域不一定必须以基因为统计单位，也可能是内含子/基因上某一段重复，也可能是基因的某一些区域重复）
错误想法
有的人会说，那把每一个碱基的测序次数，也就是深度信息统计出来，然后用（位置信息~测序深度）画图，肯定能找到测序相对较多的区域。实际上就是这么回事，如果1-100bp的碱基上每一个测序到了20X，100-200bp的碱基每一个测序到了10X，那么画图后，很明显的就能出现一组碱基在二楼（纵坐标为20），一组碱基在一楼（纵坐标为10），我们可以很好的处理数据，来得到那个二楼的区域。
但是实际测序中，并不是这样，你画出来的深度图的y轴的点必然忽上忽下，看看igv的覆盖度就能猜到。最起码你的图也是一个波动性的图，也就说明用一个碱基一个碱基的深度来统计拷贝区域是不现实的，因为你要考虑到测序reads所带来的碱基深度的波动性。
错误想法带来的提示
虽然上面的想法失败了，但是我们也能明白，我们要找的并不是一个碱基突然的覆盖率高，即测到的深度高，我们需要的是一个区域。
这样好得到一段序列的拷贝数异常，异常有重复，也有缺失。
也就是说，大部分深度测序均为1X，有的为2X、有的是0.5X（相信大家应该看过变异检测里CNV深度的那张图，或者自己项目结题报告内的覆盖度的那张图）
滑窗
因此我们可以使用滑窗的方法，滑完还可以再连接区域，就可以组成更大的区域。
滑窗也是一个算法，当然还有其他算法可以用。要么就是说，挨踢最难学的就是算法。每一个软件都有一套自己的算法~
虽然我觉得你可能看不懂，但是你还是硬着头皮看一下，其实就那么回事。接下来，来了解滑窗~
看下图出个题目：（手动画的，有点丑）
题目：假如第一条黑线是基因组的chr1染色体的1-20K区域，第二条黑线是基因组的chr2染色体的1-25K区域，橙色为覆盖度，点越高，代表所对应碱基测序到的reads条数比较多。
思考：看图，我们把橙色区域进行划分，按照中间的断点分为3个区域，
- 看区域1，只有某一段覆盖度高；
- 看区域2，整体较高；
- 看区域3，整体较低；
我们再假如:
- 区域1为大部分数据的正常深度，为1X；
- 区域2的深度为1.5X或者2X；区域3为0.5X。那么区域2下面箭头的区域就是拷贝数变异的异倍体中的多倍体；
- 区域3就是拷贝数变异中的缺失（具体是多位置的某多拷贝位置缺失还是本身这就是重复区域的某一段重复缺失等，不要深究，深究就给自己挖坑，以后慢慢学吧）。
窗口解决思路
窗口思路解决：比如我们从1bp-10k bp之间为一个bin，也就是一个窗口，那么我能统计出来这个窗口的reads的比对的每一个碱基的深度：
接着求出10k+1 bp到20k bp间所有碱基的深度和，也用公式来做，定能找到哪个平均深度数值比较高的。
疑问：这个时候，你会有一个疑问，万一某个区域，深度确实很高，但是被定义的bin拆给了两个区域，平均后，反而找不到那个高深度的测序区域了。
滑窗：所以我们可以用10k为一个窗口（即bin），5k为一个步移，滑动窗口的来统计，那么1-10k为一个bin，5k-15k...；同理我们还可以以5K或者20K为一个bin，步移也可以自定义。当我们得到了我们的区域，肯定间隔都是bin长度的，接下来我们还可以合并我们找到的异倍体区域。
滑窗弊端
弊端：我们只能按照步移的单位，来统计CNAs的区域，这样的单位都是k，我们并不能定位到某一个位置的碱基。
除非你的步移为1bp，但是也应该很慢，或者想到一个更好的算法，或者找到一个算法/软件。
但是：如果你要的知识基因的CNAs，那么就没必要精确到碱基，有一个大概的区域，把区域的基因名提取出来就可以，即使某基因部分区域落在该区域。
必须掌握
- 什么是异倍体，非单倍体，可能统计出倍数的是1.2X；1.5X；2X都是有可能的。因为我们算了平均区域，再加上测序的随机性。并不能完全测序到二倍体就必须是2X。
- 滑窗的bin是什么，步移是什么？
能用软件做同样的需求，用软件；
能用R包的函数，就用函数！
网友评论