續(xù)昨

par()函數(shù)

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.
很強(qiáng)大參數(shù)較多：
cex用于表示對默認(rèn)的繪圖文本和符號放大多少倍，需要注意一些繪圖函數(shù)如plot.default等也有一個相同名字的參數(shù)，但是此時表示在函數(shù)par()的參數(shù)cex的基礎(chǔ)上再放大多少倍，此外還有函數(shù)points等接受一個數(shù)值向量為參數(shù)
cex.axis：表示在當(dāng)前的cex設(shè)定情況下，對坐標(biāo)軸刻度值字體的放大倍數(shù)
cex.lab：表示在當(dāng)前的cex設(shè)定情況下，對坐標(biāo)軸名稱字體的放大倍數(shù)
cex.main：表示在當(dāng)前的cex設(shè)定情況下，對主標(biāo)題字體的放大倍數(shù)
cex.sub：表示在當(dāng)前的cex設(shè)定情況下，對子標(biāo)題字體的放大倍數(shù)

str()

structure：緊湊的顯示對象內(nèi)部結(jié)構(gòu)，即對象內(nèi)有什么

ls+(packages.name)

類似linux用法

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

完整的代碼流程，與上一節(jié)有部分重復(fù)

rm(list = ls()) ###清除所有的環(huán)境變量
options(stringsAsFactors = F)
suppressPackageStartupMessages(library(CLL))
# 這個數(shù)據(jù)集中有22個樣本，12625個基因；使用exprs可以查看這個數(shù)據(jù)集的表達(dá)水平；
#得到表達(dá)矩陣
data(sCLLex)
class(sCLLex)
dim(sCLLex)
colnames(sCLLex)
df = as.data.frame(sCLLex)
exprSet=exprs(sCLLex)  #exprs函數(shù)提取出來表達(dá)矩陣，賦給data_expression

samples=sampleNames(sCLLex) #查看樣本編號
pdata=pData(sCLLex)# 查看分組信息
varMetadata(sCLLex) #查看所有表型變量
data_phenotype=pData(sCLLex)#提取表型信息
featureNames(sCLLex)[1:10] #查看基因芯片編碼
library(dplyr)
featureNames(sCLLex) %>% unique() %>% length() # 查看是否有重復(fù)的編碼
data_expression <- exprs(sCLLex)
group_list=as.character(pdata$Disease)
table(group_list) #統(tǒng)計(jì)表型信息
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
}

###表達(dá)矩陣的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)
  
}
##### 繪制芯片數(shù)據(jù)的質(zhì)量值(類似上文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("表達(dá)值") + guides(fill=FALSE)
p
str(exprSet)

###制作差異比較矩陣，比較矩陣之間用“-”連接
if (T) {
  contrast.matrix<-makeContrasts(paste0(unique(group_list),
                                        collapse = "-"),levels = design)
}

#lmFit()：線性擬合模型構(gòu)建【需要兩個東西：exprSet和design】
#得到的結(jié)果再和contrast一起導(dǎo)入contrasts.fit()函數(shù)
fit <- lmFit(exprSet,design)
#contrast.fit需要fit及分組矩陣，分組矩陣由makeContrasts()得到
fit2 <- contrasts.fit(fit, contrast.matrix) 
#eBayes()：利用上一步contrasts.fit()的結(jié)果
fit2 <- eBayes(fit2)
#利用上一步eBayes()的結(jié)果，并導(dǎo)出差異分析結(jié)果
tempOutput = topTable(fit2, coef=1, n=Inf)
#去掉na
nrDEG = na.omit(tempOutput) 
head(nrDEG)
write.csv(nrDEG,"limma_notrend.results.csv",quote = F)

題外話 智通寺對聯(lián) 身前有余忘縮手 眼前無路想回頭 果然每逢長假總是機(jī)緣巧合會再看紅樓夢（在想要不要專寫一個紅樓夢系列專題額---從服飾鑒賞、飲食起居、話語機(jī)鋒到詩詞折子戲等等，總是犯懶，可能也就是想想罷了）

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