11. 对第10步得到的表达矩阵进行探索,先画第一个样本的所有基因的表达量的boxplot,hist,density , 然后画所有样本的 这些图 参考: alization_for_expression_matrixl 理解ggplot2的绘图语法,数据和图形元素的映射关系  ggplot2 library(ggplot2) p=ggplot(exprset_l,aes(x=sample,y=val e,fill=gro p))+geom_boxplot print(p) p=ggplot(exprset_l,aes(x=sample,y=val e,fill=gro p))+geom_violin print(p) p=ggplot(exprset_l,aes(val e,fill=gro p))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4) print(p) p=ggplot(exprset_l,aes(val e,col=gro p))+geom_density+facet_wrap(~sample, nrow = 4) print(p) p=ggplot(exprset_l,aes(val e,col=gro p))+geom_density print(p) p=ggplot(exprset_l,aes(x=sample,y=val e,fill=gro p))+geom_boxplot p=p+stat_s mmary(f n.y="mean",geom="point",shape=23,size=3,fill="red") p=p+theme_set(theme_set(theme_bw(base_size=20 p=p+theme(text=element_text(face= bold ),axis.text.x=element_text(angle=30,hj st=1),axis.title=element_blank) print(p) image image image image image image 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) g_mad names(g_mad) [1] "dusp5" "igfbp4" "h1fx" "enpp2" "flna" "clec2b" "tspyl2" "znf266" "s100a9" "nr4a2" "tgfbi" "arf6" "apbb2" "vcan" "rbm38" [16] "capg" "plxnc1" "rgs2" "rnase6" "vamp5" "cybb" "gnly" "ccl3" "oas1" "enpp2" "trib2" "znf804a" "h1fx" "igh" "jund" [31] "slc25a1" "pcdh9" "vipr1" "cobll1" "gusbp11" "s100a8" "hbb" "fos" "lhfpl2" "fcn1" "zap70" "iglc1" "lgals1" "hbb" "fos" [46] "slamf1" "tcf7" "dmd" "igf2bp3" "fam30a" 13.根据第12步骤得到top 50 mad值的基因列表来取表达矩阵的子集,并且热图可视化子表达矩阵.试试看其它5种热图的包的不同效果. ## heatmap library(pheatmap) choose_gene=names(tail(sort(apply(exprset,1,mad)),50)) choose_matrix=exprset[choose_gene,] choose_matrix=t(scale(t(choose_matrix pheatmap(choose_matrix) image 14.取不同统计学指标mean,median,max,mean,sd,var,mad的各top50基因列表,使用upsetr包来看他们之间的overlap情况. ## upsetr #  library(upsetr) g_all - niq e(c(names(g_mean),names(g_median),names(g_max),names(g_min), names(g_sd),names(g_var),names(g_mad) )) 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) ) pset(dat,nsets = 7) image 15.在第二步的基础上面提取cll包里面的data(scllex) 数据对象的样本的表型数据. pdata=pdata(scllex) gro p_list=as.character(pdata[,2]) gro p_list # [1] "progres." "stable" "progres." "progres." "progres." "progres." "stable" "stable" "progres." "stable" "progres." "stable" "progres." "stable" # [15] "stable" "progres." "progres." "progres." "progres." "progres." "progres." "stable" dim(exprset) # [1] 8585 22 exprset[1:5,1:5] # cll11.cel cll12.cel cll13.cel cll14.cel cll15.cel mapk3 5. 6. 5. 5. 5. tie1 2. 2. 2. 2. 2. cyp2c19 3. 3. 3. 3. 3. cxcr5 7. 7. 7. 7. 6. dusp1 5.0 7. 7. 6. 5.0 16.对所有样本的表达矩阵进行聚类并且绘图,然后添加样本的临床表型数据信息(更改样本名) ## hcl st colnames(exprset)=paste(gro p_list,1:22,sep= ) # define nodepar nodepar - list(lab.cex = 0.6, pch = c(na, 19), cex = 0.7, col = "bl e") hc=hcl st(dist(t(exprset par(mar=c(5,5,5,10)) plot(as.dendrogram(hc), nodepar = nodepar, horiz = true) image 17.对所有样本的表达矩阵进行pca分析并且绘图,同样要添加表型信息. # install.packages("ggfortify") library(ggfortify) exprset - exprs(scllex) df - as.data.frame(t(exprset)) df$gro p - gro p_list # a toplot ses ggplot2 to draw a partic lar plot for an object of a partic lar class in a single command. a toplot(prcomp(df[,1:(ncol(df)-1)]), data=df, colo r = gro p ) image 18.根据表达矩阵及样本分组信息进行批量t检验,得到检验结果表格 ## t.test dat = exprset gro p_list=as.factor(gro p_list) gro p1 = which(gro p_list == levels(gro p_list)[1]) gro p2 = which(gro p_list == levels(gro p_list)[2]) dat1 = dat[, gro p1] dat2 = dat[, gro p2] dat = cbind(dat1, dat2) pvals = apply(exprset, 1, f nction(x){ t.test(as.n meric(x)~gro p_list)$p.val e }) p.adj = p.adj st(pvals, method = "bh") avg_1 = rowmeans(dat1) avg_2 = rowmeans(dat2) 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]),] deg_t.test=as.data.frame(deg_t.test) head(deg_t.test) # avg_1 avg_2 log2fc pvals p.adj _at 7. 8. 0. 1.e-05 0. _at 6. 7. 1. 4.0e-05 0. _at 7. 5. -1. 6.e-05 0. _at 4. 2. -2. 8.e-05 0. _at 5. 7.0 1.0 9.e-05 0. _at 4. 3. -0. 2.e-04 0. 19.使用limma包对表达矩阵及样本分组信息进行差异分析,得到差异分析表格,重点看logfc和p值,画个火山图(就是logfc和-log10(p值)的散点图). # deg by limma s ppressmessages(library(limma)) design - model.matrix(~0+factor(gro p_list)) colnames(design)=levels(factor(gro p_list)) rownames(design)=colnames(exprset) design contrast.matrix-makecontrasts(paste0( niq e(gro p_list),collapse = "-"),levels = design) contrast.matrix ##这个矩阵声明,我们要把progres.组跟stable进行差异分析比较 ##step1 fit - lmfit(exprset,design) ##step2 fit2 - contrasts.fit(fit, contrast.matrix) ##这一步很重要,大家可以自行看看效果 fit2 - ebayes(fit2) ## defa lt no trend !!! ##ebayes with trend=true ##step3 tempo tp t = toptable(fit2, coef=1, n=inf) nrdeg = na.omit(tempo tp t) #write.csv(nrdeg2,"limma_notrend.res lts.csv",q ote = f) head(nrdeg) ## volcano plot deg=nrdeg logfc_c toff - with(deg,mean(abs( logfc)) + 2*sd(abs( logfc)) ) deg$change = as.factor(ifelse(deg$p.val e  0.05 & abs(deg$logfc)  logfc_c toff, ifelse(deg$logfc  logfc_c toff , up , down ), not ) ) this_tile - paste0( c toff for logfc is ,ro nd(logfc_c toff,3), nthe n mber of p gene is ,nrow(deg[deg$change == up ,]) , nthe n mber of down gene is ,nrow(deg[deg$change == down ,]) ) g = ggplot(data=deg, aes(x=logfc, y=-log10(p.val e), 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-val e") + ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hj st = 0.5))+ scale_colo r_man al(val es = c( bl e , black , red )) ## corresponding to the levels(res$change) print(g) image image 20.对t检验结果的p值和limma包差异分析的p值画散点图,看看哪些基因相差很大.  different p val es head(nrdeg) head(deg_t.test) deg_t.test=deg_t.test[rownames(nrdeg),] plot(deg_t.test[,3],nrdeg[,1]) plot(deg_t.test[,4],nrdeg[,4]) plot(-log10(deg_t.test[,4]),-log10(nrdeg[,4])) image image image rownames(exprset)=ids[match(rownames(exprset),ids$probe_id),2] exprset[1:5,1:5] exprset[ gapdh ,] exprset[ actb ,] exprset[ dleu1 ,] library(ggplot2) library(ggp br) my_comparisons - list( c("stable", "progres.") ) dat=data.frame(gro p=gro p_list, sample names(exprset[ dleu1 ,]), val es= as.n meric(exprset[ dleu1 ,])) ggboxplot( dat, x = "gro p", y = "val es", color = "gro p", add = "jitter" )+ stat_compare_means(comparisons = my_comparisons, method = "t.test") image ## heatmap library(pheatmap) choose_gene=head(rownames(nrdeg),25) choose_matrix=exprset[choose_gene,] choose_matrix=t(scale(t(choose_matrix pheatmap(choose_matrix) image 生信技能树公益视频合辑:学习顺序是lin x,r,软件安装,geo,小技巧,ngs组学!