这里是佳奥！新的篇章！

我们先从背景知识开始。

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.png

rm(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 第二个

QQ截图20220828211822.png

3 miRNA-Seq流程认知

• 数据质量控制

  • 单端，50nt足够，价格贵

• 比对到参考基因组

• 比对结果文件说明

• 已知 miRNA 表达谱构建

• 新miRNA预测

  • 主要是对未注释上任何RNA且比对上基因组外显子反义链、内含子、基因间区的sRNAs

  • 通过选用软件 Mireap(该软件适用于动植物)或mirdeep(该软件适用于动物)筛选miRNA的生物特征得到的

  • miRNA初始转录位点多位于基因间隔区、内含子以及编码序列的反向重复序列上

  • 其前体具有标志性的发夹结构，成熟体的形成是由Dicer酶的剪切实现的。

• 差异表达分析

• 靶基因预测

• 靶基因功能分析

参考推文：

https://mp.weixin.qq.com/s/fJgWNGH7SGi4H89d6prvPQ

步骤多种多样，目的是把测序样品很好的mapping到参考基因组上。

还有一些整合的流程，QuickMIRSeq、sRNAnalyzer、COMPASS等等。

下一篇我们开始正式的实战分析，从环境搭建开始。

我们下一篇再见！