最近在看circRNA的文章，检测工具从2013年到现在已经有了几款，但是每个软件的检测结果都差异很大。

在biostar上，有人在最近几天讨论了2个国人做的circRNA检测工具，其中的问题还是很大的。

现在就我使用过的circRNA detection tools 做一个简要的介绍。

1. find_circ

这款工具是在2013年随着Circular RNAs are a large class of animal RNAs with regulatory potency这篇文章发布的，从而吹起了一阵研究circRNA的风潮～。

现有的所有circRNA工具的第一步都是mapping，将RNA reads mapping到基因组上，但是普通线性RNA是可以很好的mapping回去的，对于circRNA来说，成环的剪接位点不能直接mapping回基因组上。 首先，就是用比对mapping的方法筛选出这样的序列，也就是把mapping不上的序列都体取出来。 然后，将这些mapping不上的reads取两头20bp，变成短序列，重新mapping到基因组上。

这里有个参数20bp的reads，为什么要用20bp的？是否可以用更长的，或者更短的？

我们来做个简单计算。多长的一条序列可以唯一比对到基因组的某个地方？ 首先，基因组的大小为3G，也就是三十亿字节。\(4^{15}=1073741824\)，约为十亿字节。 那么任意取出一条15bp的read，它的碱基排列情况是4^15的其中一种，任意取出两条15bp的reads 它们完全相同的概率为\(\frac{1}{4^{15}}*\frac{1}{4^{15}}\)。 从基因组上的任意取出15bp长的一段序列，有三十亿种可能性，取两次， 取出同一条15bp长的序列的可能性为\(\frac{1}{3*10^{9}}*\frac{1}{3*10^{9}}\)， 这种可能性(概率)小于任意从15bp序列的全部碱基排列情况中抽出2次相同碱基序列的概率，所以说15bp的read可以不唯一比对到基因组的某个地方。 同理计算16bp reads \(4^{16}=4294967296\)，所以任意抽出一条16bp的序列可以唯一比对到基因组上。（可能算得不对，如果有人发现错误，请立即指正，谢谢）

这样看来20bp的序列，可以唯一比对到基因组上，也可以根据测序reads的长度适当放大缩小（我觉得缩小还是没必要的，可以放大，取25bp，30bp之类的）。

接下来，将短序列mapping到基因组上后，他们开发了一个方法，检测这些短序列是否是circRNA的短序列。

需要检测的条件如下：

1. GU/AG 在剪接位点的两侧出现
2. 可以检测到清晰的breakpoint
3. 只支持2个mismatch
4. breakpoint不能再短序列（anchor）以里2 nucleotides 之外的地方出现，也就是最多2nt （这条可能理解有偏差）
5. 至少有两条reads支持这个junction
6. mapping正确的一个短序列的位置要比它mapping到其他位置的分值高35分以上。

通过这些条件，就筛选出了文章认为正确的circRNA。

另外在这篇文章里重点介绍了circRNA的一个功能，富集microRNA，某些circRNA上有很多microRNA的seed。

2. CIRCexplorer

这款工具是在2014年随着Complementary Sequence-Mediated Exon Circularization的发布而问世，这篇文章我非常欣赏，我认为它是去年做RNA领域最值得读的一篇文章。

CIRCexplorer这个工具巧妙的用fusion gene这个思路去检测circRNA。 首先，过滤出Tophat无法mapping上的reads，再把这些reads用Tophat-Fusion mapping到基因组上。那些用Tophat-Fusion mapping到基因组上的非线性候选位置的reads就是潜在的back-splied juction reads。 接下来，这些reads会在有基因注释的帮助下，确定更加精确的donor, acceptor位置。 最后，对circular RNA进行注释。

这篇文章里介绍了在intron区域的反向重复Alu序列是引起circRNA形成的一个原因（即“内含子配对驱动环化”（intron-pair-driven circularization）模型，除此之外还有一种模型是“套索驱动环化”（lariat-driven circularization）模型），这些序列在跨越外显子的内含子区间上可以形成互补序列，所以在转录过程中容易形成茎环结构。环的部分就是外显子，然后在剪接酶的参与下，被剪接成了环形。

另外，circRNA也存在很多的可变环形结构，这也是由于在基因组中广泛存在的内含子上的Alu序列造成的，不同的Alu序列互补，形成含有不同exon的环形结构。

强烈建议研究circRNA的工作者好好学习一下这篇文章。附薛宇老师对这篇文章的评价，薛宇老师的博客是我看科学网唯一的动力。

3. CIRI

这款工具是中科院北京生命科学研究院的赵老师组的工作，今年发表在genome biology上， CIRI: an efficient and unbiased algorithm for de novo circular RNA identification，是一篇方法工具类的文章。

文章中总结了目前从RNA-seq数据中识别circRNA所遇到的问题：

1. circRNA比其他RNA在细胞中的比例低，一般RNA-seq的实验步骤中，不包括circRNA富集的步骤（例如RNase R treatment，我觉得他们要表到的意思是在RNA-seq中能探测到的circRNA比用RNase R处理过的数据中少），所以在RNA-seq中的circRNA比较低，假阳性也高。
2. 目前的注释文件都是根据线性RNA进行的，所以也不适合circRNA的识别，另外非模式动物的注释信息更少，就更不用谈在那些物种里找circRNA了。
3. 由于RNA测序数据的reads长度差别大，也对检测circRNA工作带来了不便
4. 套索结构和融合基因同circRNA的reads结果类似，不好区分。

文章中总结了从2012年到2014年出现的几种检测circRNA的方法：

• Salzman Cell-type specific features of circular RNA expression 用了一个依赖注释信息的方法来检测circRNA，通过搜索已知注释的外显子边界来查找，并在最近的工作中更新了方法，加入了false discovery rate 控制比对的质量分数。这种方法存在上述描述的第二个问题。
• Memczak Circular RNAs are a large class of animal RNAs with regulatory potency （也就是第一部分介绍的那个软件）用了GT-AG信号来找寻splicing 位点，也有其他工作用类似的方法筛选micorRNA-sponge 的候选circRNA。这种方法会找不到“长外显子1-短外显子-长外显子2”形成的环形结构，这种结构中一条测序Read上会有三个部分，第一部分序列属于长外显子1，第二部分序列属于短外显子，第三部分序列属于长外显子2。 Memczak的方法只是把一条序列切成两部分，这种算法会把“长外显子1-短外显子-长外显子2”丢掉，或者识别成“长外显子1-长外显子2”。
• Jeck Circular RNAs are abundant, conserved, and associated with ALU repeats 采用了比较的方法，比较没有经过RNase处理和经过RNase处理的序列的结果，用来确定潜在的circRNA，消除假阳性。这种方法在富集circRNA阶段会有系统误差。

说了那么多其他人的工作，赵老师组的工作是采用sam格式中的CIGAR值进行分析的，从sam文件中扫描PCC信号（paired chiastic clipping signals）。 CIGAR值在junction read的特征是xS/HyM或者xMyS/H，其中x,y代表碱基数目，M是mapping上的，S是soft clipping，H是hard clipping。

对于单外显子成环，或者“长外显子1-短外显子-长外显子2”形成的环形结构，CIGAR值应该是xS/HyMzS/H以及(x+y)S/HzM或者xM(y+z)S/H，CIRI软件可以很好的将这两种情况分开。

对于paired-end reads，CIRI算法会考虑一对reads，其中一条可以mapping到circRNA上，另一条也需要mapping到circRNA的区间内。

对于splicing 信号(GT,AG) CIRI也会考虑其他弱splicing 信息例如（AT-AC），算法会从GTF/GFF文件中抽取外显子边界位置，并用已知的边界来过滤假阳性。

CIRI方法消耗的内存比较大，我跑了个12G的sam文件，内存用了20G。

从mapping工具中找backsplice junction

除了上述直接找circRNA的工具外，还可以用mapping工具直接找backsplice的junction，然后在人工判断这些junction，例如segemehl。

#Create index:
./segemehl.x -d hg19.fa -x hg18.idx
#Map you reads:
./segemehl.x -d hg19.fa -i chr1.idx -q reads.fastq -S | samtools view -bS - | samtools sort -o - deleteme | samtools view -h - > mapped.sam
#Call different splice junctions:
./testrealign.x -d hg19.fa -q mapped.sam -n 

最后一步可以生成两个文件，一个是splicesites.bed，另一个是transrealigned.bed，从splicesites.bed文件中可以得到backsplice junction，在文件中的指示标志是”C”。

从我的测试样本中可以找到4085个通过检验的backsplice junction。

方法比较

从一开始的biostar链接中就可以看到不同方法的差异非常大，我测试了一个实际数据，从第一个软件有126163行结果，筛选后剩下491行，第二个软件只有327行，第三个软件有687行。

我认为不同的mapping策略从一开始就会导致搜索circRNA的空间不一样。

下游软件

画图

https://github.com/dieterich-lab/CircTest

中文介绍

环形RNA分子揭秘——张 杨 1 ,王海滨 1,2 ,陈玲玲 1 *

（未完待续）

标题是哗众取宠，Eric Lander在新英格兰医学上发表了一篇perspective(前景，观点)， 总结了基因组研究领域的新内容（CRISPR，基因编辑）和挑战，目前我们应该关注的问题，用谨慎而发展的眼光看待基因组编辑。

原文标题是Brave New Genome （美好的新基因组），我估计这个标题是类比于brave new world（美丽新世界）。

目前我们遇到的新挑战是开发一个可行的度量体系来检测人类生殖细胞编辑（human germline editing，说生殖细胞也不太顺，反正就是可遗传的细胞DNA的编辑）。

在此Lander认为应该关注四个关键问题讨论。

1. DNA编辑技术

目前的编辑技术还没成熟到可以自由的编辑指定基因，总会有不完全编辑、不精确编辑和脱靶的现象出现。

2. 这种编辑是不是利大于弊

目前人们想到的对于可遗传细胞的编辑可产生的一个主要贡献是毁灭单基因遗传病，例如亨廷顿病。 对于这些遗传病最有必要的不是基因组剪接，而是常规的遗传学筛查，使得携带这些遗传病的夫妇可以了解自己基因中潜在的风险并且愿意在生育前做产前遗传病诊断。

另外，编辑基因可能会使后代不得某种病，单也没准会加大他们患另外一些疾病的概率。这是应为，基因在生物体内主导的作用不仅仅同某一种疾病相关，更多是同很多疾病相关。 一个基因可能正向调控疾病A，反向调控疾病B。

3. 谁能决定是否进行基因编辑

这涉及到伦理问题，假设父母为了后代的健康、智力或体力因素，做了基因编辑，但是后代出生后却不同意父母的做法，这样该怎么办？

4. 对与错，今后该如何

科学家一般都不情愿去回答一些伦理问题，但是我们必须对这些问题负责并正确的对待这些问题。我们如何生活在一个存在基因编辑的社会，基因编辑会不会改变我们的世界？ 如何看待进行过基因编辑的的后代（他们是不是人工产品）？会不会出现基因编辑的流行趋势（大家都喜欢做某项基因编辑 ，变得更漂亮更聪明之类的）？携带最优秀的基因组的人会不会有更大权力？

如果我们跨越了这些界线，很难看到挽回的方法。

现在各国对于基因编辑都还是非常保守的态度。

但Lander在最后又说，我们很容易就可以推翻这种保守观点，前提是，我们都有很好的科学知识，在道德上足够明智，并且我们的研究带来了引人入胜、不可抗拒的结果。

十年前人们才第一次读到人类的基因组，现在我们需要小心而谨慎的前进（改写基因组）。

读完文章发现同我想象的内容完全不同，Eric Lander太会说话了，严谨而不失希望。

He Jiankui老师组在去年年初总结了一下目前单细胞基因组学、转录组学分析中的生物信息学问题。 说白了，就是数据的去噪，如何在实验以及分析方法上进行改进。 做单细胞数据分析的人可以先从这篇文章入手，了解单细胞测序存在的问题，并在分析自己的数据时带着这些问题去研究，看看如何检测判断，或者发现、发明新的方法处理。

从单细胞中能看到什么

我们想知道毗连的细胞在遗传和表达层面有什么差异。在人类胚胎干细胞的不同发育阶段细胞之间有什么不同。这些都是单细胞的异质性情况。

现在什么研究中会用到单细胞测序

1. 人类胚胎干细胞分化特点
2. 稀有的转录本分析
3. 肿瘤CTC细胞
4. 肿瘤异质性和微进化
5. 在不好培养的微生物方面研究中

单细胞技术集中的两个问题

1. genome coverage低
2. 扩增有偏好性

第一部分：单细胞DNA测序

Allele dropout

测序数据在基因组覆盖率低会造成SNP脱扣（SNP Dropout）。

在某些文献中MDA的脱扣率高达65%[1]

SNP假阳性

MALBAC假阳率比MDA高。

找寻SNP的算法需要做到什么

1. 需要将SNP和扩增错误分开
2. 能在低覆盖率的数据中找寻SNP

CNV的扩增偏差

在单细胞CNV研究中要适当扩大bin的大小，这样可以减少mapping偏差带来计算不准确（用reads count来算CNV的时候）。

计算CNV的算法要注意什么

MDA方法在计算CNV时，会有序列重复问题，在染色体终端也有问题，GC高的地方也有偏差。（MDA-induced copy number biases were reported to associate with sequence repeats and proximity to chromosome ends, increased GC content and annotated CNVs）

对扩增产物的pairewise comparison可以帮助减少假的CNV。

第二部分： 单细胞RNA测序

可以用来研究CTC表达量。

RSEM这个软件可以计算expression level of TPM(transcripts per million)

在单细胞全转录本扩增时会产生的问题有：

1. 扩增无法得道完整的cDNA片段
2. 转录本的扩增效率不一致
3. 低表达量的转录本难以被检测到

文章中还说道FPKM/RPKM没有考虑转录本间的偏差，所以不适合用在单细胞的计算中。（我没看明白这是为什么）

对于单细胞转录本表达情况的定量分析，文章中提出了两点建议：由于在3’和5’端测序质量不好，所以在做标准化时不要根据转录本的长度计算而是根据覆盖度范围内的长度进行计算；用一些方法（机器学习之类的）研究扩增偏差和正常情况的区别，开发新的工具来减少计算时的扩增偏差。

第三部分：单细胞能研究哪些问题

说白了，还是开头说的那些内容，一个是肿瘤演化谱系，另一个是干细胞发育。

硬件和软件的兼容性太差，系统自动升级后，竟然不能用了，反复在开机界面重启。有钱别用它。

这篇文章比较偏，如果不是做肝癌研究的，会有很多不明白的地方，另外，对于组学数据的分析写得比较少。做的东西只是描述了一些表面现象，没有挖掘机理。

不是专门做肝癌及其组学数据分析，不建议看这篇文章。

文章的分析主要工作是由Zhao Yi老师组完成的。

Identification of prognostic biomarkers in hepatitis B virus-related hepatocellular carcinoma and stratification by integrative multi-omics analysis

上面是文章中的题目，这个题目太长了，我在网页链接里缩写了一个简略的题目（同文章中的题目不完全一致）。

Background

肝癌如何分类？

肝肿瘤分成多种，其中只有几种属癌肿性。最主要的分类法，是鉴定肿瘤属良性（害处相对较少）抑或恶性（可以由肝脏扩散至其他部位，故较为严重）。

良性肿瘤

血管瘤（Hemangioma）是最常见的良性肝肿瘤，是始于胎儿的肝脏异常血管生长。身体状况正常的人中，占10% 以上者肝脏有血管瘤。大部分有血管瘤的人均无任何征状，也不需治疗。但在较为罕见的情况下，血管瘤或会扩大并流血，如出现这情况便要进行手术割除。

肝腺瘤（Hepatic adenomas） 是良性的肝细胞肿瘤，大部分情况均无征状也不需治疗。但如果体积大，或会导致痛楚或失血，在这种情况下肿瘤便需割除。肝腺瘤较常见于女性。部分个案显示，避孕药或怀孕或是触发成因。

肝脏局部结节性增生（Focal nodular hyperplasia，简称FNH） 指数种细胞出现类似肿瘤的增长。虽然属良性，但结节性肿瘤与肝癌并不容易分辨。

恶性肿瘤

成人最常见的原发性肝癌（始于肝脏的癌症）是肝细胞癌（hepatocellular carcinoma，简称HCC），即肝细胞出现癌症。这类癌症有数种增长模式。有些在开始时是单一个肿瘤，然后渐渐增大。到较后期时，癌细胞便会蔓延至肝脏其他部位。

肝癌亦可以在肝脏多个部位增长，并演变为多个肿瘤。这种情况在肝硬化病人身上最普遍。

另一种肝癌是始于胆小管的胆管癌（cholangiocarcinoma）。胆小管是输送胆汁至胆囊的管道。

然而，大部分的肝癌均非始于肝脏，而是由身体其他部位开始的癌症，扩散至肝脏。这类癌症是以癌症开始的部位（原发部位）命名，属于继发性肝癌或转移性癌症。譬如，始于肺部的癌症扩散至肝脏，便称为扩散至肝脏的转移性肺癌。继发性肝癌较原发性肝癌高出 30 倍。

Data information

本文对两个病人的外周血，癌旁，第一位病人的原发灶，卫星灶，门静脉癌拴，第二位病人的左肝区和右肝区的进行了基因组和转录组的测序。第一位病人有肝硬化，年龄是四五十岁，第二位病人没有刚硬化，年龄是七十多岁。第一位病人的肝癌属于低分化，第二位病人的肝癌属于高分化，低分化的肝癌细胞恶性程度高。

Question 1

在本文中研究的是原发性肝癌——肝细胞癌，这种癌症在国内，有很大一部分是由于乙肝、丙肝导致的。原发性肝癌，会有多个病灶，那么这些病灶是源于肿瘤细胞的某个单克隆，还是多克隆？

这个问题可以通过对乙肝病毒同宿主的DNA整合方面来进行研究。

1. 如果，对于一个病人，在不同的肝癌组织中，我们找到的乙肝病毒插入序列都在同一个位置，那么很有可能，这位病人的不同部位的肝癌细胞可能属于同一个单克隆。
2. 反之，乙肝病毒在不同肿瘤组织中，插入到了不同的位置，那么这些肝癌细胞肯定不属于同一个单克隆，也就是说这个病人的肝癌细胞起源具有多克隆性。

根据对乙肝病毒基因整合的研究，文中的第一位病人属于情况1，第二位病人属于情况2。我觉得这只是一种可能性，没准第一个病人，乙肝病毒虽然在不同位点的癌组织中整合在基因组的同一个位置，但是这些癌组织也有可能是在不同或相同时间分别产生的。只研究乙肝病毒的整合情况，能否完善的说明单克隆和多克隆的问题？

Question 2

文章分别看了两位病人的SNV, CNV, structural variation的情况，并做了系统生发树。

本文中认为，第一位病人的肝癌演化路线是normal -> primary cancer -> portal vein tumor -> satellite intrahepatic metastases，第二个病人的肝癌两个病灶更有可能是同时出现的。

但是，我没有理解这个图，从图中看，第二位病人的肝癌演化路线，更像是normal -> left part cancer -> right part cancer。

Question 3

对于转录组的分析，要结合基因组，在基因组上有大量CNV的区域的基因，其表达量同copy number数量正相关。

在转录组部分的研究中，没有什么实质的内容，对于高表达和低表达的基因，做了功能富集图，用的软件是Cytoscape plugin Enrichment Map。我没有理解图里面的小点（基因集合）里的基因有哪些。

Question 4

最后，根据基因表达量的差异和KEGG/BioCarta通路分析，文章中找到了21个基因，这些基因属于细胞周期，p53信号，组氨酸代谢等通路中。 之后又在174个病例中进行了验证，其中6个基因存在普遍的（癌症/正常）差异性，之后他们对这些基因和病人的临床数据整合(乙肝表面抗原，谷丙转氨酶，肿瘤分级等)，进行了分析。 文章中发现TTK基因的高表达和低表达，同病人的无复发存活率、总存活率相关性高，TTK基因高表达的病人在肝癌手术后还需要尽早介入治疗。

分析软件

• Circos
• Cytoscape (plugin Enrichment Map)
• SPSS

还有转录组常用分析软件，具体需要看supplementary。

单词本