我们拿到RNA-seq的测序结果,是一些打断了的reads片段,经过比对到参考基因组及reads计数等过程后,我们得到了表达量的数据,就是一个excel文档,里面有基因名,和很多很多的数字,那个就是表达量,记得刚开始做的时候我连表达量都不知道,拿到数据的时候都不知道我测序回来的结果原来就是一些不同数值的数字,我要从哪里开始做分析?
因为RNA是基于DNA转录来的,在有些条件下RNA根本不表达。所以为什么老师让你一定要在同一时间点取材,可能隔段时间,周围的环境发生变化,RNA表达水平也不同。RNA的表达水平就是表达量。
现在开始说得到表达量数据后如何做差异分析。
一、R包安装
1.1 常用软件包。
找差异基因要用到edgeR和DEGseq这两个R包, edgeR用来对得到的reads数进行归一化处理;DEGseq用来找差异基因。
基因表达量归一化:每个样本测序的总量不一样,要把它们处理到同一个数量级。
eg:A材料测序长度为10,B材料测序长度为100的,A表达量为8,B表达量为80,不能说B的表达量就比A高,因为它们不在一个数量级上
1.2 安装。
打开Bioconductor官网http://www.bioconductor.org/,在右上角Search栏输入edgeR,会弹出如下图的搜索结果,点击第一个Bioconductor-edgeR。(DEGseq的安装方法同edgeR,这里只介绍edgeR的安装)
打开Bioconductor-edgeR后,可看到详细的软件介绍、安装方法及下载地址。我们安装R包通常有两种方法:
1.直接在rstudio里面输入命令,R会自动转到下载地址自行安装,如下命令(##为解释说明)
2.手动安装。
找到网页最下面的Package Archives,根据你的电脑系统下载edgeR;打开Rstudio,最上边找到Tools,点击installpackages,Packagearchive中选你刚刚下载的edgeR;Installto library为你的安装路径。
二、代码
2.1 edgeR用法。
library(edgeR) #library为加载包的函数
library(limma)
#read.delim读取table的函数,“R-S.txt”为要进行归一化的表达量表格
x<- read.delim(“R-S.txt”,row.names=“GENEID”)
group<-factor(c(1,2)) #为得到的标准化因子
y<- DGEList(counts=x,group=group)
y<- calcNormFactors(y)
y$sample
如图所示,得到的nom.factors就是标准化因子
2.2 在得到标准化因子后到找差异基因还有一步,就是求标准值。
标准值=数值/10*标准化因子
将新得到的标准值的table保存在刚才的路径下,标为”R-Snorm.txt”
library(DEGseq) #加载DEGseq函数
geneExpFile <- system.file("extdata", “R-Snorm.txt", package="DEGseq")
geneExpMatrix1 <- readGeneExp(file=“R-Snorm.txt", geneCol=1, valCol=c(2))
geneExpMatrix2 <- readGeneExp(file=“R-Snorm.txt", geneCol=1, valCol=c(3))
DEGexp(geneExpMatrix1=geneExpMatrix1, geneCol1=1, expCol1=c(2:3), groupLabel1=“R",geneExpMatrix2=geneExpMatrix2, geneCol2=1, expCol2=c(2:3), groupLabel2=“S", method="MARS", outputDir="E:/R/station/RNAseq/1-2-8d")
# geneExpMatrix为建立矩阵,将原来的table读取为矩阵
#outputDir为输出结果的路径,之后在这个路径里会有一个output_score.txt文件,为得到的差异基因。
# |log2(Fold_change) normalized|>1(>2,根据具体实验情况决定)为标准化以后的差异倍数;Q-value<0.05。
欢迎关注生信人
TCGA | 小工具 | 数据库 |组装| 注释 | 基因家族 | Pvalue
基因预测 |bestorf | sci | NAR | 在线工具 | 生存分析 | 热图
生信不死 | 初学者 | circRNA | 一箭画心| 十二生肖 | circos
舞台|基因组 | 黄金测序 | 套路 | 杂谈组装 | 进化 | 测序简史