技术分享
包括但不限于生物信息学分析软件、流程、分析方法、编程、算法、集群管理和运维等技术分享

肿瘤变异数据分析和可视化工具maftools:安装和文件格式要求

Maftools系列文章:

  1. maftools使用方法总结以及常见问题
  2. 肿瘤变异数据分析和可视化工具maftools:安装和文件格式要求
  3. 肿瘤变异数据分析和可视化工具maftools:突变数据下载和可视化
  4. 肿瘤变异数据分析和可视化工具maftools:突变的数据分析
  5. 肿瘤变异数据分析和可视化工具maftools:CNV的可视化

Maftools简介

Maftools是一款可以对MAF格式(Mutation Annotation Format)的变异数据进行统计、分析和可视化的R包。除了可以对TCGA来源的MAF文件以外,其他任何变异数据只要是MAF格式都可以使用这款工具进行分析。

Maftools包可主要概括为可视化和分析两大模块,流程和使用方法很简单:通过read.maf读入MAF文件(或者经过格式转换)得到MAF对象,然后将对象传递给对应的分析或者可视化函数就行了。主要模块、函数和主要的分析和可视化功能见下图:

- 阅读完整内容 -

VAF,MAF,肿瘤纯度,MCF,CCF的概念和计算方法

VAF的概念和计算方法

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)。

- 阅读完整内容 -

如何去除二代测序数据中的PCR duplication才科学?

最近听闻前公司和另一家公司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才是科学的姿势呢?虽然对很多人来说应该是一个无聊的论题,但是我还挺好奇,本文姑且在这里分析一二。

- 阅读完整内容 -

利用参考基因组进行scaffolding提高组装连续性的工具整理

在基因组de novo组装后,为了进一步提高组装连续性经常对初步组装结果进行scaffolding。这时候可以直接利用手头已有的二代PE/MP/BAC文库测序(SSPACE)、三代单分子测序(SSPACE-LongRead、LRScaf)甚至转录组数据(L_RNA_scaffolder)。但是仅靠这些通常很难在保证准确性的前提下大幅度提升基因组连续性,因此通常不得不投入更多经费加测10X genomics、BioNano光学图谱或者Hi-C数据。

如果研究的物种已有较高质量的参考基因组时,比如水稻、拟南芥、人等常见模式生物,且又不关心测序个体可能存在的结构变异时,可以直接利用参考基因组scaffolding到染色体水平(依赖参考基因组的连续性)。这里整理了几个可以利用参考基因组进行scaffolding的一些工具。

- 阅读完整内容 -

还用print在Python中debug吗?试试PySnooper吧!

print对很多人来说算是最常用的debug神器了,只需要在合适的地方插入print打印变量的值,就能判断代码在这个地方是不是还按照你的预期在运行。不过这也带来一些麻烦:首先需要找准位置,然后写出对应的print语句。当函数复杂、变量多的时候确实挺烦的,通常需要多个地方插入多个变量的print

最近一个刚刚上线的Python包“PySnooper”更优雅地解决了这一需求——只需要对关注的函数前加上一个装饰器@pysnooper.snoop(),就可以将这个函数运行的时间,行号、运行过程中变量的数值以及代码等内容输出到stderr。项目刚在GitHub上线1天左右,已经有超过1300个star,可以说是相当火爆了(更新:该项目成为2019.04.23 GitHub每日趋势榜第一名):

https://github.com/cool-RR/PySnooper

- 阅读完整内容 -

Python中列表的 += 和 .extend() 的异同

一道Python题

最近有朋友“考”了我一个Python的题:使用+=.extend()两种方法扩展元组中的列表会发生什么。虽然我对Python中的可变数据类型、不可变数据类型的概念都有较深的理解,并且也对list的++=.extend().append()做过性能分析,但是+=.extend()两者无论在表现(是否为原地址修改)以及性能上都非常近似,所以对两者的区别还没有明确的概念。为了解答这个问题,我我们先直接上代码试验一下:

# 创建一个包含列表的元组:
>>> a_tuple = (1, 2, [])
>>> a_tuple[2] += ['a', 'b']        # (1)
Traceback (most recent call last):
  File "<pyshell#14>", line 1, in <module>
    a_tuple[2] += ['a', 'b']
TypeError: 'tuple' object does not support item assignment
>>> a_tuple[2].extend(['a', 'b'])   # (2)
>>> a_tuple                         # (3)
(1, 2, ['a', 'b', 'a', 'b'])
  • (1) 通过+=的方法扩展列表出现“元组不支持元素赋值”的报错。
  • (2) 使用.extend()方法。
  • (3) 有趣的是,列表被扩展了两次。虽然+=报错,但是却成功修改了列表。

- 阅读完整内容 -

awk命令坑了多少科研工作者?

著名“错误使用”案例——Excel的单元格格式

之前看过一篇文章,讲述的是2016年发表在《Genome Biology》上的一项研究报道——Microsoft Excel中的单元格格式的自动转换导致研究的期刊中将近20%的论文存在基因名转换错误。比如会把RIKEN identifier “2310009E13”作为科学计数法转换成“2.31E+19”,或者将基因名“SEPT2”作为日期转换成“2018/09/02”。其实早在2004年该问题就已经在《BMC Bioinformatics》上有相关报道,但似乎没引起警觉。不过Excel作为一款在各种领域广泛使用的电子表格软件,微软不太可能为了这类问题而牺牲其“智能”、“易用”的特色。并且这也只能算是用户的“错误使用”,而不是Excel的bug。下面是两篇论文的原文,有兴趣可以一看:

不得不小心的awk之坑,Caveat Emptor!

以上是联想到的题外话,本文主要介绍的是“GNU awk”(即gawk)中的一个坑。和Excel的故事非常类似,都可以归咎为理解不够深入而导致的“错误使用”。而且相比Excel格式转换这种错误,这个问题隐蔽性更高,即使对编程语言的精度问题有所了解,但是还是会掉到坑里,甚至让我一度以为awk出现了bug。

- 阅读完整内容 -