分析带UMI标签的测序数据

检测癌组织的低频突变，为了提高检测低频突变的灵敏度，往往进行高深度的测序。但样本之间存在交叉污染，测序有存在一定概率的错误，这些因素会导致高深度测序过程中将假阳性的信号放到，得到假阳性的结果。解决交叉污染的方法，有公司比如IDT采用唯一配对的样本index，只有配对的index中的reads才属于特定样本。解决测序错误的方法，研究人员在建库的时候，先对分子加上UMI碱基，unique molecular identifier -> UMI，然后根据来源于同一个分子的测序数据进行测序错误修正，得到正确的分子序列。两种方法结合可以减少交叉污染提高准确性（https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5759201/）。

如图中所示，左侧一个分子被测了5次，其中第二次有一个测序错误，但该错误并没有在每个测序数据中出现，所以在后续合成一个分子的时候，测序错误被修正，只保留了真正的突变。（https://www.ncbi.nlm.nih.gov/pmc/articles/PMC5852328/）

常规的肿瘤配对测序分析，或者遗传性突变位点的分析，并不需要UMI信息，所以包含UMI的数据分析是需要不一样的分析流程来得到准确的分析结果，其中包括提取UMI分子标签，合并来自同一个分子的测序reads，低频突变检测而非胚系突变检测等。

大致流程为：

ubam是Unmapped BAM Format，是BAM文件的一个变种，里面的read是未经map的。大部分测序供应FASTQ文件，这是最常见的测序分析的起始文件。FASTQ文件的优势是，压缩比比bam文件好，解压速度快。

但与ubam文件相比，FASTQ并不是最理想的：

1）单个文件更容易分析

在双端测序中，有些软件希望配对的reads放在一个文件中，有些希望配对的文件，有些软件直接根据read在文件中的位置判断read是否配对，当然现在FASTQ通过在文件中加入/1和/2来表明配对的read解决这一问题。在生信分析的时候，一个FASTQ文件往往和另一个FASTQ文件关联配对，比如R1和R2，往往会花费更多时间来验证read是否配对。但如果通过bam文件的话，会更简单。只需要在FLAG这个地方加入对应的值，比如77和141，就能指定。单个文件更加简单，能储存更多metadata的信息。

在测序的时候，我们先拿到的是样本fastq压缩后的gzip文件，这个时候可能最关心的是数据量够不够，那么fastq.gz文件大小和测序数据量有什么关系呢。

我用Miseq测序数据（gz文件200M左右），Hiseq panel（gz文件50M左右）和WES测序数据（gz文件4G左右）进行了简单的分析。有意思的地方是，虽然R1和R2的数据量是一样的，解压出来的文件大小是一样的，但R2的gzip文件总比R1大。不管是Miseq还是Hiseq的panel测序，压缩后的R2均大于R1文件，且文件越小，差异越大。

一个全基因组或者全外显子组会产生几万、几十万乃至百万数量的突变。而目标突变往往只有几个或者几十个。我们需要大海捞针似得把我们想要分析的突变找到。突变过滤就是尽可能的将无关突变过滤掉，尽可能的保留目标突变，找出最可能致病的突变。关联性分析则不需要，因为关联性分析根据的是突变在对照和研究人群中的分布。

突变的过滤需要根据分析需要，比如需要寻找一个完全没有研究过的突变，则需要根据dbSNP来过滤等。如果需要寻找新生突变，要同时用到父母双方的突变数据。

常用的过滤条件可以分为一下几类： 1）根据突变频率进行过滤：每个位点的基因型在人群中频率是不一样的，一般来说，在人群中频率高的基因型往往没有致病性，而是背景突变，只是人群中的多态性。通常用到的频率数据库有gnomAD，ESP，ExAC，1000 Genome等，这些项目检测或搜集了从上千到成万的样本数据，并把数据开放给学术界和产业界。一般过滤的阈值有0.05，0.01，0.001等。 2）根据突变位置进行过滤：不同的分析目的，关注的突变并不一样。有人关于调控区域的突变，有人关注编码区的突变，有人关注剪切位点突变等。根据分析目的不同，需要将非关注的区域内的突变过滤掉。如果关注编码区突变，则需要将基因间、内含子、非编码区的突变过滤掉。 3）根据突变类型进行过滤：编码区的突变，有同义突变、错义突变、终止密码子突变等。同义突变虽然可能会导致疾病的发生，导致调控紊乱，但一般分析会将同意突变过滤掉（和分析目的有关，如果想要研究同义突变，则不能过滤）。往往移码突变和终止密码子突变是需要关注的。 4）根据突变危害预测进行过滤：现在有很多算法在预测位点导致的氨基酸变化的危害性，通常这些算法会分析致病的突变所处的区域、氨基酸变化类型等特征，比如PolyPhen、SIFT等，进而预测出检测出的突变危害性。有些算法会分析序列的保守型，预测突变所处区域是否保守，比如PhastCon，保守区域的突变较非保守区域的突变危害性大。此外，还可以根据氨基酸类型的带电性变化进行过滤，或者根据BLOSUM等打分矩阵分值进行过滤。不同的算法预测的危害性不一定一致。一般统计多个算法预测中，预测危害性较大的比例。 5）根据先验知识库进行过滤：注释软件可以注释关联性分析的信息，也可以注释Clinvar相关数据，这些信息在后续的突变位点解读和分析中有很重要的作用。这些信息可以告诉分析人员位点和哪些疾病有关，或者在Clinvar中的致病性评级。可以根据显著性、是否为良性突变、疾病信息等内容对突变 进行过滤。 6）根据数据库标识进行过滤：注释软件会注释突变在dbSNP、COSMIC等数据库中的编号。如果研究的疾病非常非常罕见，以往没有研究过，可以尝试寻找非dbSNP数据库中的突变。如果研究和癌症相关的突变，可以寻找COSMIC数据库中出现的突变。 7）根据遗传方式进行过滤：如果疾病呈现家族性分布，可以根据系谱图推断出的遗传类型进行过滤。比如如果是隐形突变，则只考虑病人中的纯合突变。如果有线索指示突变为新生突变，则需要根据父母双方的突变数据，将后代中同样位点的突变过滤掉。 8）根据其他条件进行过滤：突变过滤还可以有其他的方式，最终还是要根据分析目的进行选择。比如有四个散发样本的突变数据，同样的表型，往往寻找四个样本在相同基因内是否都发生突变，或者是否有相同的高危害突变（比如移码突变）。

很早就知道神器zsh的存在，一直没有用，今天手痒装了一下，哈哈，用起来还不错。

zsh是bash的增强版，可以兼容bash，可谓是shell中的战斗机，非常方便的命令提示、补全、忽略大小写等功能，一个字爽，可参见

https://www.zhihu.com/question/21418449

不过通常zsh配置起来比较复杂麻烦，于是诞生了oh-my-zsh，专门为简化zsh的配置而开发，http://ohmyz.sh，极大的简化了zsh的配置，还是一个字，爽。