环状RNA 4分文章路线及代码-GEO芯片数据
环形RNA是mRNA在剪接的过程中,上游exon的5’端与下游exon的3’端剪接到一起,从而形成的首尾相接的环状RNA分子。大多数定位于细胞质, 且序列高度保守。比线性RNA更稳定, 不易被降解。许多circRNA表达水平与线性RNA的表达水平相当, 一些circRNA的表达水平甚至超过它们的线性异构体10倍。大多数来自蛋白质编码基因的外显子。大部分是非编码RNA(noncoding RNA, ncRNA)。与miRNA相互作用(MRE)。circRNA能与蛋白质结合, 抑制蛋白质活性、募集蛋白质复合体的组分或调控蛋白质的活性。一些circRNA可作为翻译的模板指导蛋白质的合成。circRNA通过碱基互补配对直接调控其他RNA水平。

找到差异环状RNA
挑选感兴趣环状RNA,在circBase里面查阅,在什么基因上,什么外显子上,序列是什么,基因组上什么位置,名字需要转换为标准的名字,也就是这个数据库的名字,得到信息后,就可以得到基因名字,序列,基因位置,可以做一个圈图。圈图包括ORF(可以编码蛋白),带三角部分,可以和miRNA结合,另外的,可以跟蛋白结合,可以和mRNA结合。
数据下载
GEO NCBI官网,https://www.ncbi.nlm.nih.gov/geo/ ,搜索关键字,会显示出现多少结果。

选择人的标本,选择芯片数据,符合要求的,如circRNA肝癌标本芯片检测。

可以选择一个数据进行分析,也可以选择两个或者多个进行联合分析,联合分析的好处是,得到的差异circRNA较少,后续分析工作量较少。CircRNA的数据库不是很完善,可能需要进行多个数据库的分析,所以最好选择联合分析。也可以分别做差异,取交集。这里解释的是单芯片分析,多芯片联合分析则以后再讲解。
把页面上的Matrix文件下载,同时下载平台文件,放于文件夹,有了平台文件,就可以将矩阵文件转化为我们所需要的circRNA ID文件。

对数据进行注释
探针名转换为circRNA的名字,一种是直接转换完成的,可以直接用。但有时候下载的数据是Agilent形式,has_circRNA_000065,而不是has_circ_000065的circBase标准名字,就要进行两次转换(如果是别的平台,名字转换比较麻烦,可以找作者咨询)。
关于此步转换,很多人都会遇到问题,就算是找到了代码,也经常出错。我们来详细解释一下此步骤。1,ID转换,就算为了得到标准的circRNA名字,并对应表达量从而得到我们后续需要进行分析的matrix(矩阵数据,就是一个Excel表格,一点也不神秘)。2,我们需要准备参考文件,也就是说,我们有5个ID,那么参考库是什么?参考库里面包括了这5个ID并且有其对应的标准circRNA ID。3,平台文件里面,即platform,有探针ID和安捷伦自己的ID,即为第一步的参考文件。那么安捷伦自己的ID,如何再转换成标准的circRNA ID,则需要另外一个参考文件,图如下。

我们想用代码,除非自己懂代码,不然不知道哪一步会遇到问题,这里我选择用Excel做ID转换。
下图为平台文件

下图为下载的matrix文件,就是需要进行分析的文件

将文件放在同一个Excel里面的不同sheet

调用函数VLOOKUP
查找值为A2,自己根据具体情况选择。

第一步参考文件是platform,选择查找和筛选数据。如下图,从第一列查找A2的内容,从第二列筛选出来需要转换的内容。

复制黏贴公式,第一步转换就结束了。
下面进行第二步转换,继续用VLOOKUP公式。



复制公式,两部转换就全部完成了,但是上图显示,有的没有找到,则为找不到相关的信息,没有对应数据,我们可以手动删除。
ID文件可以去网盘下载,见文末最后链接。
下一步则在R里面进行分析,安装完R以后,进行以下操作。
安装R包
R包包括多种,如TCGAbiolinks,limma,这里我们用limma,根据不同需要选择不同包。
if (!requireNamespace("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("limma") #安装limma
设置参数
logFoldChange=2 #设置log2的值,如果=2,则为4倍为标准
adjustP=0.05 #设置校正后的P值,这里为0.05,如果筛选出来的基因较多,也可以选择0.01或者0.001,把标准提高。
library(limma) #载入limma包(如果没有报错,说明安装好了)
#setwd("C:\\Users\\ \\Desktop\\circRNA\\breast.diff") #设置工作目录,一定要双斜杠。此步骤在bioinformatics环境中无法设置,可以直接上传数据文件,然后进行下一步操作。如下图。
rt=read.table("geneMatrix2.txt",sep="\t",header=T,check.names=F) #读取文件,geneMatrix2.txt为文件名,包括后缀也要写进去,sep是分割,以\t分割,header=T是指有表头,T即为TURE,表示需要显示,如果为FALSE,则是要求结果中不显示,是根据我们自己的选择而定。Check.name是指不对名字进行检测。
rt=as.matrix(rt) #转化为矩阵

以下代码是指,如果一个circRNA有多行数据,则取均值。
rownames(rt)=rt[,1]
exp=rt[,2:ncol(rt)]
dimnames=list(rownames(exp),colnames(exp))
rt=matrix(as.numeric(as.matrix(exp)),nrow=nrow(exp),dimnames=dimnames)
rt=avereps(rt)
校正
rt=normalizeBetweenArrays(as.matrix(rt)) #对芯片数据进行校正
#rt=log2(rt+1) 如果数据没有取log2,则需要运行这条代码,判断数据有没有取log2,可以根据下载过程中的指示,同时如果数据为上千的值,如1205.998,则是未取log2的。
以下则开始进行差异表达
所有的差异情况出来以后,我们根据我们的需要设计条件,根据我们的过滤标准进行挑选。
modType=c(rep("con",7),rep("treat",7)) #对分组进行定义,control组和处理组分别有7个。在数据里面,也要把分组按照顺序排列,如con排在前面,treat排在后面。
design <- model.matrix(~0+factor(modType)) #设计比较方案
colnames(design) <- c("con","treat") #设计比较方案
fit <- lmFit(rt,design) #计算各分组的均值和SD
cont.matrix<-makeContrasts(treat-con,levels=design) #处理组减去对照组,就是倍数,因为已经取了log2。
fit2 <- contrasts.fit(fit, cont.matrix) #计算结果
fit2 <- eBayes(fit2) #计算结果
allDiff=topTable(fit2,adjust='fdr',number=200000) #输出数据,topTable就是要把前多少位输出来。这里定义了200000,校正方法是fdr。
write.table(allDiff,file="limmaTab.xls",sep="\t",quote=F) #输出保存数据,保存在limmaTab.xls中。
diffSig <- allDiff[with(allDiff, (abs(logFC)>logFoldChange & adj.P.Val < adjustP )), ] #前面我们已经定义了logFoldChange=2和adjustP=0.05,这里设置过滤条件为logFC>logFoldChange & adj.P.Val < adjustP ,过滤条件就设置好了
write.table(diffSig,file="diff.xls",sep="\t",quote=F) #结果输出到diff.xls。
以下代码是将高表达和地表达分开,分别输出为up.xls和down.xls。
diffUp <- allDiff[with(allDiff, (logFC>logFoldChange & adj.P.Val < adjustP )), ]
write.table(diffUp,file="up.xls",sep="\t",quote=F)
diffDown <- allDiff[with(allDiff, (logFC<(-logFoldChange) & adj.P.Val < adjustP )), ]
write.table(diffDown,file="down.xls",sep="\t",quote=F)
以下代码为输出差异基因的表达量。输出为heatmap.txt文件。
hmExp=rt[rownames(diffSig),]
diffExp=rbind(id=colnames(hmExp),hmExp)
write.table(diffExp,file="heatmap.txt",sep="\t",quote=F,col.names=F)
以下代码为绘制火山图代码。
tiff(file="vol.tiff",
width = 12, #图片的宽度
height =12, #图片的高度
units ="cm",
compression="lzw",
bg="white",
res=600) #清晰度,可以根据需要进行修改,文章有的有需要。一般为300-1000。
xMax=max(-log10(allDiff$adj.P.Val))
yMax=max(abs(allDiff$logFC))
plot(-log10(allDiff$adj.P.Val), allDiff$logFC, xlab="-log10(adj.P.Val)",ylab="logFC",
main="Volcano", xlim=c(0,xMax),ylim=c(-yMax,yMax),yaxs="i",pch=20, cex=0.8) #先把所有的差异表达打点。
diffSub=subset(allDiff, adj.P.Val<adjustP & logFC>logFoldChange) #对上调的进行打点。
points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="red",cex=0.8) #打点为红色。
diffSub=subset(allDiff, adj.P.Val<adjustP & logFC<(-logFoldChange)) #对下调的进行打点。
points(-log10(diffSub$adj.P.Val), diffSub$logFC, pch=20, col="green",cex=0.8) #下调的打点为绿色。
abline(h=0,lty=2,lwd=3)
dev.off()
大家需要修改的部分,一个是工作目录,一个是分组情况,一个是输入文件。
以下代码是画热图,可以根据自己的需要进行选择数据,比如只想画其中某些circRNA,自行选择。
install.packages("pheatmap") #安装heatmap包,需要选择镜像,选国内
setwd("C:\\Users\\Desktop\\circRNA\\06.pheatmap") #设置工作目录
rt=read.table("heatmap.txt",sep="\t",header=T,row.names=1,check.names=F) #读取数据表格,对第一列作为行名,不需对名字进行检查。
library(pheatmap) #加载包
Type=c(rep("con",7),rep("treat",7)) #修改对照和处理组样品数目,下载数据的时候知道对照和处理组的数目。
names(Type)=colnames(rt)
Type=as.data.frame(Type) #样品分组界定
以下代码为画图代码
tiff(file="heatmap.tiff",
width = 20, #图片的宽度,宽度和高度需要调整,如果样品很多,那么宽度就要宽一点,不然无法显示全。如果差异基因较多,高度要调高一点。
height =25, #图片的高度,同上。
units ="cm", #单位
compression="lzw",
bg="white",
res=600) #清晰度
pheatmap(rt, #前面定义的数据
annotation=Type,
color = colorRampPalette(c("green", "black", "red"))(50), #低表达为绿色,中表达为黑色,高表达为红色。50为颜色分为多少份。
cluster_cols =F, #聚类,如果是改为TURE,有可能分组会乱
fontsize = 8, #定义字体
fontsize_row=4, #字体大小
fontsize_col=8) #字体大小
dev.off() #关闭图形,即认为不再需要编辑。

在画热图的时候出现下列错误。正在解决。
已解决,上一步代码忘记修改每组的个数,我分析的数据是每组5个,而原代码里面的分组数是7个。
大家需要修改的部分,一个是工作目录,一个是分组情况,一个是输入文件。
circBase数据库
此数据库可以找到目标环状RNA的序列,所在基因,以及在染色体上的位置。但是我们分析出来的环状RNA太多,不可能全部分析,可以如何确定我们要深入研究的环状RNA呢?1)实验验证,拿前10个取验证,得到2-3个;2)多找几个芯片,可以ronk,可以取交集,剩下几个环状RNA,做后续的分析;3)只做前5-10个,后续做富集分析,看看哪个环状RNA跟我们的符合,比如和肝癌相关的环状RNA。
可以在circBase查到环状RNA的序列,所在基因,染色体位置,正负链,转录本,在什么样本中存在,来自于哪一项研究等信息。在下载序列时,选择spliced,即为成环的序列。
拿到序列,可以查看此环状RNA可以和什么miRNA结合,也可以看一下,可以和什么蛋白结合,查看开放阅读框,看它是否有编码蛋白的能力。
CircRNA圈图的绘制,环状RNA可能可以编码蛋白,也可能可以和蛋白结合,也可能可以和miRNA结合,利用CSCD circRNA软件绘制,工具很多,可以自由选择,可以根据文献进行选择。根据下图,可以选择是肿瘤特有的,还是正常特有的,还是共有的circRNA。选择其在哪个基因上,细胞的位置,环状RNA的位置。根据环状RNA位置,我们可以先选择肿瘤特有的,进行搜索,如果没搜索到,我们可以搜索一下common组,右下角即为我们需要的圈图,如果找不到该circRNA,则需要换一个数据库。圈图右上角可以放大,下载,红色部分为和miRNA结合部分,蓝色为和蛋白结合部分,绿色为开放阅读框,数字代表exon,6,8,10可以组成一个环状RNA,环状RNA结合miRNA的结合情况都显示了。
与环状RNA结合的蛋白,基本都是转录因子,clip-seq的研究可以显示RNA和蛋白的结合情况,上图中可以找出结合位置,以及其他的一些信息,重点关注与哪个蛋白结合。做实验比较多结合蛋白研究,生信分析较少,做富集可能得到的结果较少,但是如果实验验证了环状RNA结合蛋白并发挥功能,则可以发比较好的文章。但是环状RNA可以和miRNA结合,miRNA可以和mRNA结合,预测环状RNA功能,作为miRNA海绵来发挥功能。
开放阅读框,环状RNA是否编码蛋白,氨基酸的信息是什么,还是在圈图分析中查看结果,结果显示,得到的蛋白序列超过了环状RNA的序列,所以,没办法直接得到相关准确信息,我们继而采用ORFfinder,输入序列,最小序列选择75,genetic code选择standard,选择ATG或者其他启动密码子,提交,则显示结果,会显示几个,及在正链还是负链,长度(nt),点mark,可以显示氨基酸的序列,此部分只是预测,如果实验验证发现编码成功,则可以进一步研究该蛋白的功能。
预测miRNA的靶基因,miRTarBase、miRDB和TargetScan这三个数据库进行预测,可以通过代码进行筛选靶基因。
GO和KEGG分析,上述分析,需要转换基因ID,目前没办法直接用基因的名字进行分析,因为基因的名字可能有别名。ID是唯一的,不存在别名的问题。安装ID转换包。
install.packages("RSQLite") #安装包
source("https://bioconductor.org/biocLite.R") #定义包来源
biocLite("org.Hs.eg.db") #安装包
setwd("C:\\Users\\Desktop\\circRNA\\symbol2id") #设置工作目录
library("org.Hs.eg.db") #调用包
rt=read.table("target.txt",sep="\t",check.names=F,header=T) #读取文件
genes=as.vector(rt[,2]) #基因转化为向量
entrezIDs <- mget(genes, org.Hs.egSYMBOL2EG, ifnotfound=NA) #找基因对应的ID,如果存在,就显示ID,如果不存在,就显示NA
entrezIDs <- as.character(entrezIDs) #ID转化为字符
out=cbind(rt,entrezID=entrezIDs) #ID加入到我们的rt里面,前面已经对rt进行了定义。
write.table(out,file="id.txt",sep="\t",quote=F,row.names=F) #输出文件
GO富集分析
install.packages("colorspace") #安装包
install.packages("stringi") #安装包
source("http://bioconductor.org/biocLite.R")
biocLite("DOSE") #安装包
biocLite("clusterProfiler") #安装包
biocLite("pathview") #安装包
setwd("C:\\Users\\Desktop\\circRNA\\GO") #设置工作目录
library("clusterProfiler") #引用包
library("org.Hs.eg.db") #引用包
rt=read.table("id.txt",sep="\t",header=T,check.names=F) #读取文件
rt=rt[is.na(rt[,"entrezID"])==F,] #去掉基因ID为NA的条目
sum=rt$Sum #提取ID
gene=rt$entrezID #提取ID
names(sum)=gene #提取ID
kk <- enrichGO(gene = gene,OrgDb = org.Hs.eg.db, pvalueCutoff =0.05, qvalueCutoff = 0.05) #GO富集分析,OrgDb = org.Hs.eg.db后台数据库
write.table(kk,file="GO.txt",sep="\t",quote=F,row.names = F) #输出文件
tiff(file="barplot.tiff",width = 45,height = 25,units ="cm",compression="lzw",bg="white",res=600) #柱状图,可以设置宽度和高度,单位为cm,分辨率为600
barplot(kk, drop = TRUE, showCategory =100) #绘图,显示100个通路
dev.off()
tiff(file="dotplot.tiff",width = 45,height = 25,units ="cm",compression="lzw",bg="white",res=600) #气泡图,可以设置宽度和高度,单位为cm,分辨率为600
dotplot(kk,showCategory = 100) #绘图,显示100个通路
dev.off()
install.packages("colorspace") #安装包
install.packages("stringi") #安装包
source("http://bioconductor.org/biocLite.R")
biocLite("DOSE") #安装包
biocLite("clusterProfiler") #安装包
biocLite("pathview") #安装包
setwd("C:\\Users\\Desktop\\circRNA\\GO") #设置工作目录
library("clusterProfiler") #引用包
library("org.Hs.eg.db") #引用包
rt=read.table("id.txt",sep="\t",header=T,check.names=F) #读取文件
rt=rt[is.na(rt[,"entrezID"])==F,] #去掉基因ID为NA的条目
sum=rt$Sum #提取ID
gene=rt$entrezID #提取ID
names(sum)=gene #提取ID
kk <- enrichGO(gene = gene,OrgDb = org.Hs.eg.db, pvalueCutoff =0.05, qvalueCutoff = 0.05) #GO富集分析,OrgDb = org.Hs.eg.db后台数据库
write.table(kk,file="GO.txt",sep="\t",quote=F,row.names = F) #输出文件
tiff(file="barplot.tiff",width = 45,height = 25,units ="cm",compression="lzw",bg="white",res=600) #柱状图,可以设置宽度和高度,单位为cm,分辨率为600
barplot(kk, drop = TRUE, showCategory =100) #绘图,显示100个通路
dev.off()
tiff(file="dotplot.tiff",width = 45,height = 25,units ="cm",compression="lzw",bg="white",res=600) #气泡图,可以设置宽度和高度,单位为cm,分辨率为600
dotplot(kk,showCategory = 100) #绘图,显示100个通路
dev.off()
最后显示出GO的基因,再把基因ID转化为基因名。气泡图和柱状图一半只放一个。GO.txt即为输出的结果,可以有功能分析,基因名相关信息,基因ID转换为基因名,间接的可以知道环状RNA是与哪些功能相关,大部分还是和转录相关。
KEGG通路富集分析,除了得到气泡图和柱状图,还可以得到通路图,不同于GO,KEGG更注重于通路分析,需要用到基因ID文件,KEGG分析需要输入logFC,但是我们此处没有做基因的差异表达,仅以Sum(文件里的一栏数字)作为分析值,为了得到通路。
install.packages("colorspace") #安装包(前面如果安装过了,无需重复安装,下同),安装包已经改变安装方式了,使用BioManager,所以本教程里的安装包方法已取消,但是在本站的docker中,包几乎已经全部安装好,可以直接调用。
install.packages("stringi") #安装包安装包已经改变安装方式了,使用BioManager,所以本教程里的安装包方法已取消,但是在本站的docker中,包几乎已经全部安装好,可以直接调用。
source("http://bioconductor.org/biocLite.R") #安装包安装包已经改变安装方式了,使用BioManager,所以本教程里的安装包方法已取消,但是在本站的docker中,包几乎已经全部安装好,可以直接调用。
biocLite("DOSE") #安装包安装包已经改变安装方式了,使用BioManager,所以本教程里的安装包方法已取消,但是在本站的docker中,包几乎已经全部安装好,可以直接调用。
biocLite("clusterProfiler") #安装包安装包已经改变安装方式了,使用BioManager,所以本教程里的安装包方法已取消,但是在本站的docker中,包几乎已经全部安装好,可以直接调用。
biocLite("pathview") #安装包安装包已经改变安装方式了,使用BioManager,所以本教程里的安装包方法已取消,但是在本站的docker中,包几乎已经全部安装好,可以直接调用。
setwd("C:\\Users\\Desktop\\circRNA\\KEGG") #设置工作目录,用docker不需要设置目录。
library("clusterProfiler") #引用包
rt=read.table("id.txt",sep="\t",header=T,check.names=F) #读取文件
rt=rt[is.na(rt[,"entrezID"])==F,] #删除NA行,同上
sum=rt$Sum #提取ID
gene=rt$entrezID #提取ID
names(sum)=gene #提取ID
kk <- enrichKEGG(gene = gene, organism = "hsa", pvalueCutoff =0.05, qvalueCutoff =0.05) #kegg富集分析 ,本步骤遇到问题,由于开了VPN,导致访问远程网络,下载的时候网络不行,回到结果是no mapped,如果没有下载下来数据,是不能对比的,所以显示no mapped。
write.table(kk,file="KEGG.txt",sep="\t",quote=F,row.names = F) #输出文件
tiff(file="barplot.tiff",width = 20,height = 20,units ="cm",compression="lzw",bg="white",res=600) #柱状图
barplot(kk, drop = TRUE, showCategory = 12) #柱状图,显示多少个通路
dev.off()
tiff(file="dotplot.tiff",width = 20,height = 20,units ="cm",compression="lzw",bg="white",res=600) #气泡图
dotplot(kk, showCategory = 12) #气泡图,显示多少个通路
dev.off()
library("pathview") #通路图,调用包,从网站下载通路
keggxls=read.table("KEGG.txt",sep="\t",header=T) #读取文件
for(i in keggxls$ID){
pv.out <- pathview(gene.data = sum, pathway.id = i, species = "hsa", out.suffix = "pathview") #数据是sum列,ID是i列,物种是人即has,输出为pathview,
}
遇到问题:下载GEO数据,Series Matrix File(s),下载只有4K

!Series_title “Circular RNAs are super abundant in cervical tumor and plasma detected by high throughput microarray [cervical_cancer_circRNA]”
!Series_geo_accession “GSE90737”
!Series_status “Public on Mar 01 2018”
!Series_submission_date “Dec 01 2016”
!Series_last_update_date “May 31 2018”
!Series_pubmed_id “29415187”
经过咨询大神,查阅资料,只能通过R来下载,如下:
eSet=getGEO('GSE102686')
运行后则数据下载到后台,进一步用于数据分析。
期间也用了国内镜像下载,但是可以下载示例数据,自己想要的数据不一定下载的到,可能是作者未进行更新。
请关注“恒诺新知”微信公众号,感谢“R语言“,”数据那些事儿“,”老俊俊的生信笔记“,”冷🈚️思“,“珞珈R”,“生信星球”的支持!