CN111192636A - 一种适用于oligodT富集的mRNA二代测序结果分析方法 - Google Patents

Abstract

本发明公开了一种适用于oligodT富集的mRNA二代测序结果分析方法，其特征在于，包括测序数据过滤步骤、测序数据比对步骤、样品表达量计算步骤、差异分析步骤、富集分析步骤、外显子差异分析步骤、可变剪切分析步骤、序列差异分析步骤和转录本序列制备步骤。本发明的方法更加快速和适用于ligodT富集的mRNA二代测序。

Description

一种适用于oligodT富集的mRNA二代测序结果分析方法

技术领域

本发明涉及基因检测领域，具体涉及一种适用于oligodT富集的mRNA二代测序结果分析方法。

背景技术

转录组是某个物种或者特定细胞类型产生的所有转录本的集合。转录组研究能够从整体水平研究基因功能以及基因结构，揭示特定生物学过程以及疾病发生过程中的分子机理，已广泛应用于基础研究、临床诊断和药物研发等领域。真核生物的蛋白编码基因在3’末端有一段poly(A)尾，所以对于真核生物，提取总RNA后，可以用带有Oligo(dT)的磁珠富集mRNA，打断成为短片段，再以片断后的mRNA为模板，制备RNASeq文库，构建好的文库用Illumina测序仪进行测序。

测序后的数据需要进行生物信息学分析，获取样品的基因表达信息，推断生物学意义。通常一个样品可以获得数千万个测序reads，之前的分析方法存在计算机资源消耗大、运行时间慢等缺点。同时，不断的出现新的分析方法和软件，现有的mRNA分析流程需要优化和补充。

发明内容

为了克服现有技术的上述缺陷，本发明的目的在于提供一种适用于oligodT富集的mRNA二代测序结果分析方法。

为了实现本发明的目的，所采用的技术方案是：

一种适用于oligodT富集的mRNA二代测序结果分析方法，包括如下步骤：

步骤一，测序数据过滤步骤：

使用fastp软件去除测序结果中的接头和低质量序列：

fastp使用PE reads overlap信息自动识别接头序列，具体是：

同时以5个碱基长度为窗口，从3’端向5’端滑动，截去窗口内平均质量小于20的窗口。最后保留长度大于50的reads；

将去除接头和低质量序列后的测序数据采用fastqc软件进行质量控制：

所述质量控制为根据RNASeq的测序特点提供质量控制标准；

步骤二，测序数据比对步骤：

使用hisat2软件将过滤后的数据与参考基因组进行比对：

首先从基因组数据库中下载对应物种的参考基因组，下载下来的基因组序列使用hisat2-build构建索引；

对于构建索引后的序列进行比对；

将比对完成的测序数据通过RSeQC进行质量控制；

将比对完成的测序数据通过RSeQC进行质量控制：

所述RSeQC的质量控制为，所述RSeQC结果的reads分布在基因间区上呈现出中间高，两端低的分布形状；

步骤三，样品表达量计算步骤：

使用htseq计算样品的表达量，具体计算为：

Htseq结合测序reads在基因组上的比对结果和基因在基因组上的位置，将reads分配到各个基因上。首先需要给比对的结果进行排序，运行htseq的时候需要指定文库的链特异性信息-s yes/reverse/no。htseq默认只保留唯一比对或只有对应一个基因的reads，判断可能属于多个基因的reads将被舍弃；

步骤四，差异分析步骤：

使用DESeq软件对表达后的数据进行差异分析，具体是：

将表达矩阵导入到DESeq包中进行差异分析，保留|log2FC|>1,pvalue<0.05的作为差异表达的转录本；

步骤五，富集分析步骤：

使用topGO软件包进行GO富集分析，具体是：

根据KEGG注释信息，使用phyper进行超几何检验，计算KEGG富集显著性。获得的显著性P值，使用p.adjust进行多重校正，得到校正后的P值，通常选择P<0.05的结果为最后的富集结果；

步骤六，外显子差异分析步骤：

用DEXSeq软件包进行外显子差异分析；

步骤七，可变剪切分析步骤：

使用rMATS进行差异可变剪切分析，获取差异可变剪切数据，具体是：

通过统计模型对不同样本进行可变剪切事件的定量，然后计算P value来表示两组样品在IncLevel(Inclusion Level)水平上的差异；

然后对p value进行校正得FDR值。rMATS可识别的可变剪切事件有5种，分别是skipped exon(SE)外显子跳跃，alternative 5’splice site(A5SS)第一个外显子可变剪切，alternative 3’splice site(A3SS)最后一个外显子可变剪切，mutually exclusiveexons(MXE)外显子选择性跳跃和retained intron(RI)内含子滞留；

步骤八，序列差异分析步骤：

使用varscan软件分析SNP信息，得到样品与参考基因的序列差异；

步骤九，转录本序列制备步骤：

使用stringtie进行转录拼接，重构样品的转录本序列，具体为：

对每个样品单独使用stringtie拼接，使用stringtie merge将所有样品的转录本进行合并。合并后的转录本为最后重构出来的转录本序列。

在本发明的一个优选实施例中，所述根据RNASeq的测序特点的质量控制标准为：

所述测序数据的碱基分布的四条线为：AT平行且接近或GC平行且接近或GC整体分布应该近似于正态分布或不能出现多峰；

测序数据的Duplicate水平与建库的PCR循环数一致。

在本发明的一个优选实施例中，所述构建索引为在hisat2比对时加上链特异性参数，所述链特异性参数当中，当采用dUTP建库方法时，对应的参数为--rna-strandness FR，其他方法当中参数采用默认值。

在本发明的一个优选实施例中，所述步骤二当中的比对为：

模式生物的比对率为不低于80％，多重比对率为不高于10％；如果比对率低，通常需要检查测序样品是否有污染；如果多重比对率高，通常需要检查测序数据的核糖体比例。

在本发明的一个优选实施例中，所述可变剪切分析当中：默认保留差异大于0.01％的可变剪切。

在本发明的一个优选实施例中，所述序列差异分析具体是：

首先需要使用samtools mpileup提取样品在基因组上每一个碱基位置的比对情况；

然后varscan读取每个位置的信息，从中识别其中的错配和插入缺失，并将结果转换成VCF格式输出；

每个SNP需要满足总测序深度大于8；

杂合位点的SNP频率以40％-60进行筛选。

本发明的有益效果在于：

本发明的方法更加快速和适用于ligodT富集的mRNA二代测序。

附图说明

图1为本发明的流程图。

图2为本发明的FASTQC碱基质量分布示意图。

图3为本发明的FASTQC碱基分布示意图。

图4为本发明的表达量分布示意图。

图5为本发明的样品相关性检验示意图。

图6为本发明的差异分析结果示意图。

图7为本发明的差异基因circos示意图。

图8为本发明的富集分析示意图。

图9为可变剪切数量的统计图。

图10为本发明的SNP分布示意图。

具体实施方式

一种适用于oligodT富集的mRNA二代测序结果分析方法，包括如下步骤：

步骤一，测序数据过滤步骤：

使用fastp软件去除测序结果中的接头和低质量序列：

使用fastp软件去除测序结果中的接头和低质量序列：fastp使用PE readsoverlap信息自动识别接头序列，准确性和去除的效率更高。同时以5个碱基长度为窗口，从3’端向5’端滑动，截去窗口内平均质量小于20的窗口。最后保留长度大于50的reads。

将去除接头和低质量序列后的测序数据采用fastqc软件进行质量控制：

根据RNASeq的测序特点，测序数据的碱基分布的四条线应该为：AT平行且接近，GC平行且接近；GC整体分布应该近似于正态分布，不能出现多峰；测序数据的Duplicate水平应该与建库的PCR循环数一致。

步骤二，测序数据比对步骤：

使用hisat2软件将过滤后的数据与参考基因组进行比对：

首先需要从基因组数据库中下载对应物种的参考基因组，常用的基因组数据库为Ensembl和NCBI，下载下来的基因组序列使用hisat2-build构建索引。如果使用的是链特异性文库，hisat2比对时需要加上链特异性参数，常用dUTP建库方法对于的参数为：--rna-strandness FR，其他参数通常使用默认值。

常见的模式生物的比对率通常在80％以上，多重比对率在10％以下。如果比对率低，通常需要检查测序样品是否有污染；如果多重比对率高，通常需要检查测序数据的核糖体比例。

将比对完成的测序数据通过RSeQC进行质量控制：

对于lncRNA项目，因为很多lncRNA都还没有注释信息，所以RSeQC结果的reads分布，在基因间区上需要有reads分布。Reads在基因结构上需要呈现出中间高，两端低的分布形状。

步骤三，样品表达量计算步骤：

使用htseq计算样品的表达量：

Htseq结合测序reads在基因组上的比对结果和基因在基因组上的位置，将reads分配到各个基因上。首先需要给比对的结果进行排序，运行htseq的时候需要指定文库的链特异性信息-s yes/reverse/no。htseq默认只保留唯一比对、只有对应一个基因的reads，可能属于多个基因的reads将被舍弃。

步骤四，差异分析步骤：

使用DESeq软件对表达后的数据进行差异分析：

表达矩阵导入到DESeq包中进行差异分析，保留|log2FC|>1,pvalue<0.05的作为差异表达的转录本。

步骤五，富集分析步骤：

使用topGO软件包进行GO富集分析；根据KEGG注释信息，使用phyper进行超几何检验，计算KEGG富集显著性。获得的显著性P值，使用p.adjust进行多重校正，得到校正后的P值，通常选择P<0.05的结果为最后的富集结果。

步骤六，外显子差异分析步骤：

用DEXSeq软件包进行外显子差异分析。

步骤七，可变剪切分析步骤：

使用rMATS进行差异可变剪切分析，获取差异可变剪切数据；

rMATS可以对RNA-Seq数据进行差异可变剪切分析的软件。其通过统计模型对不同样本(有生物学重复的)进行可变剪切事件的定量，然后计算P value来表示两组样品在IncLevel(Inclusion Level)水平上的差异，然后对p value进行校正得FDR值。rMATS可识别的可变剪切事件有5种，分别是skipped exon(SE)外显子跳跃，alternative 5’splicesite(A5SS)第一个外显子可变剪切，alternative 3’splice site(A3SS)最后一个外显子可变剪切，mutually exclusive exons(MXE)外显子选择性跳跃和retained intron(RI)内含子滞留。默认保留差异大于0.01％的可变剪切。

步骤八，序列差异分析步骤：

使用varscan软件分析SNP信息，得到样品与参考基因的序列差异；

首先需要使用samtools mpileup提取样品在基因组上每一个碱基位置的比对情况，然后varscan读取每个位置的信息，识别其中的错配和插入缺失，并将结果转换成VCF格式输出。每个SNP需要满足总测序深度大于8。杂合位点的SNP频率以40％-60进行筛选。

步骤九，转录本序列制备步骤：

使用stringtie进行转录拼接，重构样品的转录本序列。

首先需要对每个样品单独使用stringtie拼接，然后使用stringtie merge将所有样品的转录本进行合并。合并后的转录本为最后重构出来的转录本序列。

Claims (5)

1.一种适用于oligodT富集的mRNA二代测序结果分析方法，其特征在于，包括如下步骤：

步骤一，测序数据过滤步骤：

使用fastp软件去除测序结果中的接头和低质量序列：

fastp使用PEreads overlap信息自动识别接头序列，具体是：

同时以5个碱基长度为窗口，从3’端向5’端滑动，截去窗口内平均质量小于20的窗口。最后保留长度大于50的reads；

将去除接头和低质量序列后的测序数据采用fastqc软件进行质量控制：

所述质量控制为根据RNASeq的测序特点提供质量控制标准；

步骤二，测序数据比对步骤：

使用hisat2软件将过滤后的数据与参考基因组进行比对：

首先从基因组数据库中下载对应物种的参考基因组，下载下来的基因组序列使用hisat2-build构建索引；

对于构建索引后的序列进行比对；

将比对完成的测序数据通过RSeQC进行质量控制；

将比对完成的测序数据通过RSeQC进行质量控制：

所述RSeQC的质量控制为，所述RSeQC结果的reads分布在基因间区上呈现出中间高，两端低的分布形状；

步骤三，样品表达量计算步骤：

使用htseq计算样品的表达量，具体计算为：

Htseq结合测序reads在基因组上的比对结果和基因在基因组上的位置，将reads分配到各个基因上。首先需要给比对的结果进行排序，运行htseq的时候需要指定文库的链特异性信息-s yes/reverse/no。htseq默认只保留唯一比对或只有对应一个基因的reads，判断可能属于多个基因的reads将被舍弃；

步骤四，差异分析步骤：

使用DESeq软件对表达后的数据进行差异分析，具体是：

将表达矩阵导入到DESeq包中进行差异分析，保留|log2FC|>1,pvalue<0.05的作为差异表达的转录本；

步骤五，富集分析步骤：

使用topGO软件包进行GO富集分析，具体是：

根据KEGG注释信息，使用phyper进行超几何检验，计算KEGG富集显著性。获得的显著性P值，使用p.adjust进行多重校正，得到校正后的P值，通常选择P<0.05的结果为最后的富集结果；

步骤六，外显子差异分析步骤：

用DEXSeq软件包进行外显子差异分析；

步骤七，可变剪切分析步骤：

使用rMATS进行差异可变剪切分析，获取差异可变剪切数据，具体是：

通过统计模型对不同样本进行可变剪切事件的定量，然后计算P value来表示两组样品在IncLevel(Inclusion Level)水平上的差异；

然后对p value进行校正得FDR值。rMATS可识别的可变剪切事件有5种，分别是skippedexon(SE)外显子跳跃，alternative 5’splice site(A5SS)第一个外显子可变剪切，alternative 3’splice site(A3SS)最后一个外显子可变剪切，mutually exclusiveexons(MXE)外显子选择性跳跃和retained intron(RI)内含子滞留；

步骤八，序列差异分析步骤：

使用varscan软件分析SNP信息，得到样品与参考基因的序列差异；

步骤九，转录本序列制备步骤：

使用stringtie进行转录拼接，重构样品的转录本序列，具体为：

对每个样品单独使用stringtie拼接，使用stringtie merge将所有样品的转录本进行合并。合并后的转录本为最后重构出来的转录本序列。

在本发明的一个优选实施例中，所述根据RNASeq的测序特点的质量控制标准为：

所述测序数据的碱基分布的四条线为：AT平行且接近或GC平行且接近或GC整体分布应该近似于正态分布或不能出现多峰；

测序数据的Duplicate水平与建库的PCR循环数一致。

2.如权利要求1所述的一种适用于oligodT富集的mRNA二代测序结果分析方法，其特征在于，所述构建索引为在hisat2比对时加上链特异性参数，所述链特异性参数当中，当采用dUTP建库方法时，对应的参数为--rna-strandness FR，其他方法当中参数采用默认值。

3.如权利要求1所述的一种适用于oligodT富集的mRNA二代测序结果分析方法，其特征在于，所述步骤二当中的比对为：

模式生物的比对率为不低于80％，多重比对率为不高于10％；如果比对率低，通常需要检查测序样品是否有污染；如果多重比对率高，通常需要检查测序数据的核糖体比例。

4.如权利要求1所述的一种适用于oligodT富集的mRNA二代测序结果分析方法，其特征在于，所述可变剪切分析当中：默认保留差异大于0.01％的可变剪切。

5.如权利要求1所述的一种适用于oligodT富集的mRNA二代测序结果分析方法，其特征在于，所述序列差异分析具体是：

首先需要使用samtools mpileup提取样品在基因组上每一个碱基位置的比对情况；

然后varscan读取每个位置的信息，从中识别其中的错配和插入缺失，并将结果转换成VCF格式输出；

每个SNP需要满足总测序深度大于8；杂合位点的SNP频率以40％-60进行筛选。

Images