引言

❝
大规模检测 circRNA 的方式有芯片测序和二代高通量测序, 芯片测序的原理是针对已知的 反向可变剪接位点 来设计好探针, 然后进行检测, 而二代测序则是对细胞内产生的所有 RNA 分子通过特定的比对软件和 反向可变剪接位点 识别来鉴定 circRNA, 这种可以鉴定出一些新的环形 RNA 分子。
对于芯片测序, 美国 Arraystar 公司在首款 circRNA 芯片基础上迅即升级版本为 V2.0。其 circRNA 来源融合了环状 RNA 研究的高水平文献，所有 cicrRNA 都经过了严谨的实验验证，以便于对不同生理及病理条件下的 circRNA 进行系统的研究。同时我们对所有差异表达的 circRNA 用高匹配值的 miRNA 靶标位点进行了标注。当然还有些其他公司。

我们 GEO 找找 circRNA 芯片的数据, 第一个最新的数据就是 Arraystar 的:

10 个样本, 5 个处理组和 5 个对照组:

数据处理方法也是比较简单:

分析

❝
我们简单分析一下这个数据, 举个例子。使用 GEOquery 来下载数据。
setwd(“~/JunJun/10.array_circRNA_proj/”)

BiocManager::install(“GEOquery”)

library(GEOquery)
library(limma)
library(dplyr)

download data

gse <- getGEO(GEO = “GSE228477”,AnnotGPL = T,destdir = “./”)
gse

$GSE228477_series_matrix.txt.gz

ExpressionSet (storageMode: lockedEnvironment)

assayData: 13617 features, 10 samples

element names: exprs

protocolData: none

phenoData

sampleNames: GSM7122512 GSM7122513 … GSM7122521 (10 total)

varLabels: title geo_accession … tissue:ch1 (34 total)

varMetadata: labelDescription

featureData

featureNames: ASCRP3000001 ASCRP3000002 … ASCRP3013617 (13617 total)

fvarLabels: ID circRNA … SEQUENCE (5 total)

fvarMetadata: Column Description labelDescription

experimentData: use ‘experimentData(object)’

Annotation: GPL21825

表型数据:

get phenoData

pdata <- pData(gse[[1]])

check

head(pdata[1:3,1:3])

title geo_accession status

GSM7122512 M1 GSM7122512 Public on Apr 20 2024

GSM7122513 M2 GSM7122513 Public on Apr 20 2024

GSM7122514 M3 GSM7122514 Public on Apr 20 2024

check disease

table(pdata$characteristics_ch1.1)

disease: benign biliary obstruction disease: malignant biliary obstruction

5 5

提取表达矩阵:

extarct expression matrix

exp <- exprs(gse[[1]])

boxplot

boxplot(exp,outline = FALSE, notch = T,col = group_list, las = 2)

Quantile normalization 和 ID 转换:

Quantile normalization

exp_norm <- normalizeBetweenArrays(object = exp,method=”quantile”)

load annotation

gpl <- getGEO(filename = “GPL21825.soft.gz”) anno <- gpl@dataTable@table %>%
dplyr::filter(ID %in% rownames(exp_norm))

check

head(anno,3)

ID circRNA TYPE BUILD

1 ASCRP3000001 hsa_circRNA_082501 circRNA HG19

2 ASCRP3000002 hsa_circRNA_407247 circRNA HG19

3 ASCRP3000003 hsa_circRNA_007418 circRNA HG19

SEQUENCE

1 AAAAAAACAACCAAAAAATGTTCAACAGCATGAGAAGGTTCAGAAAGCCAGTACAGAGGG

2 AAAAAACGAAGAAAAAGAGACACCCAGCTCACCTCCAAGTTTGCCTGCAGGAGCCGGCTC

3 AAAAAACTAGAAATGTGTTCAGAAATTAAAGGTCCACAGAAGGAGGGCCTGTCCTCCCAA

check ID order

identical(x = anno$ID,y = rownames(exp_norm))

[1] TRUE

ID transform

rownames(exp_norm) <- anno$circRNA

check

head(exp_norm,3)

GSM7122512 GSM7122513 GSM7122514 GSM7122515 GSM7122516 GSM7122517 GSM7122518

hsa_circRNA_082501 5.708676 5.731464 7.206074 5.837747 5.790620 5.446059 5.644320

hsa_circRNA_407247 10.063683 10.419590 10.672558 10.525109 10.343062 11.418045 10.387986

hsa_circRNA_007418 9.026643 8.764032 8.121701 8.417818 8.471591 8.068711 8.720491

GSM7122519 GSM7122520 GSM7122521

hsa_circRNA_082501 5.924437 5.916646 6.247908

hsa_circRNA_407247 10.027347 9.905763 10.726633

hsa_circRNA_007418 8.003765 7.867321 8.264380

下面就是熟悉的差异分析了:

==============================================================================

diff analysis

==============================================================================

group_list <- factor(rep(c(“tumor”,”normal”),each = 5),levels = c(“normal”,”tumor”))
group_list

[1] tumor tumor tumor tumor tumor normal normal normal normal normal

Levels: normal tumor

design <- model.matrix(~group_list) fit <- lmFit(exp_norm,design) fit <- eBayes(fit) deg <- topTable(fit,coef = 2,number = Inf) %>%
dplyr::mutate(type = dplyr::case_when(logFC >= 1 & P.Value < 0.05 ~ “sigUp”,
logFC <= -1 & P.Value < 0.05 ~ “sigDown”,
.default = “nonSig”))

chcek

head(deg,3)

logFC AveExpr t P.Value adj.P.Val B type

hsa_circRNA_056558 4.320917 10.141836 10.428270 4.272850e-07 0.00581834 6.204301 sigUp

hsa_circRNA_100633 2.912608 10.057436 8.666649 2.730792e-06 0.01262898 4.763709 sigUp

hsa_circRNA_105031 1.999955 8.933745 8.356502 3.900703e-06 0.01262898 4.472719 sigUp

summary

table(deg$type)

nonSig sigDown sigUp

13304 41 272

画个火山图:

==============================================================================

volcanoplot

==============================================================================

library(ggplot2)

ggplot(deg,aes(x = logFC,y = -log10(P.Value))) +
geom_point(aes(color = type)) +
geom_hline(yintercept = -log10(0.05), linetype = “dashed”) +
geom_vline(xintercept = c(-1, 1), linetype = “dashed”) +
xlim(-5,5) +
theme_bw() +
theme(panel.grid = element_blank()) +
scale_color_manual(values = c(“sigUp” = “red”,”sigDown” = “blue”,”nonSig” = “grey85”))

❝
上面转换后的 ID 大家可以去 circRNA 相关数据库查询, 比如 circBase, 来获得更具具体的信息。
结尾

❝
路漫漫其修远兮,吾将上下而求索。
欢迎加入生信交流群。加我微信我也拉你进 微信群聊 老俊俊生信交流群 (微信交流群需收取 20 元入群费用,一旦交费,拒不退还!(防止骗子和便于管理)) 。QQ 群可免费加入, 记得进群按格式修改备注哦。