美文网首页注释和富集
2021-06-29 转录数据的差异分析(DESeq2包)

2021-06-29 转录数据的差异分析(DESeq2包)

作者: 学习生信的小兔子 | 来源:发表于2021-06-29 16:29 被阅读0次

参考:https://www.jianshu.com/p/694b40961fc7

数据的生物学背景

糖皮质激素是炎症过程的有效抑制剂,由于其对气道平滑肌细胞的抗炎作用,被广泛用于治疗哮喘。但是分子机制是什么呢?本研究使用RNA-seq分析了四种不同的经地塞米松(一种合成糖皮质激素分子)处理的ASM细胞系的基因表达变化。他们发现了一些差异表达基因,并将地塞米松处理的ASM细胞与对照细胞进行了比较,但主要讨论的是一种名为CRISPLD2的基因。该基因编码一种已知参与肺发育的分泌蛋白,在以往的GWAS研究中,该基因中的snp与哮喘患者吸入皮质类固醇抵抗和支气管扩张剂反应有关。他们用qPCR证实了CRISPLD2 mRNA表达上调,用Western blotting证实了蛋白表达增加。

数据的结构


countData是基因计数矩阵

本教程中的基因计数矩阵文件是:airway_scaledcounts.csv。
一定要注意的是这里的基因计数矩阵一定是没有经处理的矩阵,因为DESeq2里面会进行差异分析统计。

colData是样本信息数据

其中id列需要与countData的第一行名字相符,顺序一样。
treatment列说明哪些是对照,哪些是处理。
sex列可有可无
本教程中样本信息数据是:airway_metadata.csv

导入数据

DESeq起作用的是一个叫做DESeqDataSet的对象。
加载tidyverse包,使用read_csv功能读入数据。这个对象包含了输入数据,中间计算像怎样均一化,还有差异表达分析的结果。
构建这个对象需要:基因表达矩阵,样本信息矩阵,还有一个指定实验是怎样设计的公式。

library(tidyverse)
mycounts <- read_csv("airway_scaledcounts.csv")
metadata <-  read_csv("airway_metadata.csv")
mycounts
mycounts
metadata
metadata

请注意此处的某些内容。元数据表中的样本 ID(SRR1039508、SRR1039509 等)与计数数据的列名完全匹配,但包含 Ensembl 基因 ID 的第一列除外

DESeqDataSet的特殊类型的对象。
这个对象包含输入数据,中间的处理数据,差异表达分析结果。
对于计数矩阵,我们需要构建DESeqDataSetFromMatrix对象。
这个对象包含内容有:
countData计数矩阵
colData 样本信息
design 实验设计
计数矩阵的列名是样本ID,必须和colData的行名字对应(或者当参数tidy = TRUE时,其第一列);
countData和colData必须是数据框格式;

转换为数据框格式

mycounts <- as.data.frame(mycounts)
metadata <- as.data.frame(metadata)
head(mycounts)
head(metadata)
class(mycounts)
class(metadata)

让我们检查一下count数据的列名(除了第一个,即ensgene)是否与colData中的id相同。

names(mycounts)[-1]
[1] "SRR1039508" "SRR1039509" "SRR1039512" "SRR1039513" "SRR1039516"
[6] "SRR1039517" "SRR1039520" "SRR1039521"
metadata$id
[1] "SRR1039508" "SRR1039509" "SRR1039512" "SRR1039513" "SRR1039516"
[6] "SRR1039517" "SRR1039520" "SRR1039521"
names(mycounts)[-1]==metadata$id
[1] TRUE TRUE TRUE TRUE TRUE TRUE TRUE TRUE
all(names(mycounts)[-1]==metadata$id)
[1] TRUE

现在计数矩阵和样品信息矩阵都有了,还需要一个design,这个指明了哪些是对照,哪些是处理。本例子中,dex这列告诉我们哪些样品是处理,哪些样品对照。我们将使用波浪号指定设计,如下所示:design=~dex。
PS:很多转录组数据中都会涉及到处理和对照。这里一定要小心,别写错了
创建的dds对象如下:

dds <- DESeqDataSetFromMatrix(countData=mycounts, 
                              colData=metadata, 
                              design=~dex, 
                              tidy=TRUE)
converting counts to integer mode
Warning message:
In DESeqDataSet(se, design = design, ignoreRank) :
  some variables in design formula are characters, converting to factors
dds <- DESeq(dds)
estimating size factors
estimating dispersions
gene-wise dispersion estimates
mean-dispersion relationship
final dispersion estimates
fitting model and testing

由于本例子的实验设计简单(单因素,两个分组, 处理组 vs 对照组)。可以使用results函数进行处理。如果更加复杂,就需要看看帮助文档?results。

res <- results(dds, tidy=TRUE)
res <- tbl_df(res)
res

tidy参数是为了输出结果中第一列的名字就是行名字

练习

使用管道符和arrange()函数对padj进行排序

res%>%arrange(padj)
或者可以这样
arrange(res,padj)
#R中的%>%类似Linux中的管道符,将上个命令的数据传递给下个命令。

可以看到默认排序结果是从小到大的

#将padj降序排序
res%>%arrange(desc(padj))

将注释信息合并到结果中

annotables_grch38.csv文件是基因的注释信息

#内连接
anno <- read_csv("annotables_grch38.csv")
View(anno)
res_anno = res %>%inner_join(anno,by=c("row"="ensgene"))
res_anno %>% filter(padj<0.05) #筛选padj小于0.05的值
res_anno %>% 
  filter(padj<0.05) %>% 
  write_csv("sigresults.csv")#将结果写入文件
anno
res
sigresults.csv

数据可视化

绘图计数

查看某个基因在感兴趣组组中表达

CRISPLD2 <- res_anno %>% filter(symbol=="CRISPLD2")
View(CRISPLD2)#CRISPLD2对的基因ID是ENSG00000103196
plotCounts(dds, gene="ENSG00000103196", intgroup="dex")

值得注意的是plotCount返回的不仅仅是图,也可以是数据

# Return the data
plotCounts(dds, gene="ENSG00000103196", intgroup="dex", returnData=TRUE)
              count     dex
SRR1039508  774.5002 control
SRR1039509 6258.7915 treated
SRR1039512 1100.2741 control
SRR1039513 6093.0324 treated
SRR1039516  736.9483 control
SRR1039517 2742.1908 treated
SRR1039520  842.5452 control
SRR1039521 6224.9923 treated

画个箱线图看看

# Plot it
plotCounts(dds, gene="ENSG00000103196", intgroup="dex", returnData=TRUE) %>% 
  ggplot(aes(dex, count)) + geom_boxplot(aes(fill=dex)) + scale_y_log10() + ggtitle("CRISPLD2")
#筛选有显著差异的基因  adj.P.Val意思是矫正后P值
foldChange=1.5 #fold change=1意思是差异是两倍
pvalue =0.01 #0.05

nrDEG = na.omit(res)
diff <- nrDEG
diffSig = diff[(diff$pvalue < pvalue & (diff$log2FoldChange>foldChange | diff$log2FoldChange<(-foldChange))),]
write.csv(diffSig, "DESeq2_diffSig_test.csv")

相关文章

网友评论

    本文标题:2021-06-29 转录数据的差异分析(DESeq2包)

    本文链接:https://www.haomeiwen.com/subject/vezhultx.html