生信入门(一)——DESeq2差异基因分析

生信(一)——DESeq2差异基因分析


记录学习过程,共勉。

一、差异基因分析原理

详见

二、代码实现

1、前提:安装DESeq2包

在这里插入图片描述

2.代码实现
setwd("D:\\RData");#设置编码位置
rt<-read.table("GSE149549_mRNA_Expression_Summary.txt",head=TRUE,row.names = 1)
head(rt)
#准备
rt<-as.matrix(rt) #将数据转换为矩阵格式
condition<-factor(c(rep("case",2),rep("control",3)),levels = c("case","control")) ## 设置分组信息,建立环境(5个样本,2组处理)
condition
coldata<-data.frame(row.names=colnames(rt),condition) 
#剔除无意义数据(不存在的基因read)
#nrow(rt)
i=1
while(i<=nrow(rt)){
  if(mean(rt[i,])<10)
    rt<-rt[-i,]
  i<-i+1;
};
#nrow(rt)

coldata #展示coldata内容

condition  #展示环境
#构建dds矩阵
dds<-DESeqDataSetFromMatrix(rt,coldata,design=~condition)
#标准化,对原始数据进行normalize
dds<-DESeq(dds)
dds
resultsNames(dds)
#使用deseq2包中的result()函数提取deseq2差异分析结果
res=results(dds,contrast = c("condition","case","control"))
head(res)
summary(res)
write.csv(res,file = "results.csv")
#分析结果可视化与导出
table(res$padj<0.05)
plotMA(res,ylim=c(-2,2))
mcols(res,use.names = TRUE)
#提取差异基因
res <- res[order(res$padj),]
resdata <-merge(as.data.frame(res),as.data.frame(counts(dds,normalize=TRUE)),by="row.names",sort=FALSE)
deseq_res<-data.frame(resdata)
up_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange > 1)) #提取上调差异表达基因
down_diff_result<-subset(deseq_res,padj < 0.05 & (log2FoldChange < -1)) #提取下调差异表达基因
write.csv(up_diff_result,"up_case_VS_control_diff_results.csv") #输出上调基因
write.csv(down_diff_result,"down_case_VS_control_diff_results.csv") #输出下调基因

#plot(res$log2FoldChange,res$padj)#绘制火山图

注:DESeq2的表达数据必须是read counts,不能使用标准化后的FPKM 或TPM
在这里插入图片描述
由图可见,表达量(基因的样本平均count)数值越大,离散度越小,差异基因倍数越少

#提取两种算法标准化后的表达量(VST和rlog,大样本通常使用VST,推荐使用)
vsd<-vst(dds,blind = FALSE)
rld<-rlog(dds,blind = FALSE)
#对照使用 log 2(n+1)的标准化方法
ntd<-normTransform(dds)
head(assay(vsd),3)
head(assay(rld),3)
head(assay(ntd),3)

在这里插入图片描述

#导出标准化后的数据(可用于GSEA等下游分析)
write.csv(as.data.frame(assay(vsd)),file="count_transformation_VST.csv")
write.csv(as.data.frame(assay(rld)),file="count_transformation_rlog.csv")

#绘制离散度散点图,离群值未做收缩
plotDispEsts(dds)

在这里插入图片描述

#对标准化后的数据进行可视化处理
#安装vsn包
#library(BiocManager)
#BiocManager::install("vsn")
#load package
#library("vsn")
#library(ggplot2)
#library(cowplot)

ntdl<-meanSdPlot(assay(ntd))
vsdl<-meanSdPlot(assay(vsd))
rldl<-meanSdPlot(assay(rld))

p1<-ntdl$gg+ggtitle("log2(n+1)")
p2<-vsdl$gg+ggtitle("VST")
p3<-rldl$gg+ggtitle("rlog")

plot_grid(p1,p2,p3,ncol = 3)

在这里插入图片描述
可见使用log2(n+1)标准化方法在低counts区域标准差较大,相反,rlog算法较大;而使用VST标准化后数据的标准差比较稳定;注意,虽然VST方法看起来效果最好,但图中的纵轴为方差的平方根(标准差),需要注意可能为由实验处理产生的真实方差。

#使用标准化后平均表达量top100的基因绘制热图
library("pheatmap")
select<-order(rowMeans(counts(dds,normalized=TRUE)),decreasing = TRUE)[1:100]
colData(dds)
df<-as.data.frame(colData(dds)[,c("condition","sizeFactor")])
mycolors=colorRampPalette(c("orange","white","purple"))(100)
pheatmap(assay(vsd)[select,],cluster_rows=FALSE,show_rownames=FALSE,border_color=NA,color=mycolors,cluster_cols=FALSE,annotation_col=df)

在这里插入图片描述

#绘制样本距离聚类热图
#计算样本距离矩阵
sampleDists<-dist(t(assay(vsd)))
#由向量转成矩阵
sampleDistMatrix<-as.matrix(sampleDists)
rownames(sampleDistMatrix)<-paste(vsd$condition,vsd$sizeFactor,seq="-")
colnames(sampleDistMatrix)<-NULL

library("RColorBrewer")
colors<-colorRampPalette(rev(brewer.pal(8,"YlGnBu")))(255)
pheatmap(sampleDistMatrix,clustering_distance_rows = sampleDists,clustering_distance_cols = sampleDists,col=colors)

在这里插入图片描述

三、小结

刚入门,继续学习,加油!

  • 19
    点赞
  • 191
    收藏
    觉得还不错? 一键收藏
  • 打赏
    打赏
  • 10
    评论
### 回答1: 好的,以下是 R 语言代码,分析单重复转录组数据并得出差异基因的示例: ``` library(edgeR) # 读入数据 counts <- read.table("gene_counts.txt", header=TRUE, row.names=1) # 建立edgeR计数对象 y <- DGEList(counts=counts) # 规范化计数数据 y <- calcNormFactors(y) # 对数据进行计数校正 y <- estimateCommonDisp(y) y <- estimateTagwiseDisp(y) # 对数据进行差异分析 design <- model.matrix(~0 + condition) fit <- glmFit(y, design) lrt <- glmLRT(fit, coef=2) # 得出差异基因 topTags(lrt) ``` 其中 `condition` 是指分析的样本的组别,`gene_counts.txt` 是基因计数数据的文件名。该代码使用了 `edgeR` 包,可以分析转录组数据并得出差异基因。 ### 回答2: 在分析单重复转录组数据的差异基因时,可以使用R语言编写代码来完成。下面是一个简单的代码示例: ```R # 导入所需包 library(edgeR) # 读取单重复转录组数据 countData <- read.csv("data.csv", header = TRUE, row.names = 1) # 创建一个DGE对象 dge <- DGEList(counts = countData, group = c("Single", "Repeat")) # 校正数据 dge <- calcNormFactors(dge) # 进行差异表达分析 dge <- estimateDisp(dge) fit <- glmFit(dge) # 设置比较组 contrast <- makeContrasts(Repeat - Single, levels = dge$groups) fit2 <- glmLRT(fit, contrast) # 筛选差异基因 diffGenes <- topTags(fit2, n = 100)$table # 打印差异基因结果 print(diffGenes) ``` 在这个示例中,我们首先导入所需的包,然后读取单重复转录组数据文件(假设为data.csv)。然后,我们创建了一个DGE对象,并对数据进行校正。接下来,进行差异表达分析,使用glmFit函数进行模型拟合,然后使用makeContrasts函数设置比较组,再使用glmLRT函数进行似然比检验。最后,我们使用topTags函数筛选出差异基因,并将结果打印出来。 请注意,这只是一个简单的示例代码,具体的分析方法和参数设置可能根据具体数据和需求而有所不同。因此,在实际分析中,可能需要根据具体情况进行适当的调整和修改。 ### 回答3: 使用R语言进行差异基因分析的一种常见方法是使用edgeR软件包。下面是一个基本的代码示例: ```R # 导入edgeR包 library(edgeR) # 读入表达矩阵数据,第一列是基因名字,后面的列是样本的表达值 expressionMatrix <- read.table("表达矩阵文件路径", header = T, row.names = 1) # 将表达矩阵转换为edgeR的数据对象 data <- DGEList(counts = expressionMatrix[, -1], genes = expressionMatrix[, 1]) # 进行标准化和归一化 data <- calcNormFactors(data) data <- estimateCommonDisp(data) # 设计实验(例如单重复转录组的对照组和实验组) design <- model.matrix(~ 0 + c("对照组","实验组")) # 进行差异基因分析 fit <- glmFit(data, design) lrt <- glmLRT(fit, coef = 1) topGenes <- topTags(lrt)$table # 根据给定的调节阈值(例如log2折变大于1,调整的p值小于0.05),选择差异基因 diffGenes <- topGenes[abs(topGenes$logFC) > 1 & topGenes$adj.PValue < 0.05, ] # 输出差异基因结果 write.table(diffGenes, file = "差异基因结果文件路径", sep = "\t", quote = F, row.names = F) ``` 以上代码首先导入edgeR包,并读入表达矩阵数据。然后,将表达矩阵转换为edgeR数据对象,并进行标准化和归一化处理。接下来,设置实验设计(对照组和实验组)并进行差异基因分析。最后,根据给定的调节阈值选择差异基因,并将结果输出至文件。 请注意,以上代码仅为基本示例,具体的数据处理步骤和参数设置应根据实际情况进行调整。
评论 10
添加红包

请填写红包祝福语或标题

红包个数最小为10个

红包金额最低5元

当前余额3.43前往充值 >
需支付:10.00
成就一亿技术人!
领取后你会自动成为博主和红包主的粉丝 规则
hope_wisdom
发出的红包

打赏作者

柚子味的羊

你的鼓励将是我创作的最大动力

¥1 ¥2 ¥4 ¥6 ¥10 ¥20
扫码支付:¥1
获取中
扫码支付

您的余额不足,请更换扫码支付或充值

打赏作者

实付
使用余额支付
点击重新获取
扫码支付
钱包余额 0

抵扣说明:

1.余额是钱包充值的虚拟货币,按照1:1的比例进行支付金额的抵扣。
2.余额无法直接购买下载,可以购买VIP、付费专栏及课程。

余额充值