scRNA---Day4

续昨

par()函数

Graphical Parameters：par can be used to set or query graphical parameters. Parameters can be set by specifying them as arguments to par in tag = value form, or by passing them as a list of tagged values.
很强大参数较多：
cex用于表示对默认的绘图文本和符号放大多少倍，需要注意一些绘图函数如plot.default等也有一个相同名字的参数，但是此时表示在函数par()的参数cex的基础上再放大多少倍，此外还有函数points等接受一个数值向量为参数
cex.axis：表示在当前的cex设定情况下，对坐标轴刻度值字体的放大倍数
cex.lab：表示在当前的cex设定情况下，对坐标轴名称字体的放大倍数
cex.main：表示在当前的cex设定情况下，对主标题字体的放大倍数
cex.sub：表示在当前的cex设定情况下，对子标题字体的放大倍数

str()

structure：紧凑的显示对象内部结构，即对象内有什么

ls+(packages.name)

类似linux用法

library("hgu95av2.db")
ls("package:hgu95av2.db")

完整的代码流程，与上一节有部分重复

rm(list = ls()) ###清除所有的环境变量
options(stringsAsFactors = F)
suppressPackageStartupMessages(library(CLL))
# 这个数据集中有22个样本，12625个基因；使用exprs可以查看这个数据集的表达水平；
#得到表达矩阵
data(sCLLex)
class(sCLLex)
dim(sCLLex)
colnames(sCLLex)
df = as.data.frame(sCLLex)
exprSet=exprs(sCLLex)  #exprs函数提取出来表达矩阵，赋给data_expression

samples=sampleNames(sCLLex) #查看样本编号
pdata=pData(sCLLex)# 查看分组信息
varMetadata(sCLLex) #查看所有表型变量
data_phenotype=pData(sCLLex)#提取表型信息
featureNames(sCLLex)[1:10] #查看基因芯片编码
library(dplyr)
featureNames(sCLLex) %>% unique() %>% length() # 查看是否有重复的编码
data_expression <- exprs(sCLLex)
group_list=as.character(pdata$Disease)
table(group_list) #统计表型信息
group_list=as.character(pdata[,2])

###得到分组矩阵
if(T) {suppressMessages(library(limma))
  design <- model.matrix(~0+factor(group_list))
  colnames(design)=levels(factor(group_list))
  rownames(design)=colnames(exprSet)
  design
}

###表达矩阵的QC检测
if (T) {par(cex = 0.7)
  n.sample=ncol(exprSet)
  if(n.sample>40) par(cex = 0.5)
  cols <- rainbow(n.sample*1.2)
  boxplot(exprSet, col = cols,main="expression value",las=2)
  
}
##### 绘制芯片数据的质量值(类似上文QC检测)
library(ggplot2)
library(reshape2)
library(gpairs)
library(corrplot)
y <- melt(as.data.frame(data_expression))
p <- ggplot(data=y,aes(x=variable,y=value))
p <- p + geom_boxplot(aes(fill=variable))
p <- p + theme(axis.text.x = element_text(angle=45,hjust=1,vjust=1))
p <- p + xlab("分组信息") + ylab("表达值") + guides(fill=FALSE)
p
str(exprSet)

###制作差异比较矩阵，比较矩阵之间用“-”连接
if (T) {
  contrast.matrix<-makeContrasts(paste0(unique(group_list),
                                        collapse = "-"),levels = design)
}

#lmFit()：线性拟合模型构建【需要两个东西：exprSet和design】
#得到的结果再和contrast一起导入contrasts.fit()函数
fit <- lmFit(exprSet,design)
#contrast.fit需要fit及分组矩阵，分组矩阵由makeContrasts()得到
fit2 <- contrasts.fit(fit, contrast.matrix) 
#eBayes()：利用上一步contrasts.fit()的结果
fit2 <- eBayes(fit2)
#利用上一步eBayes()的结果，并导出差异分析结果
tempOutput = topTable(fit2, coef=1, n=Inf)
#去掉na
nrDEG = na.omit(tempOutput) 
head(nrDEG)
write.csv(nrDEG,"limma_notrend.results.csv",quote = F)

题外话智通寺对联身前有余忘缩手眼前无路想回头果然每逢长假总是机缘巧合会再看红楼梦（在想要不要专写一个红楼梦系列专题额---从服饰鉴赏、饮食起居、话语机锋到诗词折子戏等等，总是犯懒，可能也就是想想罢了）

getwd()
D:/something/scRNA_smart_seq2-master/limma