设为首页收藏本站
开启辅助访问
切换到窄版

 找回密码
 注册

QQ登录

只需一步,快速开始

搜索
查看: 679|回复: 0

生物云计算数据中心

[复制链接]
发表于 2016-7-18 16:53:52 | 显示全部楼层 |阅读模式
联系人:钟老师 QQ:946326605 邮箱:946326605@qq.com

新一代测序自从大红大紫以来  数据分析一直是其痛点,几十到一百多BP的短读长想要通过运算拼接出几十个G的基因组序,所需的数据分析运算堪称海量、所需的运算能力和存储的设备,以及信息学领域人才的缺乏曾一度令很多生命学实验室望而却步。不过,还没修炼问生物信息学大人的生物学精英们也不能就此停下脚步,比较新技术发掘的海量数据为我们带来了浅说未有的机遇,怎么忍心错过这个好机会。如果能把这些麻烦打包出去,不就简单多了吗?对海量运算能力和存储能力的需求,催生了可以领域的生物云计算平台。


生物云计算数据中心是服务于生物产业的、基于云计算架构搭建的数据中心。由于基因组测序数据迅猛增加,使得生物产业对于计算与存储的需求呈现指数级的增长速度。这种由于不断产生的信息洪流而形成的对计算与存储能力超常规的增长,使得生物云计算数据中心区别于其他行业的数据中心,具有独特的特征和对技术的更高要求。本次生物云平台建设,主要目标是搭建生物云基础设施与基础管理平台,在按规划纳入部分现有IT 系统之外,建立起统一服务管理模式及新服务节点建设标准。北京市计算中心生物计算事业部的-生物云计算平台。
  基因盒子生物云平台
  基因盒子生物云平台(www.dnabox.cn)是由北京市计算中心、北京唐唐天下生物医学信息科技有限公司合作开发的面向生命科学及相关应用领域的云计算平台,旨在为各科研院所及检测机构提供高通量生物数据的存储及计算分析支撑服务,平台目前部署在北京市计算中心与中金数据系统计算集群之上,是国内首个跨多个公有云的生物云计算平台,目前,平台已有软件及分析流程200余种,可控计算节点2万余个,存储空间10PB,可以满足客户的任何软件及硬件需求,自2015年发布以来,注册用户单位数百家,已经成功完成数万个计算任务,通过平台分析计算,辅助用户在scientific report等期刊发表论文多篇。



平台提供的数据分析报告如下:



测序实验流程
依据Illumina TruSeq™ RNA 样品准备试剂盒进行实验步骤描述,详见文档RNA.pdf.


1. 样品RNA质量检测


此步由委托的实验部门完成:
1) 1%的琼脂糖电泳检测RNA样品是否有降解以及杂质
2) NanoPhotometer® 分光光度计检测样品纯度(IMPLEN, CA, USA)
3) Qubit® 2.0 Flurometer (Life Technologies, CA, USA)检测RNA样品浓度
4) 安捷伦2100 RNA Nano 6000 Assay Kit (Agilent Technologies, CA, USA)检测RNA样品的完整性


2. 样品RNA 转录组文库制备


每个样品取总量3ug RNA上样。所有4个样本的RIN(RNA Integrity Number)值超过8.0。依照使用手册,用TruSeq RNA试剂盒制备测序文库,并添加索引码标记4个样本各自的序列。简单的说,通过带polyT的磁珠,从总的RNA中纯化出mRNA。在Illumina专门的碎片缓冲液中,升高温度,用二价阳离子进行mRNA片段化。用随机的寡核苷酸和SuperScript II 合成cdna的第一条链。随后,用DNA聚合酶I和RNase H 合成cDNA的第二条链。剩下的突出端核酸用核酸外切酶/聚合酶转化成平端,并移除相关酶。DNA片段3’端腺苷酰化后,连接上PE接头寡核苷酸后准备杂交。用AMPure XP系统纯化文库片段,筛出200bp的cDNA片段。挑出两端都有接头的DNA片段用PCR Primer Cocktail的10个PCR循环去扩增。用高敏感性的DNA矩阵对PCR产物进行纯化和定量。



图1 原核生物转录组建库流程

3. 库检


文库构建完成后,先使用 Qubit2.0 进行初步定量,稀释文库至 1ng/ul,随后使用 Agilent 2100 对文库的 insert size 进行检测, insert size 符合预期后, 使用 Bio-RAD CFX 96 荧光定量 PCR 仪, Bio-RAD KIT iQ SYBR GRN 进行 Q-PCR,对文库的有效浓度进行准确定量(文库有效浓度> 2nM),以保证文库质量。
4. 上机测序


使用 TruSeq PE Cluster Kit v3-cBot-HS( Illumia) 试剂在 cBot 上生成簇。之后在 HiSeq 2500 测序平台上运行双端测序程序( PE100), 得到 100bp或 125bp 的双端测序 reads。
数据分析流程
针对HiSeq2500测序数据进行转录组分析,通过参考基因组比对来发现物种转录本的表达情况,发现特异的表达序列,进行功能基因的分析,转录组分析如下。实验设计允许的条件下可以进行差异表达分析。


图2 有参转录组的分析流程图

测序数据质量情况汇总
我们对原始数据按照如下步骤进行数据的整理和预处理,利用trimmomatic (v0.32)[1]进行数据质量过滤,保留质量在20以上的碱基;将过滤后的数据按照长度进行筛选,去掉长度小于25,或者只有一端的read。
下表为测序数据统计:


表2 原始测序数据质量情况汇总

Sample        Raw Reads        Bases        GC (%)        Q20        Q30        Avg.Quality
chanhou17d        18540741        0.297(GB)        51        78.76%        74.86%        93.43
chanhou17d        18540741        0.319(GB)        51        78.61%        74.74%        93.28
chanqian17d        20287141        0.373(GB)        53        80.01%        76.16%        95.32
chanqian17d        20287141        0.401(GB)        53        80.06%        76.22%        95.58
chanhou1d        18871382        0.328(GB)        53        79.84%        76.01%        94.96
chanhou1d        18871382        0.349(GB)        54        79.43%        75.60%        94.36
表注:
(1) clean reads: 指按条件(q大于20,长度大于25)筛选后的reads数目。
(1) Percentage: 指按条件(q大于20,长度大于25)筛选后的reads数目占原始数据的比例。
全基因组比对
预处理后的数据,使用HISAT2)[2]比对到参考基因组上。参考基因组和基因组信息下载自Ensembl数据库。
Tophat算法示意图如下:


图3 HISAT2的算法示意图

1 READS与参考基因组比对情况统计


表3 数据的Mapping结果统计

sample        Total reads after filtered        Mapped on reference        Unmapped        Multi-mapped
chanhou17d        18540741        10848164 (58.51%)        4935978 (26.62%)        2756599 (14.87%)
chanqian17d        20287141        12020963 (59.25%)        5777144 (28.48%)        2489034 (12.27%)
chanhou1d        18871382        12600884 (66.77%)        3757321 (19.91%)        2513177 (13.32%)
表注:
(1) Total reads after filtered:测序序列经过测序数据过滤后的数量统计(Clean data)。
(2) Mapped on reference:能定位到基因组上的测序序列的数量的统计;一般情况下,如果不存在污染并且参考基因组选择合适的情况下,这部分数据的百分比大于 70%。
(3) Mapped in genes:测序序列比对到基因区的reads和比例。
(4) Mapped in exons:测序序列比对到外显子区的reads和比例。


2 READS在参考基因组不同区域的分布情况


将比对到基因组上的reads分布情况进行统计,定位区域分为Exon(外显子)、Intron(内含子)和Intergenic(基因间区)。
在基因组注释较为完全的物种中,比对到Exon(外显子)的reads含量最高,比对到Intron(内含 子)区域的reads来源于pre-mRNA的残留及可变剪切过程中发生的内含子滞留事件导致的,而比对到Intergenic(基因间区)的reads 是因为基因组注释不完全。
Reads在基因组不同区域的分布如下图:



123
图4 Read在基因组不同区域的分布
基因表达水平分析
一个基因表达水平的直接体现就是其转录本的丰度情况,转录本丰度越高,则基因表达水平越高。在RNA-seq分析中,我们可以通过定位到基因组区域或基因外显子区的测序序列(reads)的计数来估计基因的表达水平。
TOPHAT比对到参考基因组上的结果,交给cufflinks(v2.0.2)[4]进行处理,以估计每个转录本的表达值。所有基因的位置,类型和表达量等信息请查看expression_all.xls。参考基因组和基因组信息下载自Ensembl数据库,版本号:Arabidopsis_thaliana.TAIR10。
表5 表达量FPKM的区间分布统计

Sample        0 ~ <1        1 ~ <10        10 ~ <50        50 ~ <150        150 ~ <500        500 ~ <1500        1500 ~ <5000        5000 ~



通过所有基因的RPKM分布图以及盒形图对不同实验条件下的基因表达水平进行比较。对于同一实验条件下的重复样品,最终的RPKM为所有重复数据的平均值。
结果如下:
样品基因的FPKM分布如下图所示:




图5 样品基因的FPKM分布
图注:
(1) 图一:RPKM盒形图,横坐标为样品名称,纵坐标为log10(RPKM+1),每个区域的盒形图对五个统计量(至上而下分别为最大值,上四分位数,中值,下四分位数和最小值) 。
(2) 图二:RPKM分布图,横坐标为log10(RPKM+1), 纵坐标为基因的密度。
对每组的表达基因进行比较得到如下的Venn图:


图6 样品的基因表达Venn图
图注:Venn图表示基因在样品之间的表达情况,对于设置了重复实验的样品,只要基因在一个重复中表达,即认为该基因在该样品中表达。

4基因差异表达分析
1 差异表达分析


基因差异表达的输入数据为基因表达水平分析中得到的readcount数据。对于有生物学重复的样品,我们采用edgeR(v2.2.5)[6]进行分析,采用Quantile的Normalization方法,选取Pvalue小于0.05并且FDR小于0.1的基因作为差异表达的基因列出。
图注:详细信息见附件中的differential目录下excel文件,表头说明如下:
有显著性差异表达的基因用红色点表示;横坐标代表基因在不同样本中表达倍数变化; 纵坐标代表基因表达量变化差异的统计学显著性。
(1) logFC变化倍数以2为底数的log值,即表达量比值的对数;
(2) logCPM指log-Concentration, 一种表示表达水平的方法,CPM是counts per million的缩写;
(3) Pvalue统计检验的显著值,是用负二项分布模拟后检验后的统计值;
(4) FDR 对Pvalue进行多重检验后得到的显著值;
(5) Trend 表示基因上调还是下调;
(6) 每个group对应的那一列就是FPKM值。


差异基因列表如下:
表6 差异基因列表


结果如下:
chanhou1d-v-chanqian17d
chanhou1d-v-chanqian17d
chanhou1d-v-chanhou17d
chanhou1d-v-chanhou17d
chanhou17d-v-chanqian17d
chanhou17d-v-chanqian17d
chanhou1d-v-chanqian17d
chanhou1d-v-chanqian17d
chanhou1d-v-chanhou17d
chanhou1d-v-chanhou17d
chanhou17d-v-chanqian17d
chanhou17d-v-chanqian17d
GID        Col_0        Msp40_1        logFC        logCPM<150        PValue<500        FDR<1500        trend<5000
ATCG00390        0        55060.4        18.80452        14.63804        3.89E-14        9.19E-10        UP
ATCG00390        0        16885.1        17.09926        12.93299        2.00E-12        2.36E-08        UP


表注:
(1) GID为差异基因ID;
(2) 第二列表示差异组样品RPKM值;
(3) 第三列表示差异组样品RPKM值;
(4) logFC变化倍数以2为底数的log值,即表达量比值的对数;
(5) logCPM指log-Concentration, 一种表示表达水平的方法,CPM是counts per million的缩写;
(6) Pvalue统计检验的显著值,是用负二项分布模拟后检验后的统计值;
(7) FDR 对Pvalue进行多重检验后得到的显著值;
(8) Trend 表示基因上调还是下调;



2 表达聚类分析


聚类分析用于判断差异基因在不同实验条件下的表达模式;通过将表达模式相同或相近的基因聚集成类,从而识别未 知基因的功能或已知基因的未知功能;因为这些同类的基因可能具有相似的功能,或是共同参与同一代谢过程或细胞通路。以不同实验条件下的差异基因的RPKM 值为表达水平,做层次聚类(hierarchical clustering)分析,不同颜色的区域代表不同的聚类分组信息,同组内的基因表达模式相近,可能具有相似的功能或参与相同的生物学过程。
所有样品的表达量整理之后,交给Cluster3(v3.0)[7]进行基因和样品的双向聚类。详细结果见cluster目录。
聚类的目的是基于物体的相似性将物体分为不同的组,聚类图相当于某种意义上的分类图,在此分析中,我们以表达量差异进行比较,作为衡量关系远近的标准.




图8 表达聚类图
注:聚类图非常大,建议使用TreeView软件对结果文件进行查看。
6 差异基因GO富集分析
Gene Ontology(简称 GO, http://www.geneontology.org/) 是基因功能国际标准分类体系。根据实验目的筛选差异基因后,研究差异基因在 Gene Ontology 中的分布状况以期阐明实验中样本差异在基因功能上的体现。GO注释信息下载自Ensembl的Biomart数据库。GO的富集分析使用BiNGO(v2.44)[8]进行分析,GO富集结果见GOenrichment 文件夹,并得到下面的GO_biological progress和GO_cellular component的富集结果,橘色和黄色分别表示在该词条上有富集。


1 差异基因GO富集列表


结果如下所示:
表7 差异基因GO富集列表

Class        GO.ID        Term        Annotated        Significant        Expected        classicFisher        Genelist
BP        GO:0015979        photosynthesis        401        4        0.07        9.30E-08        ATCG00080,ATCG00550,ATCG00570,ATCG00690
BP        GO:0009767        photosynthetic electron transport chain        69        1        0.01        0.012        ATCG00570


表注:完整内容见附件表格
(1) Class:该GO的类别(CC:细胞组分;BP:生物学过程;MF:分子功能)
(2) GO.ID:Gene Ontology数据库中唯一的标号信息
(3) Term:Gene Ontology功能的描述信息
(4) Annotated:该GO上的所有基因数目
(5) Sigificant:与该GO相关的差异基因数目
(6) Expected:富集分析统计学显著水平
(7) classic Fsher:Fisher检验显著水平
(8) Genelist:与该GO相关的差异基因列表


2 差异基因GO富集DAG图


有向无环图(Directed Acyclic Graph,DAG)为差异基因GO富集分析结果的图形化展示方式,分支代表包含关系,从上至下所定义的功能范围越来越小,一般选取GO富集分析的结果前 10位作为有向无环图的主节点,并通过包含关系,将相关联的GO Term一起展示,颜色的深浅代表富集程度。我们的项目中分别绘制生物过程(biological process)、分子功能(molecular function)和细胞组分(cellular component)的DAG图。
如图所示(清晰的矢量图请查看pdf目录下对应名称的文件):





123456789101112


图9 GO富集有向无环图
图注:详细结果见GO目录下的结果,其中图片是差异基因的GO分类图,excel文件是对应GO分类图的差异基因具体的数目。
每个节点代表一个GO术语,方框代表的是富集程度为TOP10的GO,颜色的深浅代表富集程度,颜色越深就表示富集程度越高,每个节点上展示了该TERM的名称及富集分析的p-value
7 PATHWAY分析
在生物体内,不同基因相互协调行使其生物学功能,通过Pathway显著性富集能确定差异表达基因参与的最主 要生化代谢途径和信号转导途径。KEGG(Kyoto Encyclopedia of Genes and Genomes)是有关Pathway的主要公共数据库[9]。Pathway显著性富集分析以KEGG Pathway为单位,应用超几何检验,找出与整个基因组背景相比,在差异表达基因中显著性富集的Pathway。我们使用KOBAS(2.0)[10]进行Pathway富集分析。
结果:详细结果见附件Pathway目录下的内容,其中子目录是每组差异基因的pathway图,XLS文件是基因的Pathway注释。


图10 代谢通路示意图(具体数据请查看Pathway目录中的结果)
9 可变剪切和新转录本分析
方法:先用cufflinks处理HISAT2的比对结果,对每个样品进行转录本的组装,然后合并所有样品的转录本结果,然后使用ASprofile[8-10]对合并后的结果进行可变剪切事件的检测。
结果:见附件中的cuffcompare/AS_events.gtf。


转录本的剪接形式统计

Sample name        XSKIP        XMIR        SKIP        IR        XIR        AE        MIR        TSS        TTS        XMSKIP        XAE        MSKIP
chanhou17d        64        1        627        94        16        596        1        11656        11195        5        341        68
chanqian17d        129        2        1351        262        43        1238        10        17502        16584        9        233        143
chanhou1d        83        0        959        157        32        976        8        15757        14938        8        389        86


表中可变剪切编码的含义如下:
AS name        AS code
Approximate SKIP (XSKIP_ON,XSKIP_OFF pair)        XSKIP
Approximate MIR (XMIR_ON, XMIR_OFF pair)        XMIR
Skipped exon (SKIP_ON,SKIP_OFF pair)        SKIP
Intron retention (IR_ON, IR_OFF pair)        IR
Approximate IR (XIR_ON, XIR_OFF pair)        XIR
Alternative exon ends (5', 3', or both)        AE
Multi-IR (MIR_ON, MIR_OFF pair)        MIR
Alternative 5' first exon (transcription start site)        TSS
Alternative 3' last exon (transcription terminal site)        TTS
Approximate MSKIP (XMSKIP_ON,XMSKIP_OFF pair)        XMSKIP
Approximate AE        XAE
Multi-exon SKIP (MSKIP_ON,MSKIP_OFF pair)        MSKIP
合并后的转录本结果,通过cuffcompare与已有的基因结构进行比较,得到的未知转录本,将被认为是新的转录本;

您需要登录后才可以回帖 登录 | 注册

本版积分规则

QQ|申请友链|小黑屋|手机版|Archiver|生物信息学论坛 ( 蜀ICP备09031721号  

GMT+8, 2017-1-17 10:52 , Processed in 0.105806 second(s), 24 queries .

Powered by Discuz! X3

© 2001-2013 Comsenz Inc.

快速回复 返回顶部 返回列表