美文网首页转录组生信
【miRNA-Seq 实战】一、数据库探索,miRNA芯片数据分

【miRNA-Seq 实战】一、数据库探索,miRNA芯片数据分

作者: 佳奥 | 来源:发表于2022-08-28 21:41 被阅读0次

    这里是佳奥!新的篇章!

    我们先从背景知识开始。

    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到参考基因组上。

    还有一些整合的流程,QuickMIRSeqsRNAnalyzerCOMPASS等等。

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

    我们下一篇再见!

    相关文章

      网友评论

        本文标题:【miRNA-Seq 实战】一、数据库探索,miRNA芯片数据分

        本文链接:https://www.haomeiwen.com/subject/rticnrtx.html