参考教程：https://4va.github.io/biodatasci/r-rnaseq-airway.html#data_needed

数据来源

Himes et al. “RNA-Seq Transcriptome Profiling Identifies CRISPLD2 as a Glucocorticoid Responsive Gene that Modulates Cytokine Function in Airway Smooth Muscle Cells.” PLoS ONE. 2014 Jun 13;9(6):e99625. PMID: 24926665.

数据下载地址

• Length-scaled count matrix (i.e., countData): airway_scaledcounts.csv
• Sample metadata (i.e., colData): airway_metadata.csv
• Gene Annotation data: annotables_grch38.csv

数据的生物学背景

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

原始文章中的分析方式：

They did their analysis using Tophat and Cufflinks.

本教程采用DESeq2这个R包来分析 DESeq2.

要处理的数据结构

image.png

countData是基因计数矩阵

本教程中的基因计数矩阵文件是：airway_scaledcounts.csv。
一定要注意的是这里的基因计数矩阵一定是没有经处理的矩阵，因为DESeq2里面会进行差异分析统计。你事先如果进行过处理了，那还要DESeq2来噶啥。

colData是样本信息数据

其中id列需要与countData的第一行名字相符，顺序一样。

treatment列说明哪些是对照，哪些是处理。

sex列可有可无

本教程中样本信息数据是：airway_metadata.csv

导入数据

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

这个和read.csv是不一样的，要注意。我自己试试read.csv或者read.table()也可以读取，不过操作多了点参数，没这个方便。

library(tidyverse)

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

image.png

上图就是我们的基因计数矩阵

下面看看我们的样本信息矩阵

metadata

image.png

上图是我们的样品信息矩阵。
PS.对于这个信息矩阵，我们也可以自己打开excel自己把信息添加进去。这个教程给出了。我们直接用就好。

DESeqDataSet的特殊类型的对象。
这个对象包含输入数据，中间的处理数据，差异表达分析结果。

对于计数矩阵，我们需要构建DESeqDataSetFromMatrix对象。
这个对象包含内容有：
countData计数矩阵
colData 样本信息
design 实验设计

计数矩阵的列名是样本ID，必须和colData的行名字对应（或者当参数tidy = TRUE时，其第一列）；
countData和colData必须是数据库格式；
原文说法：

 The column names of countData are the sample IDs, and they must match the row names of colData (or the first column when tidy=TRUE). colData is an additional dataframe describing sample metadata. Both colData and countData must be regular data.frame objects – they can’t have the special tbl_df class wrapper created when importing with readr::read_*.

让我们再次查看mycounts和metadata

mycounts
metadata
class(mycounts)
class(metadata)

我们需要将这些数据转为常规的data.frame格式，以便它们能够很好地与DESeq2一起工作。

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]
metadata$id
names(mycounts)[-1]==metadata$id
all(names(mycounts)[-1]==metadata$id)

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

创建的dds对象如下：

dds <- DESeqDataSetFromMatrix(countData=mycounts, 
                              colData=metadata, 
                              design=~dex, 
                              tidy=TRUE)

如果看见下面警告，先别慌：
some variables in design formula are characters, converting to factors” don’t worry about it.

构建万dds对象后，并不完事了，这个对象没有处理，没有计算，就没法得到我们要的差异基因。我们需要使用DEseq()。DESeq()内部做了统计分析，我们得到差异基因。所以这也是为什么要用的计数矩阵是原始的计数矩阵，未被标准化处理过。

dds <- DESeq(dds)

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

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

image.png

注意：
tidy参数是为了输出结果中第一列的名字就是行名字。
你可以用View（）函数查看res.但是你可能注意到结果里面有很多NA值

练习：

Using a %>%, arrange the results by the adjusted p-value.
使用管道符和arrange()函数对padj进行排序

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

image.png

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

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

image.png

将注释信息合并到结果中

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")#将结果写入文件

最后总表

image.png

数据可视化

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

CRISPLD2 <- res_anno %>% filter(symbol=="CRISPLD2")
View(CRISPLD2)#CRISPLD2对的基因ID是ENSG00000103196

plotCounts(dds, gene="ENSG00000103196", intgroup="dex")

image

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

# Return the data
plotCounts(dds, gene="ENSG00000103196", intgroup="dex", returnData=TRUE)

image.png

Plot it

画个箱线图看看

plotCounts(dds, gene="ENSG00000103196", intgroup="dex", returnData=TRUE) %>% 
  ggplot(aes(dex, count)) + geom_boxplot(aes(fill=dex)) + scale_y_log10() + ggtitle("CRISPLD2")

image.png

MA & Volcano plots

挑选了上调或者下调的10个基因做火山图。火山图网上放的教程很多。我参考的是科研猫的代码，据说是cell期刊的配色。反正我看着挺受用的。

image.png