PowerPoint 演示文稿

Size: px

Start display at page:

Download "PowerPoint 演示文稿"

当谈
5 years ago
Views:

1 转录组学 01/15, 2017

2 前言基因组学转录组学 from en.wikipedia 中心法则 : 遗传信息传递蛋白质组学 2

3 技术发展 Experiment-based Northern blot RT-PCR Hybridization-based Microarray Sequencing-based SAGE CAGE MPSS Advanced seq NGS 3GS Single cell 转录组学研究技术革新 3

4 应用差异表达可变剪切共表达转录调控 4

5 RNA 测序 (RNA-sequencing) from GATC Biotech 1. 试验设计 2. 测序流程 3. 数据分析 4. 验证实验 5

6 试验设计问题导向型数据导向型生物学重复 (3-5 个 ) 样本提取 ( 分类和保存 ) 测序深度 ( 简单基因表达分析需 5M 以上 reads, 小 RNA 至少 30M) 数据异质性 ( 平台个体差异 ) 确定分析流程分析工具选用文库构建 ( 链特异性非特异性 ) 测序策略 ( 单端和双末端 ) 测序平台 ( 读长通量和准确率等 ) 6

7 测序流程 mrna:poly A 富集 ncrna:rrna 移除 Griffith, M. (2015) PLoS computational biology 7

8 数据分析流程系统配置差异表达聚类分析数据获取表达定量功能富集质量控制比对组装共表达网络 RNA-seq 数据分析常规流程 8

9 系统配置 9

10 数据获取 NCBI SRA TCGA/GDC(cancer) fastq-dump (SRAToolkit) Fastq 文件格式 : EBI ArrayExpress 公共数据库测序公司 10

11 质量控制去接头 ; 过滤低质量 reads FastQC 测序质量评估 FASTX-Toolkit,Trimmomatic 质量控制 11

12 比对 (reads mapping) 非剪接比对 Bowtie,BWA ( 不考虑可变剪切 ) 剪接比对 TopHat,STAR,HISAT/GSNAP,MapSplice(SNP) TopHat 工作原理 Trapnell, C. (2009) Bioinformatics 12

13 比对结果比对结果文件 SAM(SAMtools) 比对结果可视化 IGV (local) 比对结果评估 Qualimap (summary) 13

14 表达定量 Reads counting Normalization 只保留唯一匹配 reads HTSeq-count,featureCounts 保留多重匹配 reads Cufflinks,StringTie,RSEM RPKM,FPKM,TPM 校正测序深度基因长度 DESeq/edgeR(TMM) 校正异常高表达基因 14

15 比对组装策略选择 Conesa, A. (2016) Genome Biology 15

16 差异表达分析选取样本 : 样本相关性, 大样本降维 ( 主成分分析 ) 模型选择 : 高斯分布 ( 正态 ), 泊松分布 (v=μ), 负二项分布 (v=μ+αμ 2 ) 差异检验 : 组间差异 ( 处理差异 ) 组内差异 ( 个体差异 )?= 0 筛选条件 :p value( 多重检验校正 ) & FoldChange( 差异倍数 ) 工具版本标准化方式模型假设统计检验 edger TMM/Upper quartile/rle 负二项分布 Exact test DESeq DESeq sizefactors 负二项分布 Wald test/lrt bayseq quantile/tmm/total 负二项分布 empirical Bayesian NOIseq RPKM/TMM/Upper quartile 非参数 Condition vs. null Limma TMM voom 转换 Empirical Bayes Cuffdiff Geometric/quartile/FPKM β 负二项分布 t-test EBSeq DESeq median normalization 负二项分布 empirical Bayesian 常用差异表达分析工具比较 16

17 聚类分析单基因分析 vs. 基因模块分析常用聚类方法 :K-means(K 均值 ), 层次聚类,SOM( 自组织映射 ),FCM( 模糊 C 均值 ) 基因表达聚类结果 (pheatmap) 17

18 富集分析 GO/KEGG 基因集功能集常用工具 :DAVID,agriGO, GSEA,IPA,clusterProfiler GSEA 原理 Subramanian, A. (2005) PNAS 表达量 - 样本相关性排序, 功能基因集分布, 计算富集得分超几何分布 Fisher 精确检验特定功能集 S 不属于功能集 S 总基因数目标基因 x k-x k 背景基因 M N-M N 18

19 共表达网络不同样本表达模式相似的基因功能应该也类似根据表达量计算相关性矩阵, 构建共表达网络基因集相互作用基因网络无标度网络 Interaction, 相关性系数 : Pearson,Spearman WGCNA 权重基因共表达网络分析核心基因 (Hub genes) MCODE 网络模块挖掘 ( 子网络 ) CytoScape 网络可视化 19

20 验证试验相关性 PCR, 凝胶电泳基因敲除, 敲减, 过表达因果关系 20

21 拓展全长转录本 ( 三代测序 ) 单细胞测序 (single cell) 技术革新多组学整合 : 基因组, 表观组, 蛋白组, 代谢组, 表型组非编码 :lncrna,circrna 表观转录组 :m 6 A 修饰整合应用 21

22 操作实践核心分析转录组核心数据分析在 linux 系统中完成 ( 内存时间 ), 具体步骤可参考生物信息学第 3 版书籍及视频教程表达量矩阵 : 基因 X 样本 22

23 操作实践统计分析及可视化 1. 软件安装 2. 差异分析 3. 聚类分析 4. 功能富集 23

packages( ggplot2 ) library(ggplot2) # 从 bioconductor 安装

24 软件安装 DESeq2( 差异表达分析 ) ggplot2( 作图 ) pheatmap( 聚类可视化 ) CRAN # 从 cran 安装 pheatmap,ggplot2 choosecranmirror() install.packages( pheatmap ) library(pheatmap) install.packages( ggplot2 ) library(ggplot2) # 从 bioconductor 安装 DESeq2 choosebiocmirror() #China Anhui source(" bioclite("deseq2") library(deseq2) 24

$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[-

设置为行名 counttable <- counttable[, - 1] # 删除基因 ID 列 colnames(counttable) <- c("srr3418005", "SRR3418006", "SRR3418019", "SRR3418020") # 更改数据表列名 counttable <- counttable[-

25 差异分析 # 数据预处理 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("srr ", "SRR ", "SRR ", "SRR ") # 更改数据表列名 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 = ~

FoldChange 时 mock 组作为分母 ) dds$condition # 查看更改水平后的样本处理信息 dds <- DESeq(dds) # 差异表达计算 res <- results(dds) # 生成差异表达结果

csv(resordered, "DESeq2_results_all.csv") # 将差异表达分析结果输出到 csv 文件 deg <- subset(resordered, padj <= 0.

26 差异分析 #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

27 差异分析 #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

28 聚类分析 #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("srr ", "SRR ", "SRR ", "SRR ") # 设置列名 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

29 功能富集 DAVID 物种覆盖较全 ; 数据更新慢 agrigo 植物基因富集专用 WEGO 富集结果可视化 clusterprofiler 实时抓取 ; 富集方法全面 ;R 语言 MetaScape 操作简单 ; 可视化效果好 ; 物种较少 29

30 DAVID 使用点击开始分析 30

31 DAVID 使用第一步 : 导入基因列表 / 文件第三步 : 基因列表 / 背景示例文件第四步 : 提交运行第二步 : 选择 ID 类型 31

32 DAVID 使用选择功能类型 (GO,KEGG) 功能分析 32

33 DAVID 使用功能聚类集富集结果列表 33

34 MetaScape 使用上传基因选择物种提交运行 34

35 总结什么是转录组学? RNA-seq 的研究内容? 如何分析 RNA-seq 数据? 35

36 转录组分析谢谢! 浙江大学生命科学学院生物信息学实验室 36

PowerPoint 演示文稿

PowerPoint 演示文稿转录组学数据分析冯聪 ventson@zju.edu.cn 01/19-20, 2018 前言基因组学转录组学 from en.wikipedia 中心法则 : 遗传信息传递蛋白质组学 2 转录组概述 Protein-coding RNAs (mrnas) Transcriptome Non-coding RNAs (ncrnas) Structural RNAs: rrnas/trnas