A Byte of Biology

（更新时间：2023-11-26）

物种、分类群

我们的论文在动植物中对HapHiC进行了测试和验证，列表如下：

- 阅读完整内容 -

（更新时间：2023-11-24）

写在前面

近期，我们新开发的基因组挂载工具HapHiC在GitHub上公开了代码，相关论文的预印本也发布在了bioRxiv上，并有幸收到了大家热烈的反馈。然而，在这个过程中我发现虽然GitHub上已经有一份使用说明，许多朋友仍然对HapHiC的使用和原理不太了解。这可能是因为GitHub上的说明过于简单，而HapHiC的设计和参数又较为复杂的缘故。同时，越来越多的中国人加入到了基因组学研究的队伍中，并发表了大量高水平的研究成果。因此，我写下中文版的HapHiC系列教程与方法解读，供大家查阅和参考，希望能对HapHiC的用户有所帮助！

- 阅读完整内容 -

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

- 阅读完整内容 -

一个sratoolkit的报错

想到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

先做一点简单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）。下面几个图看一看应该很容易理解：

- 阅读完整内容 -