2020年各行各业受到疫情的影响应该都比较大，搞科研的也不例外。年初做湿实验的长达一个多月无法正常返校工作，还面临各种进口试剂和实验用品减产和订购周期延长的问题。相比之下我们做生物信息学的算比较幸运的，可以在家远程工作。
在家工作，效率上或多或少会下降。本来计划完成两个基金的撰写，最终权衡精力和概率，放弃了国自然青年，选择了专攻博士后面上。由于我是第一次写基金缺乏经验，理解其中“研究对象”，“拟解决的关键科学问题”，“研究目标”等概念着实让人头大。好在反复修改了几遍本子之后，对这些内容的理解深刻了不少。另外为了这个基金，在前期还花了不少时间进行预实验，利用公开数据做了一些数据分析并进行了可视化。一方面证明研究可行，另一方面体现申请人具备完成课题的能力。因为这个基金是以博士研究方向为基础，最终只让实验室同僚们稍微修改了语言就提交了。
在Python中需要初始化字典的时候我一般会调用collections中的defaultdict来实现，这样可以减少代码行数，并且避免更复杂的逻辑判断，比如：
a_list = ['a', 'b', 'c', 'a', 'd']
# 常规方法
count_dict = {}
for a in a_list:
if a in count_dict:
count_dict[a] += 1
else:
# 如果dict中不存在该key，则需要赋值来初始化
count_dict[a] = 1
# 使用defaultdict进行初始化只需要4行代码，且不需要if...else进行逻辑判断
import collections
count_dict = collections.defaultdict(int)
for a in a_list:
count_dict[a] += 1想到Linux中~和$HOME的区别这个问题，是因为在集群上安装了sratoolkit，在使用fastq-dump时遇到的一个报错：
$ fastq-dump --split-3 A549-DRR016695.1
2019-11-30T16:49:38 fastq-dump.2.9.6 err: path not found while searching dynamic library within file system module - failed to open SRA manager
=============================================================
An error occurred during processing.
A report was generated into the file '/work/xxxxx/ncbi_error_report.xml'.
If the problem persists, you may consider sending the file
to 'sra-tools@ncbi.nlm.nih.gov' for assistance.
=============================================================同样fasterq-dump也无法正常使用：
$ fasterq-dump --split-3 A549-DRR016695.1
2019-11-30T16:50:19 fasterq-dump.2.9.6 err: cmn_iter.c cmn_get_acc_type( 'A549-DRR016695.1' ).VDBManagerMakeRead() -> RC(rcFS,rcDylib,rcSearching,rcPath,rcNotFound)增强子RNA（Enhancer RNA，eRNA）近期好像因为Nature Communications的一篇文章又火了一下，作者还甩出了一个叫eRic的数据库[1]，里面不仅整理好了TCGA各类癌症中鉴定到的eRNA，表达量，还关联上了病人的生存，甚至还有药物处理后的变化。临床和实验基础的课题组看到这篇文章说不准欢呼雀跃快马加鞭做起了功能机制研究。不过在做新方向课题之前，有没有好好想过，别人的分析真的靠谱吗？
先做一点简单eRNA背景介绍：（1）激活的增强子被发现通常会转录产生RNA，称作eRNA；（2）eRNA有个特征就是会双向转录（bidirectional transcription）；（3）根据双向转录这个特征，FANTOM5通过对808个人类组织或细胞（目前已经更新到了1829种）进行CAGE-seq（Cap Analysis of Gene Expression）鉴定了目前最“靠谱”的eRNA/激活的增强子数据集，并于2014年发表了一篇Nature [2]。为什么这里“靠谱”要用引号，本文的最后一节会有解释。
众所周知，描述基因组中位置信息的坐标系通常有1-based和0-based之分。举个例子，要表示某条序列的前5个碱基，如下图的ATGCA，用1-based坐标系就应该是1-5，而0-based坐标系应该是0-5（类似与Python中的切片）:
常见的1-based坐标系的有SAM、GFF/GTF、VCF、WIG、Axt、BLAST比对结果；0-based坐标系的文件有BAM、PSL、MAF（这里指的Multiple Alignment Format）、以及BED等。如果坐标区间范围很大，只做一些简单分析或者可视化，这一个碱基的增减可能影响不大。如果区间很小，比如TSS（Transcription Start Site，转录起始位点）；又或者是要根据编码基因coding region的DNA序列转换成对应蛋白的氨基酸序列，那这一个碱基的差别可能就会带来严重错误。
因为最近博士后入职一堆事情，有一段时间没有更新文章了，内心充满罪恶感，这篇文章也在我草稿中躺了很久了。想起这个主题的起因是前几个月给大名鼎鼎Nature子刊审稿（其实是Scientific Reports），发现文章里做各类相关性分析的时候Pearson和Spearman相关系数随缘混用，基本瞎搞。虽然在这个水平杂志上出现这种问题似乎不会让人感到意外，但是意外的是文章内容做得相当扎实，一作也不是新手，而且还在世界顶尖大学里生物信息学方向的课题组。再到后来入职博士后，发现在我理解里统计学水平应该比生物更强的医学专业学生或者老师，对这个也不是很有概念的样子，所以觉得很有必要拿出来讲讲。
首先分别简单介绍一下我们可能最常用到的Pearson和Spearman相关系数是什么东西，不写公式，尽量说得通俗易懂（其实是怕麻烦）。Pearson相关系数很简单，是用来衡量两个数据集的线性相关程度。而Spearman相关系数不关心两个数据集是否线性相关，而是单调相关，Spearman相关系数也称为等级相关或者秩相关（即rank）。下面几个图看一看应该很容易理解：
本文作为总结篇，将在这里对已经写过的几篇文章内容进行概括，可作为整个系列文章的目录，此外还有一些常见问题的解决方法。加上本篇目前一共写了5篇maftools相关文章，基本上按照整个官方文档使用TCGA-LUAD的数据都跑了一遍，并且加上了一些数据整理、重要参数的解释以及自己遇到的问题，应该可以说是目前最完整的中文资料了。其实maftools本身使用起来很简单，在读入数据之后，基本可视化和数据分析通过1~2行代码就能实现，参考官方文档的example足够了。我写得比较细、花的篇幅较多是因为自己也是初学肿瘤的数据分析，通过学习maftools这样功能丰富的分析工具可以更快速入门。
上篇文章《肿瘤变异数据分析和可视化工具maftools：突变的数据分析》介绍了使用maftools分析MAF格式的突变数据，但maftools本身功能不限于此，它也可以针对拷贝数变异进行一些分析。
还是继续之前TCGA-LUAD的例子。Maftools接受的CNV数据需要是GISTIC输出的结果，包括4个文件——all_lesions.conf_XX.txt、 amp_genes.conf_XX.txt、del_genes.conf_XX.txt以及scores.gistic（其中XX代表置信水平）。但是从TCGA官网能下载到的GISTIC的输出结果只有LUAD.focal_score_by_genes.txt，因此需要下载更上游的DNAcopy的输出结果，然后自己再根据TCGA的流程和参数跑一遍GISTIC。