RNA

2024-07-10 12:50| 来源: 网络整理| 查看: 265

============================================ 写在前面：可以参考另外一篇《得到差异基因后怎么做？》 ============================================

序走到这一步，似乎顺畅了许多，最主要时间不用花费那么多了。另外，以前曾经处理过不计其数的芯片，挑选差异表达基因，筛选关键基因，功能富集，还有基于全部数据的WGCNA(当然你也可以用差异基因来做，虽然不推荐，看不少文章也这么发），GSEA，PPI等等，这些后续我会慢慢发出来。但是，因为以前处理的芯片表达谱数据是符合正态分布，所以可以用t检验来筛选差异表达基因，但RNA-seq的read count普遍认为符合泊松分布。所以筛选DEGs的方法还是不一样 ------------要求--------------- 载入表达矩阵设置好分组信息用DEseq2进行差异分析输出差异分析结果来源于生信技能树 -------------------------参考文章----------------------------------- 2018.7月:Analyzing RNA-seq data with DESeq2 Count-Based Differential Expression Analysis of RNA-seq Data 可以下载的pdf havard,点此下载 -------------------开始------------------------ 正式载入数据之前，先看一下数据结构data structure，因为这是我们下面一切分析的起点。我们需要准备2个table：一个是countData，一个是colData two kinds of data.png 关于上面两个表的说明 countData表示的是count矩阵，行代表gene，列代表样品，中间的数字代表对应count数。colData表示sample的元数据，因为这个表提供了sample的元数据。 because this table supplies metadata/information about the columns of the countData matrix. Notice that the first column of colData must match the column names of countData (except the first, which is the gene ID column). colData的行名与countData的列名一致（除去代表gene ID的第一列） 1 载入数据（countData和colData） > library(tidyverse) > library(DESeq2) > #import data > setwd("F:/rna_seq/data/matrix") > mycounts head(mycounts) X control1 control2 treat1 treat2 1 ENSMUSG00000000001 1648 2306 2941 2780 2 ENSMUSG00000000003 0 0 0 0 3 ENSMUSG00000000028 835 950 1366 1051 4 ENSMUSG00000000031 65 83 52 53 5 ENSMUSG00000000037 70 53 94 66 6 ENSMUSG00000000049 0 3 4 5 #这里有个x，需要去除，先把第一列当作行名来处理 > rownames(mycounts) mycounts head(mycounts) control1 control2 treat1 treat2 ENSMUSG00000000001 1648 2306 2941 2780 ENSMUSG00000000003 0 0 0 0 ENSMUSG00000000028 835 950 1366 1051 ENSMUSG00000000031 65 83 52 53 ENSMUSG00000000037 70 53 94 66 ENSMUSG00000000049 0 3 4 5 # 这一步很关键，要明白condition这里是因子，不是样本名称；小鼠数据有对照组和处理组，各两个重复 > condition condition [1] control control treat treat Levels: control treat #colData也可以自己在excel做好另存为.csv格式，再导入即可 > colData colData condition control1 control control2 control treat1 treat treat2 treat 2构建dds对象,开始DESeq流程注释：dds=DESeqDataSet Object > dds dds # 查看一下dds的内容 > dds

显示为

class: DESeqDataSet dim: 6 4 metadata(1): version assays(3): counts mu cooks rownames(6): ENSMUSG00000000001 ENSMUSG00000000003 ... ENSMUSG00000000037 ENSMUSG00000000049 rowData names(21): baseMean baseVar ... deviance maxCooks colnames(4): control1 control2 treat1 treat2 colData names(2): condition sizeFactor 3 总体结果查看

接下来，我们要查看treat versus control的总体结果，并根据p-value进行重新排序。利用summary命令统计显示一共多少个genes上调和下调（FDR0.1）

> res = results(dds, contrast=c("condition", "control", "treat")) #或下面命令 > res= results(dds) > res = res[order(res$pvalue),] > head(res) > summary(res) #所有结果先进行输出 > write.csv(res,file="All_results.csv") > table(res$padj res = results(dds2, contrast=c("condition", "control", "treat")) > res = res[order(res$pvalue),] > head(res) log2 fold change (MLE): condition control vs treat Wald test p-value: condition control vs treat DataFrame with 6 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue padj ENSMUSG00000003309 548.1926 3.231611 0.2658125 12.157485 5.234568e-34 8.193146e-30 ENSMUSG00000046323 404.1894 3.067050 0.2628220 11.669687 1.820923e-31 1.425055e-27 ENSMUSG00000001123 341.8542 2.797485 0.2766499 10.112004 4.887441e-24 2.549941e-20 ENSMUSG00000023906 951.9460 2.382307 0.2510718 9.488551 2.342684e-21 9.116395e-18 ENSMUSG00000018569 485.4839 3.136031 0.3312999 9.465836 2.912214e-21 9.116395e-18 ENSMUSG00000000184 601.0842 -2.827750 0.3154171 -8.965112 3.099648e-19 8.085948e-16 > summary(res) out of 28335 with nonzero total read count adjusted p-value < 0.1 LFC > 0 (up) : 445, 1.6% LFC < 0 (down) : 625, 2.2% outliers [1] : 0, 0% low counts [2] : 12683, 45% (mean count < 18) [1] see 'cooksCutoff' argument of ?results [2] see 'independentFiltering' argument of ?results > write.csv(res,file="All_results.csv") > table(res$padj diff_gene_deseq2 1) #或 #> diff_gene_deseq2 1 | log2FoldChange < -1)) > dim(diff_gene_deseq2) > head(diff_gene_deseq2) > write.csv(diff_gene_deseq2,file= "DEG_treat_vs_control.csv")

结果显示如下：

> diff_gene_deseq2 1) > dim(diff_gene_deseq2) [1] 431 6 > head(diff_gene_deseq2) log2 fold change (MLE): condition control vs treat Wald test p-value: condition control vs treat DataFrame with 6 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue padj ENSMUSG00000003309 548.1926 3.231611 0.2658125 12.157485 5.234568e-34 8.193146e-30 ENSMUSG00000046323 404.1894 3.067050 0.2628220 11.669687 1.820923e-31 1.425055e-27 ENSMUSG00000001123 341.8542 2.797485 0.2766499 10.112004 4.887441e-24 2.549941e-20 ENSMUSG00000023906 951.9460 2.382307 0.2510718 9.488551 2.342684e-21 9.116395e-18 ENSMUSG00000018569 485.4839 3.136031 0.3312999 9.465836 2.912214e-21 9.116395e-18 ENSMUSG00000000184 601.0842 -2.827750 0.3154171 -8.965112 3.099648e-19 8.085948e-16 5 用bioMart对差异表达基因进行注释和RNA-seq(6): reads计数，合并矩阵并进行注释代码一样 library('biomaRt') library("curl") mart mart my_ensembl_gene_id mms_symbols head(mms_symbols) ensembl_gene_id external_gene_name description 1 ENSMUSG00000000120 Ngfr nerve growth factor receptor (TNFR superfamily, member 16) [Source:MGI Symbol;Acc:MGI:97323] 2 ENSMUSG00000000184 Ccnd2 cyclin D2 [Source:MGI Symbol;Acc:MGI:88314] 3 ENSMUSG00000000276 Dgke diacylglycerol kinase, epsilon [Source:MGI Symbol;Acc:MGI:1889276] 4 ENSMUSG00000000308 Ckmt1 creatine kinase, mitochondrial 1, ubiquitous [Source:MGI Symbol;Acc:MGI:99441] 5 ENSMUSG00000000320 Alox12 arachidonate 12-lipoxygenase [Source:MGI Symbol;Acc:MGI:87998] 6 ENSMUSG00000000708 Kat2b K(lysine) acetyltransferase 2B [Source:MGI Symbol;Acc:MGI:1343094] 5合并数据:res结果+mms_symbols合并成一个文件

合并的话两个数据必须有共同的列名，我们先看一下

> head(diff_gene_deseq2) log2 fold change (MLE): condition control vs treat Wald test p-value: condition control vs treat DataFrame with 6 rows and 6 columns baseMean log2FoldChange lfcSE stat pvalue padj ENSMUSG00000003309 548.1926 3.231611 0.2658125 12.157485 5.234568e-34 8.193146e-30 ENSMUSG00000046323 404.1894 3.067050 0.2628220 11.669687 1.820923e-31 1.425055e-27 ENSMUSG00000001123 341.8542 2.797485 0.2766499 10.112004 4.887441e-24 2.549941e-20 ENSMUSG00000023906 951.9460 2.382307 0.2510718 9.488551 2.342684e-21 9.116395e-18 ENSMUSG00000018569 485.4839 3.136031 0.3312999 9.465836 2.912214e-21 9.116395e-18 ENSMUSG00000000184 601.0842 -2.827750 0.3154171 -8.965112 3.099648e-19 8.085948e-16 > head(mms_symbols) ensembl_gene_id external_gene_name 1 ENSMUSG00000000120 Ngfr 2 ENSMUSG00000000184 Ccnd2 3 ENSMUSG00000000276 Dgke 4 ENSMUSG00000000308 Ckmt1 5 ENSMUSG00000000320 Alox12 6 ENSMUSG00000000708 Kat2b

可见，两个文件没有共同的列名，所以要先给'diff_gene_deseq2'添加一个‘ensembl_gene_id’的列名，方法如下:(应该有更简便的方法)

> ensembl_gene_id diff_gene_deseq2 colnames(diff_gene_deseq2)[1] diff_namehead(diff_name) #查看Akap8的情况 Akap8 ensembl_gene_id diff_gene_deseq2 colnames(diff_gene_deseq2)[1] diff_name head(diff_name) DataFrame with 6 rows and 9 columns ensembl_gene_id baseMean log2FoldChange lfcSE stat pvalue padj external_gene_name 1 ENSMUSG00000000120 434.04177 -1.293648 0.2713146 -4.768072 1.859970e-06 2.064698e-04 Ngfr 2 ENSMUSG00000000184 601.08417 -2.827750 0.3154171 -8.965112 3.099648e-19 8.085948e-16 Ccnd2 3 ENSMUSG00000000276 668.12500 -1.071362 0.2445381 -4.381168 1.180446e-05 9.603578e-04 Dgke 4 ENSMUSG00000000308 207.46719 1.944949 0.3427531 5.674489 1.391035e-08 3.819733e-06 Ckmt1 5 ENSMUSG00000000320 61.96266 1.451927 0.4637101 3.131109 1.741473e-03 4.105051e-02 Alox12 6 ENSMUSG00000000708 1070.03203 -1.046546 0.2049062 -5.107440 3.265530e-07 5.056107e-05 Kat2b description 1 nerve growth factor receptor (TNFR superfamily, member 16) [Source:MGI Symbol;Acc:MGI:97323] 2 cyclin D2 [Source:MGI Symbol;Acc:MGI:88314] 3 diacylglycerol kinase, epsilon [Source:MGI Symbol;Acc:MGI:1889276] 4 creatine kinase, mitochondrial 1, ubiquitous [Source:MGI Symbol;Acc:MGI:99441] 5 arachidonate 12-lipoxygenase [Source:MGI Symbol;Acc:MGI:87998] 6 K(lysine) acetyltransferase 2B [Source:MGI Symbol;Acc:MGI:1343094] > Akap8 Akap8 DataFrame with 1 row and 9 columns ensembl_gene_id baseMean log2FoldChange lfcSE stat pvalue padj external_gene_name 1 ENSMUSG00000024045 2281.053 -1.276329 0.2428315 -5.256028 1.471996e-07 2.775865e-05 Akap8 description 1 A kinase (PRKA) anchor protein 8 [Source:MGI Symbol;Acc:MGI:1928488] 至此，差异表达基因提取并注释完毕，下一步先进行数据可视化(Data visulization) 然后进行富集分分析及可视化后记：也可以用Y叔的ClusterProfiler或Annotation包进行基因名转换，很方便。

【本文地址】

公司简介

联系我们