这里是佳奥!新的篇章!
我们先从背景知识开始。
1 miRBase探索
miRBase数据库收录miRNA序列:
mirbase.org
在download目录可下载:
QQ截图20220828141102.png
我们在Linux下查看前体序列hairpin.fa(解压缩$ gzip -d hairpin.fa.gz
):mir表示
$ cat hairpin.fa | head
>cel-let-7 MI0000001 Caenorhabditis elegans let-7 stem-loop
UACACUGUGGAUCCGGUGAGGUAGUAGGUUGUAUAGUUUGGAAUAUUACCACCGGUGAAC
UAUGCAAUUUUCUACCUUACCGGAGACAGAACUCUUCGA
>cel-mir-1 MI0000003 Caenorhabditis elegans miR-1 stem-loop
AAAGUGACCGUACCGAGCUGCAUACUUCCUUACAUGCCCAUACUAUAUCAUAAAUGGAUA
UGGAAUGUAAAGAAGUAUGUAGAACGGGGUGGUAGU
可以看到这里面碱基是AUCG而不是ATCG。
而成熟体序列mature.fa:miR表示
$ cat mature.fa | head
>cel-let-7-5p MIMAT0000001 Caenorhabditis elegans let-7-5p
UGAGGUAGUAGGUUGUAUAGUU
>cel-miR-1-5p MIMAT0020301 Caenorhabditis elegans miR-1-5p
CAUACUUCCUUACAUGCCCAUA
##查看物种以及数量
$ grep '>' hairpin.fa | awk '{print $3,$4}' | sort | uniq -c
$ grep '>' hairpin.fa | awk '{print $3,$4}' | sort | uniq -c | wc
265 795 7120
##265个物种
##查看人类,前体有1917条
$ grep '>' hairpin.fa | awk '{print $3,$4}' | sort | uniq -c | grep 'Homo'
1917 Homo sapiens
##人类成熟体有2656条
$ grep '>' mature.fa | awk '{print $3,$4}' | sort | uniq -c | grep 'Homo'
2656 Homo sapiens
前体miRNA和成熟的miRNA的区别,前体miRNA长度一般是70–120 碱基,一般是茎环结果,也就是发夹结构,所以叫做hairpin。成熟之后,一般22 个碱基。
查看成熟体的序列长度:确保当前环境有perl
$ grep -v '>' mature.fa | perl -alne '{print length}' | head
22
22
21
22
22
21
22
23
22
22
##查看长度分布
$ grep -v '>' mature.fa | perl -alne '{print length}' | sort | uniq -c
3 15
16 16
113 17
430 18
941 19
3374 20
13652 21##
18301 22##最多分布
7694 23##
3238 24
819 25
177 26
70 27
33 28
18 29
3 30
1 32
1 33
1 34
单独查看人类的情况:生成人类的.human.fa
perl -alne '{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};next if $tmp!=1;s/U/T/g if !/>/;print }' hairpin.fa > hairpin.human.fa
perl -alne '{if(/^>/){if(/Homo/){$tmp=1}else{$tmp=0}};next if $tmp!=1;s/U/T/g if !/>/;print }' mature.fa > mature.human.fa
##查看人类成熟体序列长度分布
$ grep -v '>' mature.human.fa | perl -alne '{print length}' | sort | uniq -c
8 16
58 17
76 18
98 19
171 20
530 21
1177 22##最多
384 23
103 24
37 25
10 26
3 27
1 28
生成的hairpin.human.fa/mature.human.fa
称为参考基因集,后续的测序文件要和他们比较。
还有一些网页工具,做序列比对以及预测靶向基因等,这些等到后面的实战再演示。
2 miRNA芯片数据分析流程
miRNA芯片表达矩阵
三个公司:affymetrix、agilent、illumina。
本次实战重复的文章:miRNA expression data in breast cancer and adjacent normal tissues
https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE143564
分析:差异分析、网络分析、生存分析。
我们开始下载数据:GSE143564
在R中操作,有关代码参考:step0-step4
https://github.com/jmzeng1314/GEO/tree/master/GSE42872_main
##因为这个GEO数据集只有一个GPL平台,所以下载到的是一个含有一个元素的list
rm(list = ls())
options(stringsAsFactors = F)
library(AnnoProbe)
library(GEOquery)
gset <- GEOquery::getGEO('GSE143564',
AnnotGPL = F,
getGPL = F)
a <- gset[[1]]
dat <- exprs(a)#a现在是一个对象,取a这个对象通过看说明书知道要用exprs这个函数
dim(dat)#看一下dat这个矩阵的维度
##查看下载
> a
Annotation: GPL21572
##归一化数据
dat[1:4,1:4]#查看dat这个矩阵的1至4行和1至4列,逗号前为行,逗号后为列
boxplot(dat,las=2)##生成未归一化表达量箱线图
library(limma)
dat <- normalizeBetweenArrays(dat)
boxplot(dat,las=2)##生成归一化表达量箱线图
pd <- pData(a)##生成分组信息,三个Normal三个Tumor。通过查看说明书知道取对象a里的临床信息用pData
##挑选一些感兴趣的临床表型。
未归一化表达量:
QQ截图20220828170100.png
归一化后:
QQ截图20220828170216.png
##拿到探针和mRNA信息
library(stringr)
group_list <- rep(c('N','T'),each=3)
table(group_list)
dat[1:4,1:4]
a
##下载注释文件,结尾会用
options('download.file.method.GEOquery'='libcurl')
library(GEOquery)
gpl <- getGEO('GPL21572')
##把结果保存一下
save(dat,file='step1-output.R')
差异分析前——通过热图查看分组情况
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'step1-output.R')
group_list <- rep(c('N','T'),each=3)
table(group_list)
##每次都要检测数据
dat[1:4,1:4]
##下面是画PCA的必须操作
dat=t(dat)#画PCA图时要求是行名时样本名,列名时探针名,因此此时需要转换
dat=as.data.frame(dat)#将matrix转换为data.frame
dat=cbind(dat,group_list)#cbind横向追加,即将分组信息追加到最后一列
library("FactoMineR")#画主成分分析图需要加载这两个包
library("factoextra")
#The variable group_list (index = 54676) is removed
#before PCA analysis
dat.pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#现在dat最后一列是group_list,需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
fviz_pca_ind(dat.pca,
geom.ind = "point", # show points only (nbut not "text")
col.ind = dat$group_list, # color by groups
# palette = c("#00AFBB", "#E7B800"),
addEllipses = TRUE, # Concentration ellipses
legend.title = "Groups"
)
ggsave('all_samples_PCA.png')
QQ截图20220828190605.png
热图
rm(list = ls())
load(file = 'step1-output.R')
dat[1:4,1:4]
cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行('1'是按行取,'2'是按列取)取每一行的方差,从小到大排序,取最大的1000个
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F)#对那些提取出来的1000个基因所在的每一行取出,组合起来为一个新的表达矩阵
n=t(scale(t(dat[cg,]))) #'scale'可以对log-ratio数值进行归一化
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
group_list <- rep(c('N','T'),each=3)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #把ac的行名给到n的列名,即对每一个探针标记上分组信息
pheatmap(n,show_colnames =F,show_rownames = F,
annotation_col=ac,filename = 'heatmap_top1000_sd.png')
dev.off()/dev.new()##不出图了加这个
蓝色N红色T
heatmap_top1000_sd.pngrm(list = ls())
load(file = 'step1-output.R')
dat[1:4,1:4]
exprSet <- dat
pheatmap::pheatmap(cor(exprSet))
group_list <- rep(c('N','T'),each=3)
coID <- data.frame(group=group_list)
rownames(coID) <- colnames(exprSet)
pheatmap::pheatmap(cor(exprSet),
annotation_col = coID,
show_rownames = F)
##filename='cor_all.png'
QQ截图20220828192338.png
##取变异前500
exprSet <- exprSet[names(sort(apply(exprSet,1,mad),decreasing = T)[1:500]),]
dim(exprSet)
M<-cor(exprSet)
pheatmap::pheatmap(M,annotation_col = coID)
pheatmap::pheatmap(M,
show_rownames = F,
annotation_col = coID)
##filename='cor_top500.png'
QQ截图20220828194000.png
差异分析
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'step1-output.R')
# 每次都要检测数据
dat[1:4,1:4]
group_list <- rep(c('N','T'),each=3)
table(group_list) #table函数,查看group_list中的分组个数
#通过为每个数据集绘制箱形图,比较数据集中的数据分布
boxplot(dat[1,]~group_list) #按照group_list分组画箱线图
#定义一个函数g,函数为{}里的内容
bp=function(g){
library(ggpubr)
df=data.frame(gene=g,stage=group_list)
p <- ggboxplot(df, x = "stage", y = "gene",
color = "stage", palette = "jco",
add = "jitter")
#Add p-value
p + stat_compare_means()
}
bp(dat[1,]) ##调用上面定义好的函数,避免同样的绘图代码重复多次敲。
bp(dat[2,])
bp(dat[3,])
bp(dat[4,])
dim(dat)
library(limma)
design=model.matrix(~factor( group_list ))
fit=lmFit(dat,design)
fit=eBayes(fit)
##上面是limma包用法的一种方式
options(digits = 4) #设置全局的数字有效位数为4
#topTable(fit,coef=2,adjust='BH')
topTable(fit,coef=2,adjust='BH')
## 但是上面的用法做不到随心所欲的指定任意两组进行比较
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
head(design)
exprSet=dat
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts("T-N",
levels = design)
contrast.matrix ##这个矩阵声明,我们要把 Tumor 组跟 Normal 进行差异分析比较
deg = function(exprSet,design,contrast.matrix){
##step1
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix)
##这一步很重要,大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
return(nrDEG)
}
deg = deg(exprSet,design,contrast.matrix)
head(deg)
bp(dat[rownames(deg)[1],])
save(deg,file = 'deg.Rdata')
可以看到差异十分显著。
QQ截图20220828202042.png
绘制火山图
## for volcano
if(T){
nrDEG=deg
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(P.Value))
library(ggpubr)
df=nrDEG
df$v= -log10(P.Value) #df新增加一列'v',值为-log10(P.Value)
ggscatter(df, x = "logFC", y = "v",size=0.5)
df$g=ifelse(df$P.Value>0.01,'stable', #if 判断:如果这一基因的P.Value>0.01,则为stable基因
ifelse( df$logFC >2,'up', #接上句else 否则:接下来开始判断那些P.Value<0.01的基因,再if 判断:如果logFC >1.5,则为up(上调)基因
ifelse( df$logFC < -2,'down','stable') )#接上句else 否则:接下来开始判断那些logFC <1.5 的基因,再if 判断:如果logFC <1.5,则为down(下调)基因,否则为stable基因
)
table(df$g)
df$name=rownames(df)
head(df)
ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
label = "name", repel = T,
#label.select = rownames(df)[df$g != 'stable'] ,
label.select = head(rownames(deg)), #挑选一些基因在图中显示出来
palette = c("#00AFBB", "#E7B800", "#FC4E07") )
ggsave('volcano.png')
ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)
df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
table(df$p_c )
ggscatter(df,x = "AveExpr", y = "logFC", color = "p_c",size=0.2,
palette = c("green", "red", "black") )
ggsave('MA.png')
}
QQ截图20220828202446.png
QQ截图20220828202507.png
上下游较少,因为这里阈值选择是2,比较大。
绘制热图
## for volcano
if(T){
nrDEG=deg
head(nrDEG)
attach(nrDEG)
plot(logFC,-log10(P.Value))
library(ggpubr)
df=nrDEG
df$v= -log10(P.Value) #df新增加一列'v',值为-log10(P.Value)
ggscatter(df, x = "logFC", y = "v",size=0.5)
df$g=ifelse(df$P.Value>0.01,'stable', #if 判断:如果这一基因的P.Value>0.01,则为stable基因
ifelse( df$logFC >2,'up', #接上句else 否则:接下来开始判断那些P.Value<0.01的基因,再if 判断:如果logFC >1.5,则为up(上调)基因
ifelse( df$logFC < -2,'down','stable') )#接上句else 否则:接下来开始判断那些logFC <1.5 的基因,再if 判断:如果logFC <1.5,则为down(下调)基因,否则为stable基因
)
table(df$g)
df$name=rownames(df)
head(df)
ggscatter(df, x = "logFC", y = "v",size=0.5,color = 'g')
ggscatter(df, x = "logFC", y = "v", color = "g",size = 0.5,
label = "name", repel = T,
#label.select = rownames(df)[df$g != 'stable'] ,
label.select = head(rownames(deg)), #挑选一些基因在图中显示出来
palette = c("#00AFBB", "#E7B800", "#FC4E07") )
ggsave('volcano.png')
ggscatter(df, x = "AveExpr", y = "logFC",size = 0.2)
df$p_c = ifelse(df$P.Value<0.001,'p<0.001',
ifelse(df$P.Value<0.01,'0.001<p<0.01','p>0.01'))
table(df$p_c )
ggscatter(df,x = "AveExpr", y = "logFC", color = "p_c",size=0.2,
palette = c("green", "red", "black") )
ggsave('MA.png')
}
## for heatmap
if(T){
load(file = 'step1-output.R')
# 每次都要检测数据
dat[1:4,1:4]
group_list <- rep(c('N','T'),each=3)
table(group_list)
x=deg$logFC #deg取logFC这列并将其重新赋值给x
names(x)=rownames(deg) #deg取probe_id这列,并将其作为名字给x
cg=c(names(head(sort(x),100)),#对x进行从小到大排列,取前100及后100,并取其对应的探针名,作为向量赋值给cg
names(tail(sort(x),100)))
library(pheatmap)
pheatmap(dat[cg,],show_colnames =F,show_rownames = F) #对dat按照cg取行,所得到的矩阵来画热图
n=t(scale(t(dat[cg,])))#通过“scale”对log-ratio数值进行归一化,现在的dat是行名为探针,列名为样本名,由于scale这个函数应用在不同组数据间存在差异时,需要行名为样本,因此需要用t(dat[cg,])来转换,最后再转换回来
n[n>2]=2
n[n< -2]= -2
n[1:4,1:4]
pheatmap(n,show_colnames =F,show_rownames = F)
ac=data.frame(g=group_list)
rownames(ac)=colnames(n) #将ac的行名也就分组信息 给到n的列名,即热图中位于上方的分组信息
pheatmap(n,show_colnames =F,
show_rownames = F,
cluster_cols = F,
annotation_col=ac) #列名注释信息为ac即分组信息
##filename = 'heatmap_top200_DEG.png'
}
write.csv(deg,file = 'deg.csv')
dev.new()
QQ截图20220828203022.png
接下来,我们使用文章阈值重新作图一次。
我们需要先前的注释文件:
options('download.file.method.GEOquery'='libcurl')
library(GEOquery)
gpl <- getGEO('GPL21572')
##保存一下
save(gpl,file='GPLdownload.R')
load(file = 'GPLdownload.R')
QQ截图20220828211113.png
gpl真大。
##查看一下对象
> class(gpl)
[1] "GPL"
attr(,"package")
[1] "GEOquery"
##注释信息
tmp=gpl@dataTable@table
文章中的显著基因:miR-15b-5p 第二个
分析先到这里,后面还有其他芯片的分析,不一一列举。
3 miRNA-Seq流程认知
-
数据质量控制
- 单端,50nt足够,价格贵
-
比对到参考基因组
-
比对结果文件说明
-
已知 miRNA 表达谱构建
-
新miRNA预测
-
主要是对未注释上任何RNA且比对上基因组外显子反义链、内含子、基因间区的sRNAs
-
通过选用软件 Mireap(该软件适用于动植物)或mirdeep(该软件适用于动物)筛选miRNA的生物特征得到的
-
miRNA初始转录位点多位于基因间隔区、内含子以及编码序列的反向重复序列上
-
其前体具有标志性的发夹结构,成熟体的形成是由Dicer酶的剪切实现的。
-
-
差异表达分析
-
靶基因预测
-
靶基因功能分析
参考推文:
https://mp.weixin.qq.com/s/fJgWNGH7SGi4H89d6prvPQ
步骤多种多样,目的是把测序样品很好的mapping到参考基因组上。
还有一些整合的流程,QuickMIRSeq、sRNAnalyzer、COMPASS等等。
下一篇我们开始正式的实战分析,从环境搭建开始。
我们下一篇再见!
网友评论