美文网首页做图R plotIMP research
吨吨吨的生信技能树r语言习题(1)

吨吨吨的生信技能树r语言习题(1)

作者: 吨吨吨_e0c9 | 来源:发表于2022-02-24 01:19 被阅读0次

    生信技能树高级20题

    1.安装一些R包
    数据包: ALL, CLL, pasilla, airway
    软件包:limma,DESeq2,clusterProfiler
    工具包:reshape2
    绘图包:ggplot2
    可以在bioconductor里面下载

    2. 了解ExpressionSet对象,比如CLL包里面就有data(sCLLex) ,找到它包含的元素,提取其表达矩阵(使用exprs函数),查看其大小
    用limma对芯片数据做差异分析:
    http://www.bio-info-trainee.com/bioconductor_China/software/limma.html
    ExpressionSet 对象简单讲解:
    https://github.com/bioconductor-china/basic/blob/master/ExpressionSet.md

    3.了解 str,head,help函数,作用于 第二步提取到的表达矩阵
    str(): show the structure of an Arbitrary R Object

    4.安装并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后 显示的那些变量
    正则表达式

    if (!require("BiocManager", quietly = TRUE)) 
      install.packages("BiocManager")
    

    5.理解 head(toTable(hgu95av2SYMBOL)) 的用法,找到 TP53 基因对应的探针ID
    hgu95av2SYMBOL中提供的是l两列probe_id和symbol名。

    ids=toTable(hgu95av2SYMBOL)
    ids[grep('TP53',ids$symbol),] #筛选内容和从哪筛选
    subset(ids,ids$symbol=='TP53')#从哪筛选然后筛选啥
    

    grep函数和subset函数类似,都可以用来筛选,dplyr包的filter也可以
    \color{green}{grep函数}R 函数笔记 | grep()函数与R语言中的正则表达式 - 简书 (jianshu.com)
    \color{green}{subset函数}(10条消息) R语言subset函数_Christina-CSDN博客_r语言subset函数

    6.理解探针与基因的对应关系,总共多少个基因,基因最多对应多少个探针,是哪些基因,是不是因为这些基因很长,所以在其上面设计多个探针呢?

    length(unique(ids$symbol)) #查看有多少基因
    dim(exprSet) #查看多少探针
    head(sort(table(ids$symbol),decreasing=T))
    

    7.第二步提取到的表达矩阵是12625个探针在22个样本的表达量矩阵,找到那些不在 hgu95av2.db 包收录的对应着SYMBOL的探针。

    b2 <- exprSet[!(rownames(exprSet) %in% ids$probe_id),] #!表示非
    

    8.过滤表达矩阵,删除那1165个没有对应基因名字的探针。

    exprSet= exprSet[(rownames(exprSet) %in% ids$probe_id),]
    

    9.✨✨✨整合表达矩阵,多个探针对应一个基因的情况下,只保留在所有样本里面平均表达量最大的那个探针
    提示,理解 tapply,by,aggregate,split 函数 , 首先对每个基因找到最大表达量的探针。
    然后根据得到探针去过滤原始表达矩阵

    法一:
    \color{green}{by函数}: by函数的功能是整理与整合大致理念为by(要整理的内容, 分类依据, summary等function)

       tmp = by(exprSet,ids$symbol, function(x) rownames(x)[which.max(rowMeans(x))] )
    #这个操作首先是根据symbol对exprSet进行分类;然后对同一类数据计算其行平均值;最后找出平均值最大的行,取其行名。这里即探针名
       probes = as.character(tmp) #重新命名
       dim(exprSet)
       exprSet=exprSet[rownames(exprSet) %in% probes ,]#整理exprSet表留下取得探针
       dim(exprSet)
    

    法二:我之前还问z老师这个是甚么,现在我只想为我的愚蠢自罚一杯

      identical(ids$probe_id,rownames(exprSet))
      dat=exprSet
      ids$median=apply(dat,1,median) #ids新建median这一列,列名为median,同时对dat这个矩阵按行操作,取每一行的中位数,将结果给到median这一列的每一行
      ids=ids[order(ids$symbol,ids$median,decreasing = T),]#对ids$symbol按照ids$median中位数从大到小排列的顺序排序,将对应的行赋值为一个新的ids
      ids=ids[!duplicated(ids$symbol),]#将symbol这一列取取出重复项,'!'为否,即取出不重复的项,去除重复的gene ,保留每个基因最大表达量结果s
      dat=dat[ids$probe_id,] #新的ids取出probe_id这一列,将dat按照取出的这一列中的每一行组成一个新的dat
      rownames(dat)=ids$symbol#把ids的symbol这一列中的每一行给dat作为dat的行名
      dat[1:4,1:4]  #保留每个基因ID第一次出现的信息
      dim(dat)
    

    10.把过滤后的表达矩阵更改行名为基因的symbol,因为这个时候探针和基因是一对一关系了。

    这里针对法一继续做了下去
    rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
    

    11.对第10步得到的表达矩阵进行探索,先画第一个样本的所有基因的表达量的boxplot,hist,density , 然后画所有样本的这些图
    我只画了所有的

    library(reshape2)
    exprSet_L=melt(exprSet)
    colnames(exprSet_L)=c('probe','sample','value')
    exprSet_L$group=rep(group_list,each=nrow(exprSet))
    p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
    print(p)
    p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
    print(p)
    p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
    print(p)
    



    12. 理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标,最后按照mad值排序,取top 50 mad值的基因,得到列表。

    g_mean <- tail(sort(apply(exprSet,1,mean)),50)
    g_median <- tail(sort(apply(exprSet,1,median)),50)
    g_max <- tail(sort(apply(exprSet,1,max)),50)
    g_min <- tail(sort(apply(exprSet,1,min)),50)
    g_sd <- tail(sort(apply(exprSet,1,sd)),50)
    g_var <- tail(sort(apply(exprSet,1,var)),50)
    g_mad <- tail(sort(apply(exprSet,1,mad)),50)
    

    MAD(Median absolute deviation, 中位数绝对偏差)是单变量数据集中样本差异性的稳健度量。mad是一个健壮的统计量,对于数据集中异常值的处理比标准差更具有弹性,可以大大减少异常值对于数据集的影响。
    对于单变量数据集 X={X1,X2,X3,...,Xn}X={X1,X2,X3,...,Xn},mad的计算公式为:
    mad(X)=median(|Xi−median(X)|)

    13.根据第12步骤得到top 50 mad值的基因列表来取表达矩阵的子集,并且热图可视化子表达矩阵。试试看其它5种热图的包的不同效果。
    5种包分别是heatmap函数、ggplot2包、gplot包、lattice包
    (没有lattice包,朕倦了)
    (1)heatmap包

    choose_gene=names(tail(sort(apply(exprSet,1,mad)),50))
    choose_matrix=exprSet[choose_gene,]
    choose_matrix=t(scale(t(choose_matrix)))
    heatmap(choose_matrix)
    

    (2)pheatmap包

    library(pheatmap)
    annotation_col =data.frame(group_list)
    row.names(annotation_col) <- colnames(choose_matrix)
    pheatmap(choose_matrix,
             color=colorRampPalette(c("navy","white","firebrick3"))(100),
             border_color = NA,
             annotation_col = annotation_col)
    

    (3)gpolots

    factor=as.factor(group_list)
    plot_color = c('orange','green')[factor]
    gplots::heatmap.2(choose_matrix,ColSideColors = plot_color)
    

    (4)ggplots2

    library(ggplot2)
    library(reshape2)
    c <- as.data.frame(melt(choose_matrix))
    p1<-ggplot(c,aes(x=Var2 ,y=Var1)) #热图绘制
    p2 <- p1+scale_color_gradientn(values = seq(0,1,0.2),colours = c('#6699CC','#FFFF99','#CC3333'))+
      theme_bw()+
      geom_point(aes(size=`value`,
                     color=`value`))+
      theme(panel.grid = element_blank(),axis.text.x =element_text(angle =90,hjust =0.5,vjust = 0.5))+
      xlab('sample') + ylab('gene')+
      theme(panel.border = element_rect(fill=NA,color="black", size=1, linetype="solid"))+
      geom_vline(xintercept=c(-3,3.5,8.5,11.5,15.5),size=.8)
    p2
    

    14.取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表,使用UpSetR包来看他们之间的overlap情况。(一个结合可视化包)

    对UpsetR的一个介绍:
    UpSetR:集合可视化神包 - 简书 (jianshu.com)

    library(UpSetR)
    g_all <- unique(c(names(g_mean),names(g_median),names(g_max),names(g_min),
                      names(g_sd),names(g_var),names(g_mad) ))#看所有前50里面所有的不重复的基因有哪些等于取交集
    dat=data.frame(g_all=g_all,
                   g_mean=ifelse(g_all %in%  names(g_mean) ,1,0),
                   g_median=ifelse(g_all %in%  names(g_median) ,1,0),
                   g_max=ifelse(g_all %in%  names(g_max) ,1,0),
                   g_min=ifelse(g_all %in%  names(g_min) ,1,0),
                   g_sd=ifelse(g_all %in%  names(g_sd) ,1,0),
                   g_var=ifelse(g_all %in%  names(g_var) ,1,0),
                   g_mad=ifelse(g_all %in%  names(g_mad) ,1,0))#看交集在不在各个集合里面在的1,不在0
    upset(dat,nsets = 7)#nset参数显示最多展示多少个集合数据
    

    上图



    如何看图:
    黑色表示该位置有数据,灰色的点表示没有。不同点连线表示存在交集。具体数据可以看上面的条形图。不同类型的数据的总量看左边的条形图。所以本图说明了,首先mean,median,max,mean,sd,var,mad列表里都有50个基因(这好像是废话),所以左面的条长度相同。可以看到同时拥有交集最多的是mean,median,max,min,即有36个基因他的这些指标都在前50里面。其余的同。这个表还可以各种改变得炫酷。
    这个方法的目的应该是更方便大家看到哪些基因表达的比较厉害,表达的厉害可能和这个病就有点关系了(猜的)

    15.在第二步的基础上面提取CLL包里面的data(sCLLex) 数据对象的样本的表型数据。

    data("sCLLex")
    pdata=pData(sCLLex) #这个其实一开始就提取了
    

    这里面是sample ID 和disease

    16.对所有样本的表达矩阵进行聚类并且绘图,然后添加样本的临床表型数据信息(更改样本名)

    pdata=pData(sCLLex)
    colnames(exprSet)=paste(group_list,1:22,sep='')
    #这一步我本来想用pData操作建一个新的列但是失败了,因为他是一个function而不是一个dataframe
    cluster<- t(exprSet)
    hc <- hclust(dist(cluster))
    plot(as.dendrogram(hc))
    #jimmy老师并不是这样聚类的
    nodePar <- list(lab.cex = 0.6, pch = c(NA, 19), 
                    cex = 0.7, col = "blue")
    hc=hclust(dist(t(exprSet)))
    par(mar=c(5,5,5,10)) 
    plot(as.dendrogram(hc), nodePar = nodePar,  horiz = TRUE)
    #嗯区别不大
    
    普通层次聚类
    炫酷聚类

    17.对所有样本的表达矩阵进行PCA分析并且绘图,同样要添加表型信息。

    pc <- prcomp(t(exprSet),scale=TRUE)
    pcx=data.frame(pc$x)
    pcr=cbind(samples=rownames(pcx),group_list, pcx) 
    library(psych)
    library(ggplot2)
    p  <- ggplot(pcr, aes(x=PC1, y=PC2,color=pcr$group_list )) +
      geom_point(size=4,alpha=0.5)+geom_text(aes(label=samples),hjust=-0.1, vjust=-0.3)+
      theme(axis.text= element_text(size=20))+
      theme(panel.grid.major =element_blank(),
            panel.grid.minor = element_line(color="steelblue"),
            panel.background = element_blank(),
            axis.line = element_line(colour = "black"))+
      stat_ellipse(lwd=1,level = 0.8) 
    ##再贴上jimmy老师的代码
    法一
    # BiocManager::install('ggfortify')
    library(ggfortify)
    df=as.data.frame(t(exprSet))
    df$group=group_list 
    autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')
    法二
    library("FactoMineR")#画主成分分析图需要加载这两个包
    library("factoextra") 
    df=as.data.frame(t(exprSet))
    dat.pca <- PCA(df, graph = FALSE)#现在dat最后一列是group_list,需要重新赋值给一个dat.pca,这个矩阵是不含有分组信息的
    fviz_pca_ind(dat.pca,
                 geom.ind = "point", # show points only (nbut not "text")
                 col.ind = group_list, # color by groups
                 # palette = c("#00AFBB", "#E7B800"),
                 addEllipses = TRUE, # Concentration ellipses
                 legend.title = "Groups"
    )
    
    
    我画的pca
    jimmy老师法一
    jimmy老师法二

    有没有霸霸解释一下为什么这两种方法一二轴的解释量不一样啊😂😂

    18.据表达矩阵及样本分组信息进行批量T检验,得到检验结果表格

    dat = exprSet
    group_list=as.factor(group_list)  ##将grouplist化为factor
    group1 = which(group_list == levels(group_list)[1])##progress组
    group2 = which(group_list == levels(group_list)[2])##stable组
    dat1 = dat[, group1]##progress的基因表达
    dat2 = dat[, group2]#stable的基因表达
    dat = cbind(dat1, dat2)#二合一
    pvals = apply(exprSet, 1, function(x){
      t.test(as.numeric(x)~group_list)$p.value
    })#每一行progress数据和stable数据产生的p值
    p.adj = p.adjust(pvals, method = "BH")#bh调整p值
    avg_1 = rowMeans(dat1)#progress平均数
    avg_2 = rowMeans(dat2)#stable平均数
    log2FC = avg_2-avg_1#二者差
    DEG_t.test = cbind(avg_1, avg_2, log2FC, pvals, p.adj)
    DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),]#将pvalue一列从小到大排序
    DEG_t.test=as.data.frame(DEG_t.test)
    head(DEG_t.test)
    

    19.使用limma包对表达矩阵及样本分组信息进行差异分析,得到差异分析表格,重点看logFC和P值,画个火山图(就是logFC和-log10(P值)的散点图。)。
    开始的时候已经学习了limma包
    其中有很多东西我没明白,先保留一些资料
    对于R语言中的model.matrix的一些简单理解 - R - TONX (treeh.cn)
    model matrix - 简书 (jianshu.com)
    \color{green}{limma包 }
    那些常用的limma操作 - 简书 (jianshu.com)

    suppressMessages(library(limma)) 
    design <- model.matrix(~0+factor(group_list))
    colnames(design)=levels(factor(group_list))
    rownames(design)=colnames(exprSet)
    design#制作分组矩阵,将progress和stable用0,1表示
    
    ## 下面的 contrast.matrix 矩阵非常重要,制定了谁比谁这个规则
    contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
    ##这个矩阵声明,我们要把progres.组跟stable进行差异分析比较
    ##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)
    ##火山图
    DEG=nrDEG
    logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )#设置cutoff
    DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                                  ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
    )
    this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                        '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                        '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
    )
    this_tile
    head(DEG)
    g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
      geom_point(alpha=0.4, size=1.75) +
      theme_set(theme_set(theme_bw(base_size=20)))+
      xlab("log2 fold change") + ylab("-log10 p-value") +
      ggtitle( this_tile  ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
      scale_colour_manual(values = c('blue','black','red'))  ## corresponding to the levels(res$change)
    print(g)
    
    

    20.对T检验结果的P值和limma包差异分析的P值画散点图,看看哪些基因相差很大。

    DEG_t.test=DEG_t.test[rownames(nrDEG),]
    plot(DEG_t.test[,3],nrDEG[,1]) ## 对比logFC,相反
    plot(DEG_t.test[,4],nrDEG[,4]) # 对比p,可以看到使用limma包和t.test本身的p值差异尚可接受
    plot(-log10(DEG_t.test[,4]),-log10(nrDEG[,4]))
    new=nrDEG[DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,]
    choose_gene=rownames(head(new[order(-new$abslog),],25))
    choose_matrix=exprSet[choose_gene,]
    choose_new=nrDEG[DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,]
    matrix=t(scale(t(choose_matrix))) 
    pheatmap(choose_matrix)
    
    logFC
    p value

    取P.Value < 0.05 & abs(logFC) > logFC_cutoff的基因按照abslogFC排序作图,这些基因差异比较大

    相关文章

      网友评论

        本文标题:吨吨吨的生信技能树r语言习题(1)

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