这个文章是华大做的，华大从12年开始连续发了好几篇关于肿瘤单细胞相关的文章，他们用的测序方法是MDA， 几篇文章都主要围绕测序方法的准确性、寻找肿瘤中的driver（驱动）基因、以及重要突变、单克隆多克隆这几个问题展开讨论。

在这篇文章中，他们主要发现了结肠癌的双克隆起源，以及一个新的致癌基因（就是标题）。 文中发现肿瘤主克隆中在致癌基因APC,TP53基因上有突变，而在另一个克隆上有CDC27和PABPC1的突变。 他们还发现SLC12A5在单细胞中有高突变率(frequency of mutation)而在群体中的普遍程度却比较低，在实验验证后发现这个基因是一个潜在的致癌基因。

从这篇文章中学习一下文章的结构。

第一部分：说明测序方法的可行性（准确性）

• 检测乘数，覆盖度
• 检测假阳性，假阴性
• 检测组织样本和单细胞样本的突变频率相关性

第二部分：肿瘤细胞可分为两个群体

• PCA 分析区分group
• 基于UPGMA的phylogenetic tree分析
• 随机模拟单克隆到5个克隆的情况，发现实际的频率分布同双克隆相似

第三部分：从两个分组中发现了不同的驱动事件

• 引入21个结肠癌病人数据，发现有overlap的驱动基因可以将单细胞分成三组（正常，克隆1，克隆2）。
• 主克隆里有APC，TP53，副克隆里有CDC27和PABPC1。

第四部分：细胞起源异质性分析

• 用21个样本检测是不是多克隆是个普遍的现象。

第五部分：检测癌症驱动基因

• 检测非同义突变、同义突变、错义突变、无义突变在一些重要癌症基因的出现频率。
• 找出enrichment的基因，这些基因是驱动基因。

第六部分：对于新致癌基因的功能验证

• 做了PCR看这个基因是不是在肠癌细胞系中表达水平高。

材料与方法部分有很多信息

1. 找somatic mutation的原则。
2. 如何检测allele dropout。
3. 在研究异质性问题是，所用的SNV要排除那些在CNV异常区域的。

单词本

Adevances and Applications of Single-Cell Sequencing Technologies 这篇review很值得一读，它总结了现有的单细胞领域技术、研究内容、现有问题以及今后方向。 下面简单总结一下比较感兴趣的内容。

单细胞DNA测序技术的缺陷包括: 测序覆盖度不均匀、allelic dropout（等位片段测序丢失）、假阳性、假阴性等。

multiple-displacement-amplification(MDA)方法由于测序覆盖不均匀，在CNV研究中不值得采用这种方法。

multiple annealilng- and looping-based amplification cycles(MALBAC) 方法假阳性高，但是研究CNV比较好。

NUC-SEQ/single nucleus exome sequencing(SNES) 方法可减少技术误差。

单细胞RNA测序技术

WAT会造成3’mRNA 的偏差。

区分技术误差与生物学多样性(Biological Variation)

为了解释一个突变或者转录物的确在样本里时杂合的，需要对这些结果做必要的实验验证。 在DNA层面用深度测序，digital droplet PCR（微滴式数字PCR）。 在RNA层面，可用单细胞RT-qPCR以及微滴式数字PCR。

Tissue Mosaicism的研究结果争议性很大（可以参考上一篇文章）

总结下来就是每个研究中都可以看到copy number mosaicism现象，但是这些现象的作用时正向还是负向目前还尚未有定论。

癌症方面的研究

研究主要集中在瘤内异质性、单/多克隆演化、CTC细胞携带mutation特点(肿瘤细胞转移问题)。

研究CTC细胞是为了追溯肿瘤和转移灶的mutation情况。

###小结

单细胞研究的问题涵盖：

1. 解读群体差异性
2. 追踪细胞系
3. 细胞类型分类
4. 罕见细胞的基因组学信息刻画

###展望

1. 从单细胞入手，研究单细胞其实是研究单细胞所在的组织的情况。
2. 如何连接并研究单细胞的phenotype到genotype。
3. 如何研究单细胞的多层面数据（例如，某个单细胞的DNA和RNA，某个单细胞的RNA和表观）
4. 在技术上，如何能同时测序上千个单细胞，并降低费用。

这篇去年在PNAS上的文章用Sigma GenomePlex Single Cell Whole Genome Amplification Kit进行单细胞测序，发现异倍体在正常细胞组织中出现的不多。 这个结论同之前的多项研究中说明异倍体是正向选择的结果正好相反。文章只有短短的6页，推翻了一个之前在多篇文章中都有涉及的结论。 整个内容就是找不同的正常人的肝脏和脑神经细胞，进行单细胞测序，计算CNV（500kb一个窗口），看是不是有CNV的变化。

有用的知识

1. 异倍体在正常组织和癌组织中都有

2. 在人类BUB1B基因上的突变回导致 mosaic variegated aneuploidy(MVA，斑花叶非整倍体综合征)，BUB1B编码 BubR1蛋白，这个蛋白在纺锤体assembly checkpoint过程中染色质分离时起到作用。 MVA疾病的病人在一条或多条染色体上会出现异倍体的情况。MVA病人发育迟缓，会在儿童时期得癌症。

3. 前人用的FISH技术和spectral karyotyping (SKY)技术都过高估计了异倍体的出现，所以用single cell来解决这个问题（一个一个细胞测得准）。

单词本

merotelic:一个着丝粒同时收到来自相反方向的纺锤丝牵引

Tang Fuchou 老师实验室做胚胎发育分析以及单细胞RNA-seq非常有经验，最近（7月23日）他们发表了一个新单细胞转录组测序方法———— single-cell universal poly(A)-independent RNA sequencing (SUPeR-seq)，主要特点是在cDNA合成时用有固定anchor序列的随机primer代替传统的oligo(dT) primer。

对于实验部分我理解的不多，主要看信息分析部分有什么关于线性、环状RNA的新结果。

在小鼠胚胎细胞不同发育阶段的研究中：

• 该研究中大多数circRNA含有多个外显子，少数（9%）circRNA只含有一个外显子。
• 只有0.9%的circRNA有多于20个潜在的miRNA绑定位点，在功能上绝大多数可能没有起到miRNA海绵的作用。
• circRNA在卵母细胞阶段就已经可以检测到，到受精卵四个-八个细胞阶段开始下降，囊胚细胞阶段降得很低(为啥？)。
• 出现circRNA多的基因的线性转录本表达量也高。
• circRNA两侧的intron要比不环化的RNA intron长。
• 一个基因可以产生多个circRNA，另外这些circRNA倾向于使用同一个5’端的exon，而3’端的exon差异大。
• 对于使用同一个5’端exon的circRNAs，在下游有长intron的circRNA的数量相对同个基因产生的其他circRNA来说显著的多。反之，使用同一个3’端的exon，上游有长intron的circRNA更富集。
• repeat元件在circRNA两侧的intron中不明显富集，但由于circRNA两侧的intron长度相对较长，所以repeat序列的总数要比其他的intron多。（从这点来看repeat序列其作用不在于排列的密度，在于绝对数量）
• 有互补序列的circRNA数量更加富集。
• 相对距离5kb左右的互补序列才有作用。
• 有强splicing motif时不容易形成circRNA。
• circRNA的上下游motif对circRNA的表达起到相反的作用，上游motif强，表达的circRNA多，下游motif弱，表达的circRNA多。

另外，我认为下面这句话看着不太地道。

Gene ontology(GO) analysis of all the 1316 host genes producing these circRNA transcripts showed strong enrichment for terms related to chromatin organization, cell division, and response to DNA damage stimulus, suggesting potential roles of these circRNAs in these functional areas.

最近在看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 *

（未完待续）