写在前面
最近在用的包经常涉及到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
获取吧!点个在看吧各位~ ✐.ɴɪᴄᴇ ᴅᴀʏ 〰
<center> <b>📍 往期精彩 <b> </center>
📍 <font size=1>🤩 WGCNA | 值得你深入学习的生信分析方法!~</font>
📍 <font size=1>🤩 ComplexHeatmap | 颜狗写的高颜值热图代码!</font>
📍 <font size=1>🤥 ComplexHeatmap | 你的热图注释还挤在一起看不清吗!?</font>
📍 <font size=1>🤨 Google | 谷歌翻译崩了我们怎么办!?(附完美解决方案)</font>
📍 <font size=1>🤩 scRNA-seq | 吐血整理的单细胞入门教程</font>
📍 <font size=1>🤣 NetworkD3 | 让我们一起画个动态的桑基图吧~</font>
📍 <font size=1>🤩 RColorBrewer | 再多的配色也能轻松搞定!~</font>
📍 <font size=1>🧐 rms | 批量完成你的线性回归</font>
📍 <font size=1>🤩 CMplot | 完美复刻Nature上的曼哈顿图</font>
📍 <font size=1>🤠 Network | 高颜值动态网络可视化工具</font>
📍 <font size=1>🤗 boxjitter | 完美复刻Nature上的高颜值统计图</font>
📍 <font size=1>🤫 linkET | 完美解决ggcor安装失败方案(附教程)</font>
📍 <font size=1>......</font>
本文由mdnice多平台发布
网友评论