• CLC  RNA-seq 数据分析案例

CLC RNA-seq 数据分析案例

演示数据本教程将会使用三个不同品系雌雄黑腹果蝇的RNA-seq数据,这些果蝇是在一年中的不同时间点收集的,并使用了3种不同的RNA-seq方法。目的是,控制其它因素,找到不同品系间差异表达的基因。本次演示仅使用原726个样本中的10个样本,并且每个样本的测序Reads数量降低到仅mapping于染色体2R的250,000条。Excel表中metadata备注信息:-SSR_ID为每个样本的SRA编

演示数据

本教程将会使用三个不同品系雌雄黑腹果蝇的RNA-seq数据,这些果蝇是在一年中的不同时间点收集的,并使用了3种不同的RNA-seq方法。目的是,控制其它因素,找到不同品系间差异表达的基因。

本次演示仅使用原726个样本中的10个样本,并且每个样本的测序Reads数量降低到仅mapping于染色体2R250,000条。

Excel表中metadata备注信息:-SSR_ID为每个样本的SRA编号,是唯一的;- DGRP_Number为果蝇的品系;-Sex为果蝇的雌雄(M F);-Environment为果蝇的收集时间,以23表示不同的时间;-RNA_Prep_MethodRNA制备方法,都使用QIAGEN RNeasy kit,但后续使用离心或者真空方法两种不同方法;-Lane为使用的测序仪。

演示数据下载链接:http://download.clcbio.com/testdata/RNA_Seq_Droso2.zip,下载后解压。

1. 数据导入

1.1  File | Import |Standard Import,选择文件“RNA-Seq statistics workflow”和“gene_association.fb”。

1.2  File | Import | Illumina,选择10个测序文件,均为*.fastq格式文件,不勾选“Paired reads”,如图。点击Next,选择保存位置。

1499411963495037.jpg

1.3 File | Import | Import Metadata,选择文件“DrosophilaMetadata.xls”为样本信息文件,选择10个测序文件为要被关联的文件,匹配模式Exact,完成,保存。

1499411991118376.jpg

1499412048559843.jpg

1.4  打开刚刚保存的metadata table,选择所有的行,单击下方的“Find Associated Data”,将会显示所有和这些行相关的相应测序Reads文件。

1499412077373343.jpg

2. 下载参考基因组

分析之前还需要下载Drosophila melanogaster的参考基因组序列。

单击菜单栏的Download,选择“Download Reference Genome Data”,在对话框中选择“Animals (others)”,下拉选项中选择“Drosophila melanogaster”,下一步中,选择“Download genome sequence”,点击Next,在“Select annotations”窗口中,勾选“Genome Annotations”,本次分析不需要变异信息。选择保存位置,完成。

1499412147389523.jpg

1499412166213609.jpg

最终所有文件如图。

1499412195833783.jpg

3. RNA-seq 计算表达值

将测序Reads map到参考序列上,并计算基因/转录本的表达值。

Toolbox | Transcriptomics Analysis | RNA-Seq Analysis | RNA-Seq Analysis,先要勾选左下角的Batch选项,然后选择10个测序文件,下一步,检查一下选择的文件是否有误,点击Next,勾选“"Genome annotated with genes and transcripts”,定义好已下载的Gene track and mRNA track文件,“Also map to inter-genic regions”选项建议用于真核生物。下一步的mapping参数使用默认值,点击Next。在“Expression level options”对话框,勾选“Use EM estimation”,其他参数默认。选择保存位置,完成。

1499412219184692.jpg1499412245126358.jpg

任务开始以后,可以在左下角看到任务进度,在标准笔记本电脑上,这10个样本的mapping大约需要20min。当任务完成后,在Metadata Elements table中单击下方的Refresh,可以看到分析完成以后,有更多的结果文件关联到样本上。

1499412267516381.jpg

4.  PCA plot

下面使用PCA plot来检查离散样本,并查看metadata中的哪个因素是(导致基因差异表达检出的)品系因素的混淆因素。

Metadata Elements table中,单击“Role”列表头来对所有文件进行排序,选中role为“Gene expression”的十个文件;Toolbox | Transcriptomics Analysis | RNA-Seq Analysis | PCA for RNA-Seq,所有的Gene Expression files已经被预先选中了,保存,完成。

1499412292463495.jpg

打开PCA plot图,样本是依据它们的DGRP_Number号显示的,在右边的setting panel中,对Metadata | Symbol color顺序选择不同的值,可以发现,样本的主要影响因素(PC1)源于它们的性别,第二主要影响因素(PC2)为实验中的所使用的RNA-seq准备方法。由RNA-seq准备方法导致的差异是一种批次效应(Batch Effects),Batch Effects可以由多种原因产生,如使用了同一测序仪的不同lanes,或者样本为不同时间准备的。所以,记录更多的样本信息,可以去除统计学分析中可能的Batch effects,有助于解读我们的数据。

1499412329119886.jpg很明显,果蝇的性别影响很显著,所以有必要在做后续的统计学分析之前,去除雌性果蝇样本。对于RNA 提取方法对本实验的影响,在后续的分析会展示如何去除改影响。

另外你也可以选择3D PCA plot显示方式来调研更多的信息,如图。

1499412350534497.jpg

5. 使用Workflow实现其它统计学分析

双打打开导入的Workflow

1499412370121474.jpg

可以直接从文件导航栏运行该Workflow,但是,为了使用metadata table功能中的文件选择功能,建议将workflow安装到workbench中。

安装方法:右下角单击“installation”,在弹出的窗口中点击Next,选择“Install the workflow on your local computer”,单击完成。

在运行Workflow之前,首先要去除两个雌性果蝇样本。在metadata table中,单击右上角的“Filter”,选择Sex = M,单击Filter,选中保留下来的8个样本,单击下方的“find associated data”,现在只有雄性果蝇的数据保留了下来,依据“Role”排序,仅选择role为“Gene expression”的文件,右键,选择Toolbox | Workflows | RNA-Seq statistics workflow,如图。

1499412397807663.jpg

在弹出的对话框中,点击运行。下面是参数设定:要分析的数据文件已经被默认选中;在"Select GO annotation" 对话框中,选择文件“gene_association.fb”;在“Differential Expression for RNA-Seq”窗口中,参数设定如图,我们要找出那些因品系导致的差异表达基因,而控制那些由RNA提取方法导致的假阳性差异基因,另外,要对不同的品系进行两两比较,所以选择“All group pairs”,单击Next。保存,完成。

1499412420235082.jpg

6 结果解读和解释

Workflow运行完成以后,得到以下文件。

1499412449565394.jpg

6.1 韦恩图 Venn diagram

韦恩图可以展示检测到有多少个转录本在每组比较中是差异表达的。

选中任意两个样本的交叉处,这些基因在这两个样本中是差异表达的,在这个状态下,点击左下角第二个图标,切换到table视图,可以看到那些差异表达的基因仍然是被选中的状态,然后单击下方的“"Create from selection”将其选取出来,这样会打开一个新的窗口仅含有那些被选中的基因。

1499412470864250.jpg

6.2 统计学比较 Statistical comparison

Track List中查看那些统计学意义上发生显著变化的基因。

首先要创建Track ListFile | New | Track list,勾选统计学比较结果和参考序列,点击完成。

1499412499358196.jpg

双击打开其中一个统计学比较结果(statistical comparison track),再重新打开韦恩图,选择两个样本的交叉区域,然后切换到table视图中,点击“Select Genes in Other views”,由于diagramtable都是关联的,韦恩图中选中的那些基因在table中也是自动被选中的,单击“Select genes in other view”,这样那些被选中的基因,同样会在统计学比较结果中选中(statistical comparison track)。而此时,这些被选中的基因在Track List中也会以红线选中。在统计学比较结果(statistical comparison track)中,点击左下角的第二个图标,可以切换到火山图视图,如图。

1499412532117429.jpg

此时,那些在韦恩图中被选中的差异基因,也会在火山图中,以红色着重显示。同样,你也可以在火山图中选中某些基因, Track List中也会跳转到这些基因,并以红线着重显示。

6.3 表达值查看 Expression Browser

Expression Browser可以查看基因和转录本的表达水平和多个样本间的统计学结果,此外,还可以查看每个基因或转录本的GO生物学过程注释信息。

在这里,在韦恩图中那些选中的差异表达的基因,同样会在Expression Browser 中自动选中。

打开Expression Browser,并再次打开韦恩图,选中某个交叉区域,再切换到table视图,单击“Select Genes in Other views”,然后打开Expression Browser,那些差异表达的基因也同时”在Expression Browser中自动选中了,单击下方的“Copy Gene Names to Clipboard”,在右上方打开Filter选项,选择Name,逻辑为 is in list,把刚才复制的基因名子粘贴到第三项中,单击Filter。此时,Expression Browser中仅显示那些之前选中的那些基因。如图。

1499412582296876.jpg你可以单击蓝色的生物学过程名字,会打开一个相应的Gene Ontology网页,同样如果有这个基因的MID RefSeq HGNC 或者 UniProt 等编号,也可单击打开相应介绍。

6.4 热图 Heat map

最后打开RNA-Seq的热图结果,聚类分析的参数已经预设为基于变异系数(coefficient of variation/ relative standard deviation)选出25个最有意义的基因,并且这些样本同样在水平轴上进行聚类(in unsupervised clustering)。

在右边的panel中,轮流选择metadata layers,你会发现这25个最有意义的基因是由于RNA提取方法导致的差异表达,这也是我们能预料到的,因为在PCA plot中我们已经发现RNA提取方法要比品系差异的影响大。

那么,如何在热图中用那些由菌株品系导致的差异表达基因对样本进行聚类分析呢?

回到“samplesmetadata table中,确保这8个雄性样本的“Gene expression”文件为选中状态;右键,运行“Create Heat Mapfor RNA-Seq”工具,在打开的对话框中,那些样本已经被默认选中,点击Next,参数默认,点击Next,在“Set filtering”窗口中,选择“Filter by statistics”,在“Statistical comparison”中选择其中一个样本的两两比较结果,点击OK,其他参数默认,完成。

1499412610260686.jpg1499412631579767.jpg1499412648125721.jpg

打开新出现的热图结果,很明显这些样本是依据品系进行了聚类。

1499412671763814.jpg

以上是一个比较完整的转录组测序数据分析案例,在这里,我们不仅可以检测出Batch effect对实验的影响,在后续的分析中,我们还可以控制这些因素的影响,从而找出品系特异的差异基因。

关注微信