上一篇文章《肿瘤变异数据分析和可视化工具maftools：突变数据下载和可视化》主要以TCGA-LUAD为例介绍突变数据和临床数据的下载、处理以及简单的可视化。这篇文章更详细介绍可以利用maftools对肿瘤MAF格式的突变数据做哪些分析。
具体的数据下载和处理方法这里就不再赘述了，请移步上篇文章：《肿瘤变异数据分析和可视化工具maftools：突变数据下载和可视化》。这篇文章里，读入的临床数据终于可以派上用场了。
library(maftools)
luad <- read.maf(maf="TCGA.LUAD.maf", clinicalData="clinical.tsv")本文接上次的内容：《肿瘤变异数据分析和可视化工具maftools：安装和文件格式要求》，本文以TCGA肺腺癌（LUAD）的数据为例介绍突变数据下载和可视化。
在TCGA官网下载TCGA-LUAD项目中mutect2输出的MAF格式的突变以及临床信息：
$ ls -lhrt
total 50M
-rw-r--r-- 1 xiaofei xiaofei 50M May 14 22:44 TCGA.LUAD.maf.gz
-rw-r--r-- 1 xiaofei xiaofei 157K May 15 00:12 clinical.tsvMaftools是一款可以对MAF格式（Mutation Annotation Format）的变异数据进行统计、分析和可视化的R包。除了可以对TCGA来源的MAF文件以外，其他任何变异数据只要是MAF格式都可以使用这款工具进行分析。
Maftools包可主要概括为可视化和分析两大模块，流程和使用方法很简单：通过read.maf读入MAF文件（或者经过格式转换）得到MAF对象，然后将对象传递给对应的分析或者可视化函数就行了。主要模块、函数和主要的分析和可视化功能见下图：
VAF的全称是Variant Allele Frequency（变异等位基因频率）或Variant Allele Fraction（变异等位基因分数）。简单来说就是在基因组某个位点支持alternate/mutant allele的reads覆盖深度占这个位点总reads覆盖深度的比例。以VCF文件中的字段为例，其中DP代表Total Depth，AD代表Allele Depth，因此VAF的计算就是：
$$VAF = \frac{Allele\ Depth}{Total\ Depth} = \frac{AD}{DP}$$
VAF用得比较多的地方是在二倍体germline的genotyping中，杂合位点的VAF在高深度（比如depth>80）情况下应该接近50%；如果VAF接近0.25/0.75说明基因组上可能还有另一份拷贝。另一个应用场景就是癌症基因组的somatic genotyping。肿瘤组织、cfDNA、ctDNA、CTC genotyping的结果中会包含正常的allele（与正常体细胞一致）以及突变的allele，其中突变allele的所占的比例就是VAF。VAF可以用于推断肿瘤的异质性和肿瘤纯度，此外VAF的高低可能会影响癌症的预后。
线性回归问题即已知一系列样本的自变量和因变量的值，求解以下方程中的各θj：
$$h_θ(x) = θ_0 + θ_1x_1 + θ_2x_2 + \cdots + θ_nx_n$$
设样本数量为m，评估拟合的直线与实际样本之间差异的代价函数（Cost Function）为：
$$J(θ_0,θ_1,\cdots,θ_n) = \frac{1}{2m}\sum_{i=1}^{m}(h_θ(x^{(i)})-y^{(i)})^2$$
因此，寻找最佳拟合的线性回归模型则转化为求解该代价函数的最小值，常用方法为最小二乘法（Least Squares Method）和梯度下降法（Gradient Descent Method）。
最近听闻前公司和另一家公司A因为一个项目中二代测序数据PCR duplication的问题在扯皮。前公司比较两对Paired-end reads的第8~113bp区间内碱基是否完全相同（PE150）；而A公司把1~150bp都完全相同的才当作PCR duplication。最终参数的分歧使得两家公司的统计结果相差接近一倍，关乎到数据量是否的满足合同要求。
而前公司这部分QC流程的参数就是我设置的。选择截取reads中间的区域进行比较的原因很简单，因为Illumina测序前面几个碱基和末端的区域通常测序质量较差，测序错误较多。如果要求两对reads完全一致，相当于600bp（150*2*2）不能出现任何测序错误。具体为什么设置成8~113bp我并没有详细去分析和论证，只是因为当时测序都外包给N公司做，而N公司内部QC流程去PCR duplication设置的范围设置的就是这个，我也就顺理成章用了同样的参数。
那么怎样去除二代测序数据中的PCR duplication才是科学的姿势呢？虽然对很多人来说应该是一个无聊的论题，但是我还挺好奇，本文姑且在这里分析一二。
在基因组de novo组装后，为了进一步提高组装连续性经常对初步组装结果进行scaffolding。这时候可以直接利用手头已有的二代PE/MP/BAC文库测序（SSPACE）、三代单分子测序（SSPACE-LongRead、LRScaf）甚至转录组数据（L_RNA_scaffolder）。但是仅靠这些通常很难在保证准确性的前提下大幅度提升基因组连续性，因此通常不得不投入更多经费加测10X genomics、BioNano光学图谱或者Hi-C数据。
如果研究的物种已有较高质量的参考基因组时，比如水稻、拟南芥、人等常见模式生物，且又不关心测序个体可能存在的结构变异时，可以直接利用参考基因组scaffolding到染色体水平（依赖参考基因组的连续性）。这里整理了几个可以利用参考基因组进行scaffolding的一些工具。
print对很多人来说算是最常用的debug神器了，只需要在合适的地方插入print打印变量的值，就能判断代码在这个地方是不是还按照你的预期在运行。不过这也带来一些麻烦：首先需要找准位置，然后写出对应的print语句。当函数复杂、变量多的时候确实挺烦的，通常需要多个地方插入多个变量的print。
最近一个刚刚上线的Python包“PySnooper”更优雅地解决了这一需求——只需要对关注的函数前加上一个装饰器@pysnooper.snoop()，就可以将这个函数运行的时间，行号、运行过程中变量的数值以及代码等内容输出到stderr。项目刚在GitHub上线1天左右，已经有超过1300个star，可以说是相当火爆了（更新：该项目成为2019.04.23 GitHub每日趋势榜第一名）：