环状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来下载，如下：