CN105483210A

CN105483210A - 一种rna编辑位点的检测方法

Info

Publication number: CN105483210A
Application number: CN201410525846.7A
Authority: CN
Inventors: 刘栋兵; 陈健华; 李欣玥
Original assignee: BGI Shenzhen Co Ltd
Current assignee: BGI Shenzhen Co Ltd
Priority date: 2014-09-30
Filing date: 2014-09-30
Publication date: 2016-04-13

Abstract

本发明提供了一种RNA编辑位点的检测方法，所述方法包括步骤：对群体样本进行测序，获得群体样本的DNA和RNA数据；分析所获得的所述群体样本的DNA和RNA数据，得到RNA编辑位点数据。使用上述方法获得的RNA编辑位点数据更为准确，得到的RNA编辑位点更为全面。使用本方法获得的RNA编辑位点与SNP位点差异显著。

Description

一种RNA编辑位点的检测方法

技术领域

本发明属于生物技术领域，具体地说，本发明涉及一种RNA编辑位点的检测方法。

背景技术

RNA(核糖核酸)编辑是指在RNA水平上发生的突变，包括插入、缺失和单碱基替换等。RNA编辑过程会导致蛋白质编码区的非同义替换、发生可变剪切等，发生RNA编辑可使碱基序列发生改变，甚至可使氨基酸序列乃至蛋白质的高级结构发生改变，它对基因表达调控、基因表达产物的多样化起重要的作用，同时还可能跟某些疾病的致病机理有关。哺乳动物中，发生RNA编辑最常见的是A(adenosine)到I(inosine)突变。

要确定整个基因组范围内的RNA编辑事件，可以通过对RNA样本的高通量测序以及生物信息学分析来实现。随着二代测序技术日益成熟，测序成本不断降低，研究人员可以通过对全转录组进行测序进行RNA编辑方面的研究。当前国内外在此方面的研究主要在个体水平上，对个体水平的RNA编辑研究，可以找到一些跟疾病发病有关的基因靶点，比如跟很多遗传病相关的ADAR基因家族，但是目前由于缺乏标准的比对RNA组，导致很难获得准确的RNA编辑事件。本发明主要涉及一种准确的从RNA高通量测序数据中检测RNA编辑事件的方法。

发明内容

本发明的目的在于提供一种一种RNA编辑位点的检测方法。

本发明提供了一种RNA编辑位点的检测方法，所述方法包括步骤：对群体样本进行测序，获得群体样本的DNA和RNA数据；分析所获得的所述群体样本的DNA和RNA数据，得到RNA编辑位点数据；

其中，所述群体样本中的样本数量≥10个；优选地所述群体样本中样本数量≥20个，更优选地所述群体样本中样本数量≥50个，最优选地所述群体样本中样本数量≥100个。

在另一优选例中，所述样本为正常组织样本和/或肿瘤组织样本。

在另一优选例中，所述样本选自：正常人和/或肿瘤患者。

在另一优选例中，所述方法包括步骤：

(1)比对

(1.1)获得群体样本的高通量DNA和RNA测序的原始数据；

(1.2)对原始数据进行初步过滤得有效数据；

(1.3)将有效数据比对到参考基因组上；

(1.4)校正比对结果的碱基质量值；

(1.5)去除比对结果中的重复序列；

(1.6)分割含N的片段(所述含N的片段是指即包含选择性剪接位点的片段)；

(2)检测突变

检测突变，得原始RNA编辑位点数据；

(3)过滤突变

(3.1)过滤原始RNA编辑位点数据中的假阳性位点；

(3.2)排除基因组上的单核苷酸多态性(SNP)位点；

(3.3)移除所检测RNA的插入或缺失位点左右各30bp内的位点；

(3.4)滤除假编辑位点；

(3.5)滤除FS(Phred-scaledp-valueusingFisher'sexacttesttodetectstrandbias)＞20的位点；

(3.6)滤除基因间区以及处在剪切位点左右2bp内的位点；

获得所述RNA编辑位点数据。

在另一优选例中，所述步骤(1.2)中，初步过滤去除的内容有：

(i)含接头的片段；

(ii)N的比例较高(优选地比例≥10％)的片段；

(iii)低质量片段(优选地，所述低质量片段是指平均质量值小于15的片段)。

在另一优选例中，所述步骤(1.1)中的所述高通量DNA和RNA测序的原始数据由Illumina测序平台产生。

在另一优选例中，所述步骤(1.3)中，比对工具为基于Bowtie的RNA-seq比对工具tophat。

在另一优选例中，所述步骤(1.4)中，使用GATK(GenomeAnalysisToolkit)工具校正比对结果的碱基质量值。

在另一优选例中，所述步骤(1.5)中，利用Picard工具包去除比对结果中存在的重复序列。

在另一优选例中，所述步骤(1.6)中，使用GATK分割比对结果中含N的片段。

在另一优选例中，所述步骤(2)中，使用GATK的UnifiedGenotyper工具检测突变。

在另一优选例中，所述步骤(2)中，分别对样本的正常RNA、肿瘤RNA、正常DNA和肿瘤DNA四组bam文件进行突变检测。

在另一优选例中，所述步骤(3.1)中，使用GATK对检测出来的原始RNA编辑位点数据做VQSR(Variantqualityscorerecalibration)，过滤掉假阳性位点。

在另一优选例中，所述步骤(3.2)中的单核苷酸多态性(SNP)位点包括：基因组上的单核苷酸多态性(SNP)位点与RNA检测的共有位点，以及dbSNP数据库中的SNP位点与RNA检测的共有位点。

在另一优选例中，所述步骤(3.4)中，以深度大于2且突变支持数大于1作为一个可信的发生编辑的样本，若该组样本的编辑样本支持数少于2个，则视为假编辑位点。

应理解，在本发明范围内中，本发明的上述各技术特征和在下文(如实施例)中具体描述的各技术特征之间都可以互相组合，从而构成新的或优选的技术方案。限于篇幅，在此不再一一累述。

附图说明

图1左图和右图分别显示了darned和dbsnp位点的自由能分布情况。

图2显示了darned编辑位点自由能分布。

图3显示了正常组织自由能分布。

图4显示了肿瘤组织自由能分布。

具体实施方式

本发明人通过广泛而深入的研究，获得一种RNA编辑位点的检测方法，实验结果表明，使用所述的方法获得的RNA编辑位点数据更为准确，得到的RNA编辑位点更为全面。使用本方法获得的RNA编辑位点与SNP位点差异显著。

测序

在本发明中，可用常规的测序技术和平台进行测序。优选的测序方法包括：LifeTechnologies的proton或PGM,IlluminaHiSeq,ABISOLiD,Roche454等测序仪器。

在本发明中，特别适合对本发明构建的PCR－free文库进行测序的方法是IonProton法。在一优选例中，将符合上机测序标准的文库片段，使用TheIonProton^TMSystem进行测序。

数据处理

在本发明的优选例中，数据处理通常包括以下步骤：以NCBI数据库中公布的人基因组为参考标准。将测序的reads转换为fastq格式，并与人基因组序列比对，确定匹配的读序(即比对上的读序)。

在另一优选例中，在数据分析之前，还包括对所获得的读序进行长度选择，选取读长为80-220bp，较佳地90-200bp，更佳地100-190bp的读序。

数据处理可以用本领域采用的方法或软件进行，包括市售的软件、公开的软件(尤其是全部开源的软件)进行。

检测方法

本发明主要针对Illumina测序平台生产的高通量测序数据，所包括的检测方法步骤如下：

(1)比对

(1.1)获得原始测序数据，所述原始测序数据为群体样本的测序数据；

在本发明的一个较佳实施方式中，所述原始测序数据包括正常DNA、肿瘤DNA、正常RNA、肿瘤RNA的高通量测序数据；

(1.2)原始数据过滤，目的是过滤掉一些含有接头或者质量值比较低的片段，获得“干净的”数据；主要内容有：

(i)去除含接头的片段；当片段被接头污染时，可能测到接头序列，所以要

除接头；

(ii)去除N的比例较高(优选地比例≥10％)的片段，N含量过高会引起比

对错误；

(iii)去除低质量片段,测序时存在测错的概率，低质量的片段有可能存在

测错的碱基。

(1.3)比对，利用基于Bowtie的RNA-seq比对工具tophat将测序数据比对到参考基因组上，生成bam格式的文件。

(1.4)使用GATK(GenomeAnalysisToolkit)对比对结果的碱基质量值校正。illumina测序结果在给定每个碱基质量值的时候存在偏差，需要根据整个文库所有测序reads的质量值分布进行校正。

(1.5)利用Picard工具包去除比对结果中存在的PCR重复序列。

(1.6)使用GATK分割比对结果中存在剪切的序列(含N的片段)。

(2)使用GATK的UnifiedGenotyper工具检测突变，分别对正常RNA、肿瘤RNA、正常DNA和肿瘤DNA四组bam文件进行突变检测，得到正常RNA、肿瘤RNA、正常DNA和肿瘤DNA一共4个vcf文件，作为原始RNA编辑位点数据(原始SNP)。

(3)过滤突变

(3.1)使用GATK对检测出来的SNP做VQSR(Variantqualityscorerecalibration)，对vcf(VariantCallFormat)文件中的一些高质量的位点作为可信位点构建高斯混合模型(Gaussianmixturemodel)，并对所有位点进行评估，从而过滤其中的假阳性位点，具体操作可以参考软件说明。

(3.2)分别移除DNA和RNA、RNA和dbSNP数据库共有的位点，因这些位点并不是在转录过程中发生的突变，不属于RNA编辑事件，需要排除。

(3.3)移除RNA检测的indel(插入或缺失)位点左右各30bp(basepair)内的位点，由于indel附近容易发生比对错误，造成较高的假阳性，因此将INDEL附近的位点排除掉。

(3.4)以深度大于2且突变支持数大于1作为一个可信的发生编辑的样本，若该组样本的编辑样本支持数少于2个，则视为假编辑位点滤掉。

(3.5)过滤掉FS(Phred-scaledp-valueusingFisher'sexacttesttodetectstrandbias)大于20的位点。

(3.6)移除基因间区以及处在剪切位点左右2bp内的的位点，由于处在这些区域的位点的突变并不会对基因表达产物产生直接影响，因此也需要过滤掉。

最终得到高质量的在基因区的RNA编辑位点。

其中，步骤1.3)中基于公开的开源Bowtie的RNA-seq比对工具tophat(下载地址如：http://ccb.jhu.edu/software/tophat/index.shtml)，进行比对，命令行如下：

tophat--solexa1.3-quals--read-mismatches2--read-gap-length3--read-edit-dist3--library-typefr-unstranded-p6-r30--b2-fast--rg-centerbgi--rg-platformillumina--no-novel-juncs--no-novel-indels-odirreferencesequence.fq1sequence.fq2

步骤1.4)中，使用公开的开源GATK(GenomeAnalysisToolkit)软件(下载地址如：https://www.broadinstitute.org/gatk/)，校正参数为-knownSites-nct-U–BQSR，GATK的软件使用可以参考产品使用说明。

步骤1.5)中，公开的开源Picard工具包(下载地址如：http://picard.sourceforge.net/)去除比对结果中存在的PCR重复序列，设置如下：

java-Xmx4g-jarMarkDuplicates.jarINPUT＝in.bamOUTPUT＝out.bamMETRICS_FILE＝rmdup.metREMOVE_DUPLICATES＝trueVALIDATION_STRINGENCY＝SILENTASSUME_SORTED＝trueCREATE_INDEX＝true。

步骤1.6)中，GATK的设置如下：

java-Xmx512M-jarFilterBadCigar.jarin.bamout.bamjava-Xmx6g-jarGenomeAnalysisTK.jar-TSplitNCigarReads-Iin.bam-oout.bam-UALL–Rreference.fa。

步骤(1.2)中，GATK的UnifiedGenotyper工具的设置如下：

java-Xmx6g-jar-Djava.io.tmpdir＝tmpGenomeAnalysisTK.jar-TUnifiedGenotyper-lINFO-Ibam.list-Rreference.fa--dbsnpdbsnp_138-stand_call_conf30-stand_emit_conf4-dcov200-GStandard-nt6-glmBOTH-UALLOW_N_CIGAR_READS-Lchr-metricsmetrics-ochr.vcf

步骤3.1)中，VQSR(Variantqualityscorerecalibration)是指:对vcf(VariantCallFormat)文件中的一些高质量的位点作为可信位点构建高斯混合模型(Gaussianmixturemodel)，并对所有位点进行评估，从而过滤其中的假阳性位点；

主要步骤：(i)对vcf(VariantCallFormat)文件中的一些高质量的位点作为可信位点构建高斯混合模型(Gaussianmixturemodel)，并对所有位点进行评估；(ii)将建立的高斯混合模型参数应用到输入的VCF文件，对每一个变异位点进行VQSLOD值的注释，从而过滤其中的假阳性位点。

VQSR通过机器学习的方法根据已知的变异位点训练出一组变异位点集，并会给每个位点赋一个VQSLOD值，变异位点越接近集合的中心其值就会越高；然后根据模型在对新检测出的变异位点进行打分，如果分值在训练集合内就认为是一个质量高的变异位点，否则认为是一个假阳性位点。

步骤(3.5)，FS(Phred-scaledp-valueusingFisher'sexacttesttodetectstrandbias)使用Fish检验的方法，检测比对在某一位点片段是否存在链的偏好性。

与现有的RNA编辑研究方法相比，本发明结合新一代测序技术和生物信息学分析方法，可同时从多个样本中检测RNA编辑事件，并且超过90％的位点处在Alu(一种转座子，是哺乳动物基因组中SINE家族的一员，每个长约300bp，单个成员的每个末端上有Alu切割点区)，在Alu区95％左右的位点属于“A->G”突变，突变的类型及相关比例都跟单个样本的检测结果相当吻合，说明本分析流程的检测结果是比较准确的。

本发明的主要优点在于：

(1)首次给出了一种一种RNA编辑位点的检测方法，采用该方法能够简单便捷的处理高通量的核酸数据。

(2)本发明的方法检测结果准确，能够有效的区分SNP位点与RNA编辑位点。

(3)能够对群体的RNA编辑事件进行检测。

下面结合具体实施例，进一步陈述本发明。应理解，这些实施例仅用于说明本发明而不用于限制本发明的范围。下列实施例中未注明详细条件的实验方法，通常按照常规条件如Sambrook等人，分子克隆：实验室手册(NewYork:ColdSpringHarborLaboratoryPress,1989)中所述的条件，或按照制造厂商所建议的条件。除非另外说明，否则百分比和份数按重量计算。

实施例1

为了使本发明的用法和效果更加易于理解和掌握，下面将举一个实例进行进一步的阐述。

1.样本/数据来源

本实施例中群体样本为65例前列腺癌患者，对每个患者的正常DNA、肿瘤DNA、正常RNA、肿瘤RNA分别进行高通量测序，读长为90bp。

2.分析步骤

2.1下机数据过滤，过滤条件为,过滤掉含“N”比例大于10％和平均质量值小于15的片段，过滤掉质量值小于10的碱基比例大于40％的片段。

2.2比对及突变检测

a)利用tophat将测序数据比对到人类基因组(hg19)上,参数为：

--read-mismatches2--read-gap-length3--read-edit-dist3-r30--no-novel-juncs--no-novel-indels

其中“--read-mismatches”为片段容许的最大错配数；“--read-gap-length”为比对结果中容许开间隙的最大总长度；“--read-edit-dist”为测序片段对参考基因组的最大编辑距离；“-r”为双末端测序的两条片段的平均距离；“--no-novel-juncs”表示不预测新的剪切位点，“--no-novel-indels”表示不检测indel，这两个参数主要是为了减少程序运行的时间；其他参数设置为默认。

b)对于测了多个单泳道(lane)的样本，需要用picard将它们的比对结果合并成一个文件,随后给每个合并的bam文件添加相应的样本信息并做去重处理。

c)利用DATK的“BaseRecalibrator”对每个bam文件做处理，得到包含校正信息的表格文件，再用“PrintReads”工具结合该文件可得到校正后的bam文件，最后用“SplitNCigarReads”将bam文件中含有“N”的片段分割开。

d)用GATK的“UnifiedGenotyper”分别对正常RNA、肿瘤RNA、正常DNA和肿瘤DNA四组bam文件进行突变检测，得到正常RNA、肿瘤RNA、正常DNA和肿瘤DNA一共4个vcf文件。

2.3突变过滤

根据标准操作流程，对正常和肿瘤RNA的vcf文件分别进行过滤，每一步得到的数目统计如下表：

步骤	正常组织	肿瘤组织
			原始SNP	10825607	10417648

VQSR	6230390	6255138
			过滤DNA和dbsnp	1667490	1691491
过滤INDEL附近位点	1594692	1627980
			样本支持数过滤	477545	460967
FS过滤	470670	455615
			过滤基因间区/剪切位点	399385	379833

3.编辑位点评估

3.1A->G以及Alu区位点统计

从上表看，这组数据跟之前的研究结果(RamaswamiGetal.Naturemethods,579–583(2012))很接近，说明结果是比较准确的。

3.2利用RNAfold对本次分析结果、dbsnp和darned数据库的位点进行二级结构自由能预测并比较，可看到dbSNP(寡核苷酸多态性数据库，网址：(http://hgdownload.cse.ucsc.edu/goldenPath/hg19/database/snp138.txt.gz))与darned(RNA编辑位点数据库，网址：http://darned.ucc.ie/)的自由能分别有较大区别(图1)，本发明人的分析结果则跟数据库结果较接近(图2-4)，说明检测结果准确性较好，能够有效的区分SNP位点与RNA编辑位点。

群体RNA编辑检测，对于每个位点，可以综合利用所有样本在该位置的测序信息，相比较单样本而言，能有效的鉴别可能因单样本数据量少而造成的漏检情况。

在本发明提及的所有文献都在本申请中引用作为参考，就如同每一篇文献被单独引用作为参考那样。此外应理解，在阅读了本发明的上述讲授内容之后，本领域技术人员可以对本发明作各种改动或修改，这些等价形式同样落于本申请所附权利要求书所限定的范围。

参考文献

1.McKennaA,HannaM,BanksE,etal.TheGenomeAnalysisToolkit:aMapReduceframeworkforanalyzingnext-generationDNAsequencingdata[J].Genomeresearch,2010,20(9):1297-1303.

2.TrapnellC,WilliamsBA,PerteaG,etal.TranscriptassemblyandquantificationbyRNA-Seqrevealsunannotatedtranscriptsandisoformswitchingduringcelldifferentiation[J].Naturebiotechnology,2010,28(5):511-515.

3.TrapnellC,PachterL,SalzbergSL.TopHat:discoveringsplicejunctionswithRNA-Seq[J].Bioinformatics,2009,25(9):1105-1111.

4.LangmeadB,SalzbergSL.Fastgapped-readalignmentwithBowtie2[J].Naturemethods,2012,9(4):357-359.

5.LiH,HandsakerB,WysokerA,etal.Thesequencealignment/mapformatandSAMtools[J].Bioinformatics,2009,25(16):2078-2079.

6.RamaswamiG,LinW,PiskolR,etal.AccurateidentificationofhumanAluandnon-AluRNAeditingsites[J].Naturemethods,2012,9(6):579-581.

7.PengZ,ChengY,TanBCM,etal.ComprehensiveanalysisofRNA-SeqdatarevealsextensiveRNAeditinginahumantranscriptome[J].Naturebiotechnology,2012,30(3):253-260.

8.MaasS,KawaharaY,TamburroKM,etal.A-to-IRNAeditingandhumandisease[J].RNABiol,2006,3(1):1-9.

Claims

1.一种RNA编辑位点的检测方法，其特征在于，所述方法包括步骤：对群体样本进行测序，获得群体样本的DNA和RNA数据；分析所获得的所述群体样本的DNA和RNA数据，得到RNA编辑位点数据；

2.如权利要求1所述的方法，其特征在于，所述样本为正常组织样本和/或肿瘤组织样本。

3.如权利要求1所述的方法，其特征在于，所述方法包括步骤：

(1)比对

(1.1)获得群体样本的高通量DNA和RNA测序的原始数据；

(1.2)对原始数据进行初步过滤得有效数据；

(1.3)将有效数据比对到参考基因组上；

(1.4)校正比对结果的碱基质量值；

(1.5)去除比对结果中的重复序列；

(1.6)分割含N的片段；

(2)检测突变

检测突变，得原始RNA编辑位点数据；

(3)过滤突变

(3.1)过滤原始RNA编辑位点数据中的假阳性位点；

(3.2)排除基因组上的单核苷酸多态性(SNP)位点；

(3.3)移除所检测RNA的插入或缺失位点左右各30bp内的位点；

(3.4)滤除假编辑位点；

(3.6)滤除基因间区以及处在剪切位点左右2bp内的位点；

获得所述RNA编辑位点数据。

4.如权利要求3所述的方法，其特征在于，所述步骤(1.2)中，初步过滤去除的内容有：

(i)含接头的片段；

(ii)N的比例较高(优选地比例≥10％)的片段；

5.如权利要求3所述的方法，其特征在于，所述步骤(1.1)中的所述高通量DNA和RNA测序的原始数据由Illumina测序平台产生；和/或

所述步骤(1.3)中，比对工具为基于Bowtie的RNA-seq比对工具tophat；和/或

所述步骤(1.4)中，使用GATK(GenomeAnalysisToolkit)工具校正比对结果的碱基质量值；和/或

所述步骤(1.5)中，利用Picard工具包去除比对结果中存在的重复序列；和/或

所述步骤(1.6)中，使用GATK分割比对结果中含N的片段。

6.如权利要求3所述的方法，其特征在于，所述步骤(2)中，使用GATK的UnifiedGenotyper工具检测突变。

7.如权利要求6所述的方法，其特征在于，所述步骤(2)中，分别对样本的正常RNA、肿瘤RNA、正常DNA和肿瘤DNA四组bam文件进行突变检测。

8.如权利要求3所述的方法，其特征在于，所述步骤(3.1)中，使用GATK对检测出来的原始RNA编辑位点数据做VQSR(Variantqualityscorerecalibration)过滤掉假阳性位点。

9.如权利要求3所述的方法，其特征在于，所述步骤(3.2)中的单核苷酸多态性(SNP)位点包括：基因组上的单核苷酸多态性(SNP)位点与RNA检测的共有位点，以及dbSNP数据库中的SNP位点与RNA检测的共有位点。

10.如权利要求3所述的方法，其特征在于，所述步骤(3.4)中，以深度大于2且突变支持数大于1作为一个可信的发生编辑的样本，若该组样本的编辑样本支持数少于2个，则视为假编辑位点。