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