Day5-小龙

4.检查数据质量

# 清空环境 
rm(list = objects( all = TRUE )) 
if (!is.null( dev.list() )) 
dev.off()
 
if (!dir.exists('./fig/')) {  
                                  dir.create("./fig/") 
}
 
load(file = './data/final_AssayData.Rdata')#加载之前保存的结果
 
#画PCA图检查数据
pca_data <- as.data.frame(t(AssayData)) 
library("FactoMineR") 
library("factoextra") 
dat.pca <- PCA(pca_data, graph = FALSE) 
get_eig(dat.pca) 
fviz_screeplot(dat.pca,
 addlabels = TRUE,
 ylim = c(0, 50)) 

## Contributions of variables to PC1 
fviz_contrib(dat.pca, choice = "var", axes = 1, top = 10) 
fviz_pca_ind( dat.pca,
              geom.ind = "point",
              col.ind = group_list,
palette = c("#00AFBB", "#FC4E07"),
              addEllipses = TRUE,              
legend.title = "Groups" )
 
library("ggfortify") 
pca_data$group = group_list 
autoplot(   prcomp( pca_data[, 1:(ncol( pca_data ) - 1)] ),
   data = pca_data,
   colour = 'group') + theme_bw()

PCA图

5.差异分析与可视化

library( "limma" ) #著名的limma包，做差异分析使用，还有其他的，这一段不明白怎么写出来的，我也在看说明书。
{
  design <- model.matrix( ~0 + factor( group_list ) )  
  colnames( design ) = levels( factor( group_list ) )  
  rownames( design ) = colnames( AssayData ) 
}
design#检查一下
 
contrast.matrix <- makeContrasts( "TNBC-control", levels = design ) #用乳腺癌的数据比上正常对照数据，一般都是用癌症：正常组织
contrast.matrix
 
{  fit <- lmFit( AssayData, design )
  fit2 <- contrasts.fit( fit, contrast.matrix )
   fit2 <- eBayes( fit2 )
  nrDEG <-  topTable( fit2, coef = 1, n = Inf )  save(nrDEG, file = './data/nrDEG.Rdata') 
} 
head(nrDEG)
 
# 差异分析到此结束，取上调和下调前50的基因画一个热图看看效果 
library( "pheatmap" ) 
{  
nrDEG_Z = nrDEG[ order( nrDEG$logFC ), ]  nrDEG_F = nrDEG[ order( -nrDEG$logFC ), ]  
choose_gene = c( rownames( nrDEG_Z )[1:50], rownames( nrDEG_F )[1:50] )  
choose_matrix = AssayData[ choose_gene, ]  choose_matrix = t( scale( t( choose_matrix ) ) )    choose_matrix[choose_matrix > 2] = 2  choose_matrix[choose_matrix < -2] = -2    
annotation_col = data.frame( CellType = factor( group_list ) )  
rownames( annotation_col ) = colnames( AssayData )  
pheatmap( fontsize = 6, 
choose_matrix, 
annotation_col = annotation_col,             show_rownames = T, 
show_colnames = F,           
annotation_legend = T, 
cluster_cols = F,             
filename = "./fig/heatmap_top100_logFC.png") }

直接贴图，后续可以分享代码讨论~

热图

火山图

气泡图

小tips:
1.setwd()很重要;
2.学会使用包很重要；
3.代码逻辑很重要;
4.了解你的数据最重要~jimmy说的，四脚赞同
从我学R的经历来看，先把基础的知识掌握好，然后再谈跑代码，这些就算放在这咱也看不懂是不。我也是按照代码一个一个查，一遍又一遍试才真正跑出来的结果，我也试过跑自己想做的方向，同样问题多多，但总归也算有结果！结果出来那一刻很开心，但学习永不止步，发现花花的公众号更简单一下，jimmy永远的神，每每看公众号文章，总能有所启发（快点学R和linux，需要伙伴一起学习讨论进步~）
由于每天都是实验室搬砖，我电脑和数据都没带回来，纯粹是靠笔记来copy和paste，后面再贴剩下的代码吧，好累，过七夕脚都要逛断了：)