🤠 SE | 哇哦！让人不断感叹真香的数据格式！~

2023-02-24 本文已影响0人生信漫卷

写在前面

最近在用的包经常涉及到SummarizedExperiment格式的文件，不知道大家有没有遇到过。🤒

一开始觉得这种格式真麻烦，后面搞懂了之后发现真是香啊，爱不释手！~😜

什么是SummarizedExperiment

这种class主要包括了以下几个内容：👇

assay(), 包含一个或多个矩阵, 如行为基因名, 列为样本;

colData(), 对列的注释，格式为DataFrame, 如样本信息的描述;

rowData() 和（或）rowRanges(), 对列的注释，如rowRanges()描述基因坐标，rowData()描述每个基因差异分析后的logFC和pvalue;

metadata()，描述整个object的list;

用到的包

rm(list = ls())
library(SummarizedExperiment)
library(tidyverse)
library(RColorBrewer)
library(airway)
library(edgeR)

处理SummarizedExperiment对象

4.1 示例数据

这里我们用到airway包内的示例数据，讲解一下如何操作。🧐
这个data的基本研究设计是，用地塞米松处理人气道平滑肌细胞后进行RNA-seq。🤠

data(airway, package="airway")
se <- airway
se

4.2 查看colData

colData包含了样本或表型信息，返回的格式为DataFrame。🥰

colData(se)

4.3 提取colData的指定列

se$cell

4.4 查看列名和行名

我们看一下行名和列名。😉

colnames(se)

head(rownames(se))

4.5 查看表达矩阵名

assayNames(se)

4.6 查看表达矩阵

一个SummarizedExperiment格式的object是可以包含多个assay的。

assays(se)

4.7 查看指定assay

head(assay(se, "counts"))

4.8 rowRanges或granges

接下里是重中之重了，SummarizedExperiment允许代表不同特征的rowRanges (或granges)数据。🤩

length(rowRanges(se))

dim(se)

这里我们可以看到行特征对应了很多注释信息，这样我们在操作的时候就更加方便调取了。🤩

rowRanges(se)

4.9 获取start信息

start(rowRanges(se))

对于这种IRanges 对象，你也可以直接使用start()函数获取，其他常见的函数还有end，width。🤩

start(se)

4.10 提取制定对象

如果我们只想获取制定条件下的SummarizedExperiment对象，可以用subsetByOverlaps()函数，或者直接使用GRanges[List]。🤓

gr <- GRanges(seqnames = "1", ranges = IRanges(start = 1, end = 10^7))

subsetByOverlaps(airway, gr)

手动创建SummarizedExperiment

5.1 读入数据

这里我准备了样本数据和counts矩阵两个文件，大家跟着我一起试一下吧。

pdata <- read.csv("./SummarizedExperiment/airway-sample-sheet.csv") 
counts <- read.csv("./SummarizedExperiment/airway-read-counts.csv")

5.2 整理数据并创建SummarizedExperiment

pdata <- column_to_rownames(pdata, "Run")
counts <- column_to_rownames(counts, "Run")
se_juan <- SummarizedExperiment(t(counts), colData = pdata, rowRanges = )
se_juan

5.3 准备rowData

我们再试着把rowData加进之前的SummarizedExperiment里。🤠

这里我们用一下EnsDb.Hsapiens.v86包来获取基因的各种信息，如染色体位置、起止位点、类型、id等等，这个包以后我们再具体讲怎么用。

输出的文件为Granges，完美匹配。😁

library(EnsDb.Hsapiens.v86)

edb <- EnsDb.Hsapiens.v86
filter <- rownames(se_juan)
genes <- genes(edb)
genes <- genes[genes$gene_id %in% filter]
head(genes)

5.4 添加rowData

这里需要说一下，有的基因没有具体的位点信息等，可能和版本有关系，以后我们再讲怎么处理。🥰

rowData(se_juan) <- genes
se_juan
rowData(se_juan)

小练习

我们做个小练习，试试画个基因平均表达的boxplot吧，还要取一下log哦。😏

assay(se_juan)  %>% 
  log() %>% 
  boxplot(col = colorRampPalette(brewer.pal(8, "Set2"))(8))

<center>最后祝大家早日不卷!~</center>

需要示例数据的小伙伴，在公众号回复SummarizedExperiment获取吧！

点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰

本文由mdnice多平台发布