小鼠胚胎单细胞RNA-seq线性、环状RNA分析
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 检测工具介绍
最近在看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的短序列。
需要检测的条件如下:
- GU/AG 在剪接位点的两侧出现
- 可以检测到清晰的breakpoint
- 只支持2个mismatch
- breakpoint不能再短序列(anchor)以里2 nucleotides 之外的地方出现,也就是最多2nt (这条可能理解有偏差)
- 至少有两条reads支持这个junction
- 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所遇到的问题:
- circRNA比其他RNA在细胞中的比例低,一般RNA-seq的实验步骤中,不包括circRNA富集的步骤(例如RNase R treatment,我觉得他们要表到的意思是在RNA-seq中能探测到的circRNA比用RNase R处理过的数据中少),所以在RNA-seq中的circRNA比较低,假阳性也高。
- 目前的注释文件都是根据线性RNA进行的,所以也不适合circRNA的识别,另外非模式动物的注释信息更少,就更不用谈在那些物种里找circRNA了。
- 由于RNA测序数据的reads长度差别大,也对检测circRNA工作带来了不便
- 套索结构和融合基因同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在新英格兰医学杂志上安利了一下CRISPR
标题是哗众取宠,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老师组在去年年初总结了一下目前单细胞基因组学、转录组学分析中的生物信息学问题。 说白了,就是数据的去噪,如何在实验以及分析方法上进行改进。 做单细胞数据分析的人可以先从这篇文章入手,了解单细胞测序存在的问题,并在分析自己的数据时带着这些问题去研究,看看如何检测判断,或者发现、发明新的方法处理。
从单细胞中能看到什么
我们想知道毗连的细胞在遗传和表达层面有什么差异。在人类胚胎干细胞的不同发育阶段细胞之间有什么不同。这些都是单细胞的异质性情况。
现在什么研究中会用到单细胞测序
- 人类胚胎干细胞分化特点
- 稀有的转录本分析
- 肿瘤CTC细胞
- 肿瘤异质性和微进化
- 在不好培养的微生物方面研究中
单细胞技术集中的两个问题
- genome coverage低
- 扩增有偏好性
第一部分:单细胞DNA测序
Allele dropout
测序数据在基因组覆盖率低会造成SNP脱扣(SNP Dropout)。
在某些文献中MDA的脱扣率高达65%[1]
SNP假阳性
MALBAC假阳率比MDA高。
找寻SNP的算法需要做到什么
- 需要将SNP和扩增错误分开
- 能在低覆盖率的数据中找寻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)
在单细胞全转录本扩增时会产生的问题有:
- 扩增无法得道完整的cDNA片段
- 转录本的扩增效率不一致
- 低表达量的转录本难以被检测到
文章中还说道FPKM/RPKM没有考虑转录本间的偏差,所以不适合用在单细胞的计算中。(我没看明白这是为什么)
对于单细胞转录本表达情况的定量分析,文章中提出了两点建议:由于在3’和5’端测序质量不好,所以在做标准化时不要根据转录本的长度计算而是根据覆盖度范围内的长度进行计算;用一些方法(机器学习之类的)研究扩增偏差和正常情况的区别,开发新的工具来减少计算时的扩增偏差。
第三部分:单细胞能研究哪些问题
说白了,还是开头说的那些内容,一个是肿瘤演化谱系,另一个是干细胞发育。
整个周末都耗费在修理小米手机上了
硬件和软件的兼容性太差,系统自动升级后,竟然不能用了,反复在开机界面重启。有钱别用它。