转录组学 01/15, 2017 ventson@zju.edu.cn
前言 基因组学 转录组学 from en.wikipedia 中心法则 : 遗传信息传递 蛋白质组学 2
技术发展 Experiment-based Northern blot RT-PCR Hybridization-based Microarray Sequencing-based SAGE CAGE MPSS Advanced seq NGS 3GS Single cell 转录组学研究技术革新 3
应用 差异表达 可变剪切 共表达 转录调控 4
RNA 测序 (RNA-sequencing) from GATC Biotech 1. 试验设计 2. 测序流程 3. 数据分析 4. 验证实验 5
试验设计 问题导向型 数据导向型 生物学重复 (3-5 个 ) 样本提取 ( 分类和保存 ) 测序深度 ( 简单基因表达分析需 5M 以上 reads, 小 RNA 至少 30M) 数据异质性 ( 平台 个体差异 ) 确定分析流程 分析工具选用 文库构建 ( 链特异性非特异性 ) 测序策略 ( 单端和双末端 ) 测序平台 ( 读长 通量和准确率等 ) 6
测序流程 mrna:poly A 富集 ncrna:rrna 移除 Griffith, M. (2015) PLoS computational biology 7
数据分析流程 系统配置 差异表达 聚类分析 数据获取 表达定量 功能富集 质量控制 比对组装 共表达网络 RNA-seq 数据分析常规流程 8
系统配置 9
数据获取 NCBI SRA TCGA/GDC(cancer) fastq-dump (SRAToolkit) Fastq 文件格式 : EBI ArrayExpress 公共数据库 测序公司 10
质量控制 去接头 ; 过滤低质量 reads FastQC 测序质量评估 FASTX-Toolkit,Trimmomatic 质量控制 11
比对 (reads mapping) 非剪接比对 Bowtie,BWA ( 不考虑可变剪切 ) 剪接比对 TopHat,STAR,HISAT/GSNAP,MapSplice(SNP) TopHat 工作原理 Trapnell, C. (2009) Bioinformatics 12
比对结果 比对结果文件 SAM(SAMtools) 比对结果可视化 IGV (local) 比对结果评估 Qualimap (summary) 13
表达定量 Reads counting Normalization 只保留唯一匹配 reads HTSeq-count,featureCounts 保留多重匹配 reads Cufflinks,StringTie,RSEM RPKM,FPKM,TPM 校正测序深度 基因长度 DESeq/edgeR(TMM) 校正异常高表达基因 14
比对组装策略选择 Conesa, A. (2016) Genome Biology 15
差异表达分析 选取样本 : 样本相关性, 大样本降维 ( 主成分分析 ) 模型选择 : 高斯分布 ( 正态 ), 泊松分布 (v=μ), 负二项分布 (v=μ+αμ 2 ) 差异检验 : 组间差异 ( 处理差异 ) 组内差异 ( 个体差异 )?= 0 筛选条件 :p value( 多重检验校正 ) & FoldChange( 差异倍数 ) 工具 版本 标准化方式 模型假设 统计检验 edger 3.18.1 TMM/Upper quartile/rle 负二项分布 Exact test DESeq2 1.16.1 DESeq sizefactors 负二项分布 Wald test/lrt bayseq 2.10.0 quantile/tmm/total 负二项分布 empirical Bayesian NOIseq 2.20.0 RPKM/TMM/Upper quartile 非参数 Condition vs. null Limma 3.32.10 TMM voom 转换 Empirical Bayes Cuffdiff2 2.2.1 Geometric/quartile/FPKM β 负二项分布 t-test EBSeq 1.16.0 DESeq median normalization 负二项分布 empirical Bayesian 常用差异表达分析工具比较 16
聚类分析 单基因分析 vs. 基因模块分析常用聚类方法 :K-means(K 均值 ), 层次聚类,SOM( 自组织映射 ),FCM( 模糊 C 均值 ) 基因表达聚类结果 (pheatmap) 17
富集分析 GO/KEGG 基因集 功能集 常用工具 :DAVID,agriGO, GSEA,IPA,clusterProfiler GSEA 原理 Subramanian, A. (2005) PNAS 表达量 - 样本相关性排序, 功能基因集分布, 计算富集得分 超几何分布 Fisher 精确检验 特定功能集 S 不属于功能集 S 总基因数 目标基因 x k-x k 背景基因 M N-M N 18
共表达网络 不同样本表达模式相似的基因功能应该也类似 根据表达量计算相关性矩阵, 构建共表达网络 基因集 相互作用 基因网络 无标度网络 Interaction, 相关性系数 : Pearson,Spearman WGCNA 权重基因共表达网络分析 核心基因 (Hub genes) MCODE 网络模块挖掘 ( 子网络 ) CytoScape 网络可视化 19
验证试验 相关性 PCR, 凝胶电泳 基因敲除, 敲减, 过表达 因果关系 20
拓展 全长转录本 ( 三代测序 ) 单细胞测序 (single cell) 技术革新 多组学整合 : 基因组, 表观组, 蛋白组, 代谢组, 表型组 非编码 :lncrna,circrna 表观转录组 :m 6 A 修饰 整合应用 21
操作实践 核心分析 转录组核心数据分析在 linux 系统中完成 ( 内存 时间 ), 具体步骤可参考 生物信息学 第 3 版书籍及视频教程 表达量矩阵 : 基因 X 样本 22
操作实践 统计分析及可视化 1. 软件安装 2. 差异分析 3. 聚类分析 4. 功能富集 23
软件安装 http://www.bioconductor.org/ DESeq2( 差异表达分析 ) ggplot2( 作图 ) pheatmap( 聚类可视化 ) CRAN # 从 cran 安装 pheatmap,ggplot2 choosecranmirror() install.packages( pheatmap ) library(pheatmap) install.packages( ggplot2 ) library(ggplot2) # 从 bioconductor 安装 DESeq2 choosebiocmirror() #China Anhui source("http://www.bioconductor.org/bioclite.r") bioclite("deseq2") library(deseq2) 24
差异分析 # 数据预处理 setwd( F:/rnaseq/data ) # 设置工作目录 ( 根据自己存放数据的目录修改 ) library(deseq2) # 载入 DESeq2 包 #reads 计数数据表操作 counttable <- read.table("count.txt", sep = "\t", header = FALSE) # 读入 reads 计数矩阵 tail(counttable) # 查看数据表最后的 小尾巴 counttable <- counttable[- c(33611:33615), ] # 去除描述行 rownames(counttable) <- counttable$v1 # 将基因 ID 设置为行名 counttable <- counttable[, - 1] # 删除基因 ID 列 colnames(counttable) <- c("srr3418005", "SRR3418006", "SRR3418019", "SRR3418020") # 更改数据表列名 counttable <- counttable[- which(rowsums(counttable) < 4), ] # 过滤 count 总数小于 4 的基因 nrow(counttable) # 查看数据表行数 ( 基因个数 ) tail(countable) # 查看修改后的数据表末尾六行 # 设置样本处理信息 ( 实验 vs. 对照 ) coldata <- data.frame(row.names = colnames(counttable), condition = c("aba", "mock", "ABA", "mock")) 25
差异分析 #DESeq2 操作 # 生成 DESeqDataSet 数据集 dds <- DESeqDataSetFromMatrix(countData = counttable, coldata = coldata, design = ~ condition) dds # 查看数据集 dds$condition # 查看样本处理信息 dds$condition <- relevel(dds$condition, "mock") # 更改 mock 水平 ( 使 DESeq 计算 FoldChange 时 mock 组作为分母 ) dds$condition # 查看更改水平后的样本处理信息 dds <- DESeq(dds) # 差异表达计算 res <- results(dds) # 生成差异表达结果 summary(res) # 查看总结信息 ( 表达上调, 下调等 ) resordered <- res[order(res$padj), ] # 按照校准后 p 值排序 write.csv(resordered, "DESeq2_results_all.csv") # 将差异表达分析结果输出到 csv 文件 deg <- subset(resordered, padj <= 0.01 & abs(log2foldchange) >= 2) # 筛选显著差异表达基因 (padj 小于 0.01 且 FoldChange 绝对值大于 4) summary(deg) # 查看筛选后的总结信息 write.csv(deg, "DESeq2_results_significant.csv") # 将差异表达显著的结果输出到 csv 文件 26
差异分析 #volcano plot 火山图 setwd("f:/rnaseq/data") library(ggplot2) volcano_data <- read.csv("deseq2_results_all.csv", row.names = "X") # 读入差异表达结果 volcano_data <- na.omit(volcano_data) # 删除含 NA 的行 significant <- as.factor(abs(volcano_data$log2foldchange) >=2 & volcano_data$padj <= 0.01) # 设置显著性阈值 ggplot(volcano_data, aes(x = log2foldchange, y = - log10(padj))) + geom_point(aes(shape = significant, color = significant)) + xlim(c(-10, 10)) + labs(x = log2foldchange, y = -log10 padj ) + scale_y_continuous(limits = c(0, 20), expand = c(0, 0)) + scale_shape_discrete(labels =c ( no, yes )) + scale_color_discrete(labels = c( no, yes )) #ggplot2 命令 27
聚类分析 #heatmap 聚类热图 setwd("f:/rnaseq/data") library(pheatmap) deseq_results_significant <- read.csv("deseq2_results_significant.csv", row.names = "X") # 读入显著差异表达结果 significant_genes <- rownames(deseq_results_significant) # 提取显著差异基因 fpkm_gtf <- read.table("fpkm.gtf") # 读入 FPKM 注释文件 fpkm_gtf <- fpkm_gtf[-which(fpkm_gtf$v2 == "-" fpkm_gtf$v2 == "."), ] # 删除 Gene ID 未知的行 rownames(fpkm_gtf) <- fpkm_gtf$v2 # 将 Gene ID 设置为行名 fpkm_significant_genes <- fpkm_gtf[significant_genes, 9:12] # 提取显著差异基因的 FPKM 值 colnames(fpkm_significant_genes) <- c("srr3418005", "SRR3418006", "SRR3418019", "SRR3418020") # 设置列名 fpkm_significant_genes <- na.omit(fpkm_significant_genes) # 删除含 NA 值的行 pheatmap(log2(t(fpkm_significant_genes + 1)), show_colnames = FALSE) # 所有差异基因热图 pheatmap(log2(t(fpkm_significant_genes[1:30, ] + 1))) # 差异基因 top30 热图 28
功能富集 DAVID 物种覆盖较全 ; 数据更新慢 agrigo 植物基因富集专用 WEGO 富集结果可视化 clusterprofiler 实时抓取 ; 富集方法全面 ;R 语言 MetaScape 操作简单 ; 可视化效果好 ; 物种较少 29
DAVID 使用 点击开始分析 30
DAVID 使用 第一步 : 导入基因列表 / 文件 第三步 : 基因列表 / 背景 示例文件 第四步 : 提交运行 第二步 : 选择 ID 类型 31
DAVID 使用 选择功能类型 (GO,KEGG) 功能分析 32
DAVID 使用 功能聚类集 富集结果列表 33
MetaScape 使用 上传基因 选择物种 提交运行 34
总结 什么是转录组学? RNA-seq 的研究内容? 如何分析 RNA-seq 数据? 35
转录组分析 谢谢! 浙江大学生命科学学院生物信息学实验室 http://bis.zju.edu.cn 36