【直播】我的基因组21:为什么我算出的染色体覆盖度与公司差异甚大
我们在第19讲留下了一个悬念:为什么我算的染色体覆盖度与公司给我的报告不符合!!
公司算出每条染色体的覆盖度都接近于100%了,而我是结果是:
很明显,chrY覆盖度严重偏低,21,22染色体也略堪忧。然而我对自己脚本能力及其有信心,反复检查觉得不可能写错。所以就理直气壮的写邮件咨询了给我测序的公司项目负责人~
他们反应挺迅速的,隔天就给我了答复:
曾老师:
您好!
根据与信息同事沟通,分析您提出的问题,现解答如下:
老师您统计的length_of_chromosome的长度是全部的碱基的长度,其中包含N,而我们统计的染色体长度是去掉N的(如下图所示),所以老师统计出来的覆盖率偏低(如21,22号染色体);
Y染色体之所以更低是因为Y染色体中含有很多重复的序列,而我们做比对分析的时候是把这些mask掉了(如果把这些也算进来会造成multimap的比例较高,不利于后面的分析),所以用来计算覆盖率的Y染色体的长度比实际的短很多。
老师可以参考我们给出的染色体长度(下图)重新计算一下就可以得出统计图中的结果了。
也就是说,他们在算chrY的覆盖度的时候,不是把Y染色体的全长58M来做分母,而是只用了22.98M来做分母,这样就算出了超高的覆盖度,接近于100%!!!
先不论公司这样做对不对,我首先重复一下他们的分析结果再说,不就是去除参考基因组的N碱基长度嘛,我还专门写了一个帖子:统计各种参考基因组的各条染色体的N含量()
的确,这样统计数据马上就漂亮很多,也就有了公司给我的测序深度和覆盖度的图。数据统计如下:
如下图,可以直接在excel表格里面做出来的。
好了,这个问题说清楚了,继续前进!
文:Jimmy、阿尔的太阳
图文编辑:吃瓜群众