实际上写完了这个全网最好的差异分析代码:免费的数据分析付费的成品代码 我就可以收工用来,但是永远不能低估粉丝的疑惑数量,任何一个细节都会被拿出来剖析。
比如代码里面我挑选了top1000的sd基因绘制热图,然后就可以分辨出来自己处理的数据集里面的样本分组是否合理啦。其实这个热图差不多等价于PCA分析的图,被我称为表达矩阵下游分析标准3图!详见:你确定你的差异基因找对了吗? ,就是下面的3张图: ![](http://image109.360doc.com/DownloadImg/2021/07/1412/226251755_1_20210714123854129_wm)
左边的热图,说明我们实验的两个分组,normal和npc的很多基因表达量是有明显差异的 中间的PCA图,说明我们的normal和npc两个分组非常明显的差异 右边的层次聚类也是如此,说明我们的normal和npc两个分组非常明显的差异
PS:如果你的转录组实验分析报告没有这三张图,就把我们生信技能树的这篇教程甩在他脸上,让他瞧瞧,学习下转录组数据分析。 为什么挑选top1000的sd基因绘制热图我这个热图是为了说明本分组是否合理,就是看样本的距离,这个时候你如果需要理解距离,那么你需要学习非常多细节知识。不仅仅是一个函数那么简单: r 语言中使用 dist ( x, method = “ euclidean ”, diag = FALSE, upper = FALSE, p = 2 ) 来计算距离。其中x是样本矩阵或者数据框。method 表示计算哪种距离。method 的取值有: r 语言中使用 hclust (d, method = “complete”, members=NULL) 来进行层次聚类。其中 d 为距离矩阵。method 表示类的合并方法,有: single 最短距离法 complete 最长距离法 median 中间距离法 mcquitty 相似法 average 类平均法 centroid 重心法 ward 离差平方和法
r 语言中主要使用 kmeans(x, centers, iter.max = 10, nstart = 1,algorithm =c(“Hartigan-Wong”, “Lloyd”,”Forgy”, “MacQueen”))来进行聚类。
也就是说,看起来非常简单的3张图,背后是几十年的统计学知识的基础建设。 当然了,也不要气馁哦,反正你只需要会看图就好!再次强调:你确定你的差异基因找对了吗? 里面的3张图: ![](http://image109.360doc.com/DownloadImg/2021/07/1412/226251755_1_20210714123854129_wm)
左边的热图,说明我们实验的两个分组,normal和npc的很多基因表达量是有明显差异的 中间的PCA图,说明我们的normal和npc两个分组非常明显的差异 右边的层次聚类也是如此,说明我们的normal和npc两个分组非常明显的差异
为什么选择top1000的sd基因绘制热图其实就是个人爱好,你可以探索top500,1000,2000,5000是否有区别。 比较不同的top基因聚类的差异同样的是一个表达矩阵和分组,如下: > dat[1:4,1:4] GSM312896 GSM312897 GSM312898 GSM312899 ZZZ3 8.286337 8.427789 9.315658 8.601464 ZZEF1 8.097163 8.419503 8.843223 8.963619 ZYX 8.187251 7.919804 8.082291 7.094361 ZYG11B 8.718149 8.195136 8.478548 9.054829 > dim(dat) [1] 18938 41 > table(group_list) group_list normal npc 10 31
然后检查不同top基因集的层次聚类情况 library(pheatmap) cg=names(tail(sort(apply(dat,1,sd)),1000)) p1=pheatmap(dat[cg,]) cg=names(tail(sort(apply(dat,1,sd)),500)) p2=pheatmap(dat[cg,]) cg=names(tail(sort(apply(dat,1,sd)),2000)) p3=pheatmap(dat[cg,]) cg=names(tail(sort(apply(dat,1,sd)),5000)) p4=pheatmap(dat[cg,]) tmp=data.frame(top1000=cutree(p1$tree_col,2), top500=cutree(p2$tree_col,2), top2000=cutree(p3$tree_col,2), top5000= cutree(p4$tree_col,2), group_list=group_list)
这个时候,你会发现,好像不一样,我修改层次聚类的类别数量 tmp=data.frame(top1000=cutree(p1$tree_col,4), top500=cutree(p2$tree_col,4), top2000=cutree(p3$tree_col,4), top5000= cutree(p4$tree_col,4), group_list=group_list)
然后你会发现,不管是top多少个基因,如果分成4类,那么1,2类都是normal,3,4类都是NPC,所以效果都挺好的。 ![image-20200205092709995](http://image109.360doc.com/DownloadImg/2021/07/1412/226251755_2_20210714123854301_wm) 其实基础知识点都在生信技能树了
|