2021-06-29 转录数据的差异分析(DESeq2包)
数据的生物学背景
糖皮质激素是炎症过程的有效抑制剂，由于其对气道平滑肌细胞的抗炎作用，被广泛用于治疗哮喘。但是分子机制是什么呢?本研究使用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")
