在学习bedtools的过程中偶遇一个非常棒的网站，就是其作者 Aaron R. Quinlan的实验室网站http://quinlanlab.org/
bedtools这个软件工具支持对BAM、BED、GFF/GTF、VCF等和基因组区域有关文件的intersect, merge, count, complement, and shuffle等操作。这个解释有点抽象，不过总的来说就是一个非常棒的工具，曾老师曾经强烈推荐初学者跟着学习的。
2010年这个工具被发在了Bioinformatics上，作者当时还在美国弗吉尼亚大学医学院，现在工作在犹他大学，他的博客有篇比较新的博文解释他们实验室最新发现http://quinlanlab.org/blog/2018/12/20/constrained-coding-regions.html
对基因组突变位点的解读，我觉得非常有意思，就记录了下来。其中有个小工具也非常棒，对bed文件等（BAM, VCF, GTF, or BED files(indexed BAMs, other files bgzipped and tabixed)）用网页工具进行IGV可视化的。https://s3.us-east-2.amazonaws.com/ccrs/ccr.html
这篇博客的一开始引用了2个历史人物的语录，因为翻译的水平有限，就直接摘录原文吧
Once in a while you get shown the light
In the strangest of places if you look at it right.
- Jerry Garcia and Robert Hunter, Scarlet Begonias
The unseen enemy is always the most fearsome.
- George R.R. Martin, A Clash of Kings
为什么作者一开始会引用这2个语录呢？往后看了才明白他的意思，在博文中他先引入了一张图后面又讲到一个来自维基百科的故事。
在这篇博文中，作者提供了他们发表在eLife的手稿的背景和概述，题目是"A map of constrained coding regions in the human genome"。简单地说，作者研究了从基因组聚合数据库(GnomAD)的2.0.1版中检测到的超过120,000个人类外显子的遗传变异，以揭示由于非典型变异的缺乏而受到限制的焦点编码区(例如，在上图所示的区域)。这些“受限编码区”(CCR，constrained coding regions)被推断为处于强烈的纯化选择（purifying selection）之下，并富集了已知的致病变异。也许这张CCR图谱最耐人寻味的方面是这样一个事实，即许多最受限制的区域位于缺乏先前疾病关联的基因内。这些区域有望在发育障碍（developmental disorders）的背景下发现新的疾病基因。
首先是解释人类基因组中最受限制的区域是什么？
人类遗传学的一个长期兴趣是确定我们基因组中对生命和正常发育最重要的子集。一般来说，这样的区域应该处于最高的纯化选择之下，因此应该表现出较低的核苷酸多样性。在蛋白质编码基因的情况下，应特别注意对蛋白质改变(i.e., missense, stop-gain, frameshift, etc.)的strong “constraint”的variants。事实上，这个概念强调了最近的“gene-wide” constraint指标背后的动机，比如Residual Variation Intolerance Score (RVIS)和较新的probability of Loss-of-function Intolerance (pLI) score。虽然这些指标已被证明对罕见疾病的研究非常有用，但单一的、全基因的指标本质上不能描述每个蛋白质编码基因中存在的constraint的区域差异。确定constraint的焦点区域是作者们的动机。下面是作者结合部分已发表文章的历史、主要成果和数据文件，以及他们对未来重要研究的看法。
背景部分
作者讲到2010年末的某个时候，当他正紧张地准备成立自己的研究小组时(他说自己甚至无法表达他有多害怕)，他看到了一场演讲，描述了NHLBI外显子测序项目对大约6000个人类外显子进行测序的目标。作者(可能还有其他许多人)很快意识到，他们可以利用在这些外显子中发现的变异，从这些样本中没有变异的情况下推断出处于纯化选择之下的编码区。与才华横溢的比尔·皮尔森进行了多次有趣的讨论，将研究集中在一起。不久之后，吉姆·哈夫里拉(Jim Havrilla)加入了他的实验室，很快就在他的博士课题中对这个想法进行深入的研究。
在接下来的8年中这个想法的核心催生了该课题组在人类基因组中受限编码区图谱的创建。这项工作从头到尾都由吉姆·哈夫里拉(Jim Havrilla)领导，同时得益于布伦特·彼得森(Brent Pedersen)和瑞安·莱尔(Ryan Layer)的聪明想法、巧妙的代码应用。
题外话：
从背景部分我们看到了一个年轻的研究员初创自己的团队在助理和博士研究生的共同协作下奋战了8年的一个结果。科研之路非常艰辛，但是团队的协作总能完成最后的工作，而且几位年轻人可以说是白手起家用公共数据完成学术研究，并有了独到的见解，非常值得我辈学习，我们国内的博士培养很多都是类似的模式，可是真正和导师一样为了同目标非常有兴趣的进行研究工作的还是少。非常羡慕这种纯粹的有意义的研究工作。这可能与很多制度问题相关，这里就不展开吐槽了。
作者研究的这个想法是基于survival bias的概念，这个概念在科学中很普遍，在二战期间亚伯拉罕·沃尔德和统计研究小组(SRG)的工作中得到了最著名的证明。这是一个来自维基百科的故事：盟军飞机被击落，军方领导层显然想阻止飞机损耗速度。然而，金属是稀缺的。此外，在添加金属进一步保护飞机的同时，也会降低飞机的机动性和燃油效率。这是一个经典的优化问题--他们如何在最大化防御的同时使用最少的金属？SRG收到了描述从他们出动归来的数百架飞机上观察到的弹孔图案的数据。据说军方领导层对这一数据的解读是，装甲应该放在弹孔最密集的地方(这就是我们被枪击的地方！)。沃尔德不同意这一点。他认为，观测到的数据是有偏差的，因为它们只来自返回(幸存)的飞机。他认为，装甲应该放在没有子弹的地方，因为这些地区很可能是被击落的飞机受损的地方。换句话说，这些都是受限制的飞行区域。
就是从上面这个小故事激发了作者的idea。
作者团队使用survival bias来识别人类基因组中的限制性(即在强烈的纯化选择下)编码区(CCR)，因为在大于12万个人类外显子中没有蛋白质变化的变异很多。
关键性结果展示
正如手稿中详细描述的那样，作者将限制性编码区确定为蛋白质编码基因片段，在基因组聚合数据库(GnomAD)中的120,000个外显子中甚至没有一个蛋白质variant 的基因。虽然gnomAD中这种变异出现的的平均密度约为每7个coding bases中包含1个，但最受限制的编码区(例如，在第99个百分位数或以上)通常缺少超过100个碱基的蛋白质改变variant 。例如，下面的红色区域反映了我们在KCNQ2和TNNT2中确定的第95个百分位数和更高的CCR。
作为一个阳性对照证明，最受限制的编码区富集了ClinVar中已知的罕见人类疾病表型背后的致病变异。例如，最受限制的区域之一是在SCN8A中没有蛋白质改变变异的274编码碱基对区域。组成该CCR的4个外显子编码大部分离子运输结构域。下面是使用IGV.js构建的CCR浏览器的屏幕截图。图像有点难以辨认，因此可以通过此链接直接查看该区域。暗红色区域反映第99个百分位数或以上的CCR。
作者发现CCR补充了其他不同的优先排序工具，用于解释罕见疾病背景下的
de novo突变。他们认为，位于最受限(例如，第99个百分位数或更高)编码区的de novo突变可能与发育表型有关。事实上，虽然它没有被写入手稿，但在他们团队最近对早期婴儿癫痫性脑病的研究中发现的几乎所有致病突变都位于CCR内，位于第95个百分位数或以上(大多数在第99个百分位数以上)。
因此，作者认为，在发育障碍的背景下，位于regions of the highest constraint内的de novo突变特别令人感兴趣。然而，一个重要的警告是，不能简单地忽略less constrained regions的突变，因为许多已知的致病等位基因位于表现密集变异的区域(例如，BRCA1)。
另一个有趣的结果是，观察到物种内的constraint往往与物种间的conservation measurements相辅相成；也就是说，conservation measurements并不总是预测物种内的constraint。此外，我们还确定了具有最大限制的蛋白质结构域家族的子集。从高水平上看，这些结构域通常与DNA相互作用或修饰染色质。关于这些领域中的constraint的大部分细节在Boukas等人的预印本文章中有更详细讲述。
是否会提示新的致病基因呢？
考虑到已知致病变异的高CCRs的富集，最令作者和我们兴奋的结果是许多高度(>99%)受限制的区域位于缺乏先前疾病关联的基因内。当然，其中一些是假阳性。然而，可以预想假设其中一些区域反映了强烈的净化选择，当它们发生突变时，会导致发育表型，甚至胚胎死亡。
作者说他们很高兴在未来的研究中探索这些区域，并希望受限编码区的地图的绘制将有助于指导我们未来的研究，并使罕见疾病研究中的突变解释成为可能。此外已经有不少报道已经证实了一想法(参见Jensen et al, Wray et al, and Boukas et al)！也是我们可以研究的热点。
作者强调他们的识别受限编码区(CCRs)的方法的优雅之处在于它非常简单。然而，它故意非常严格，因为他们希望将误报降至最低。诚然，基于gnomAD中存在的单个改变蛋白质的variant而“破坏”constrained regions可能会导致假阴性；也就是说，被单个variant破坏的更大的constrained regions。我们强调，我们创建的map揭示了dominant model下的constrained regions，而不是很适合于recessive constraint。最后，尽管gnomAD很强大，但它主要由欧洲血统的个体variation而成。因此，我们的map模型在多大程度上恰当地模拟了其他人种中的constrained regions，目前还不清楚。
在这些年里，大量的人类基因组将被测序。此外，由于gnomAD和其他努力，人们对人类基因组学中的数据共享做出了令人兴奋的commitment。因此，作者希望这项研究，以及Kaitlin Samocha的类似想法，能够启动新的方法来分离我们基因组的关键区域。预计，来自更多人类基因组的变异将提高预测处于强烈纯化选择下的区域的分辨率和准确性。同样，多亏了gnomAD、Genomics England、TopMED和Centers for Common Disease Genomics，他们正在迅速分析数据集，包括超过10万个基因组的全基因组变异。这些数据集蕴含着一个有趣的前景，即在我们基因组的非编码区对人类的constraint进行建模。
最后作者强调：我们将继续使用未来版本的gnomAD以及其他资源来更新我们的CCR map。我们也期待着将出现的新方法，并渴望继续我们在这一领域的研究。敬请关注。
