1.安裝一些R包：
if (!require("ALL")) {BiocManager::install("ALL")}
if (!require("CLL")) {BiocManager::install("CLL")}
if (!require("pasilla")) {BiocManager::install("pasilla")}
if (!require("airway")) {BiocManager::install("airway")}
if (!require("limma")) {BiocManager::install("limma")}
if (!require("DESeq2")) {BiocManager::install("DESeq2")}
if (!require("clusterProfiler")) {BiocManager::install("clusterProfiler")}
if (!require("reshape2")) {BiocManager::install("reshape2")}
以及對(duì)這些包的一些介紹:
1、ALL： DFCI麗茲實(shí)驗(yàn)室的T細(xì)胞和B細(xì)胞急性淋巴細(xì)胞白血病數(shù)據(jù)（包括2004年4月的版本）
2.CLL: CLL軟件包包含慢性淋巴細(xì)胞性白血病（CLL）基因表達(dá)數(shù)據(jù)。CLL數(shù)據(jù)包含24個(gè)樣本，這些樣本在疾病進(jìn)展方面被分類為進(jìn)行性或穩(wěn)定。
3、pasilla：該軟件包提供了從RNA-seq數(shù)據(jù)中選擇的基因計(jì)算出的每個(gè)外顯子和每個(gè)基因的讀數(shù)計(jì)數(shù)，這些數(shù)據(jù)在Brooks AN，Yang L，Duff MO，Hansen的文章“果蠅和哺乳動(dòng)物之間的RNA調(diào)節(jié)圖的保守性”中發(fā)表。
4、airway：該程序包提供了一個(gè)RangedSummarizedExperiment對(duì)象，該對(duì)象是基因的讀取計(jì)數(shù)，用于在用地塞米松治療的四種人呼吸道平滑肌細(xì)胞系上進(jìn)行RNA-Seq實(shí)驗(yàn)。小插圖包裝中提供了有關(guān)基因模型和讀數(shù)計(jì)數(shù)程序的詳細(xì)信息。
5、limma： 數(shù)據(jù)分析，線性模型和微陣列數(shù)據(jù)的差異表達(dá)。
6、DESeq2：估計(jì)高通量測(cè)序分析中計(jì)數(shù)數(shù)據(jù)的方差均值依賴性，并基于使用負(fù)二項(xiàng)式分布的模型測(cè)試差異表達(dá)。
7、clusterProfiler：該軟件包實(shí)現(xiàn)了分析和可視化基因和基因簇的功能概況（GO和KEGG）的方法。
8、reshape2：使用兩個(gè)函數(shù)靈活地重組和聚合數(shù)據(jù):melt和'dcast'(或'acast')。
2.了解ExpressionSet對(duì)象，比如CLL包里面就有data(sCLLex) ，找到它包含的元素，提取其表達(dá)矩陣(使用exprs函數(shù))，查看其大小
1.什么是ExpressionSet對(duì)象?
在Biobase基礎(chǔ)包中，ExpressionSet是非常重要的類，因?yàn)?Bioconductor設(shè)計(jì)之初是為了對(duì)基因芯片數(shù)據(jù)進(jìn)行分 析，而ExpressionSet正是Bioconductor為基因表達(dá)數(shù)據(jù)格式所定制的標(biāo)準(zhǔn)。它是所有涉及基因表達(dá)量相關(guān)數(shù)據(jù)在 Bioconductor中進(jìn)行操作的基礎(chǔ)數(shù)據(jù)類型，比如affyPLM, affy, oligo, limma, arrayMagic等等
ExpressionSet的組成：assayData--頭文件--experimentData
suppressPackageStartupMessages(library(CLL))#加載包
data(sCLLex)
exprSet=exprs(sCLLex) ##提取表達(dá)矩陣,它列是各個(gè)樣本，行是每個(gè)探針I(yè)D
samples=sampleNames(sCLLex) ##提取名字
pdata=pData(sCLLex) ##提取臨床信息
group_list=as.character(pdata[,2]) ##強(qiáng)制轉(zhuǎn)換類型,轉(zhuǎn)換為字符型
dim(exprSet) ##查看維度.
[1] 12625 22 #行列
exprSet[1:5,1:5] #查看
###再做一個(gè)QC檢測(cè)
par(cex = 0.7) #控制圖片的文字和點(diǎn)的大小.
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)做一個(gè)箱形圖,可以查看一些特異的點(diǎn)
3.了解 str,head,help函數(shù)，作用于 第二步提取到的表達(dá)矩陣
> help('str') #str（）函數(shù)來(lái)查看這個(gè)數(shù)據(jù)框當(dāng)中有哪些變量，以及變量的類型
> help('head') #查看前面部分
> help('help') #幫助函數(shù)
str(exprSet)
head(exprSet)
4.安裝并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后 顯示的那些變量
Affymetrix人類基因組U95設(shè)置注釋數(shù)據(jù)（芯片hgu95av2），使用來(lái)自公共存儲(chǔ)庫(kù)的數(shù)據(jù)組裝而成
BiocManager::install("hgu95av2.db")
library('hgu95av2.db')
ls("package:hgu95av2.db")
str("package:hgu95av2.db")
5.理解 head(toTable(hgu95av2SYMBOL)) 的用法，找到 TP53 基因?qū)?yīng)的探針I(yè)D
ls("package:hgu95av2.db") #查看hgu95av2有什么變量
head(toTable(hgu95av2SYMBOL)) #查看hgu95av2SYMBOL的前面部分的數(shù)據(jù)
ids <- toTable(hgu95av2SYMBOL) #換為表格
TP53 <- ids[ids$symbol=='TP53',] #選擇需要的內(nèi)容
6.理解探針與基因的對(duì)應(yīng)關(guān)系，總共多少個(gè)基因，基因最多對(duì)應(yīng)多少個(gè)探針，是哪些基因，是不是因?yàn)檫@些基因很長(zhǎng)，所以在其上面設(shè)計(jì)多個(gè)探針呢？
①.探針和基因的對(duì)應(yīng)關(guān)系:
在用基因芯片中,是通過(guò)一系列已知的核酸探針,然后和帶有熒光標(biāo)記的核酸序列互補(bǔ)配對(duì),通過(guò)測(cè)定熒光強(qiáng)度最強(qiáng)的位置,獲取一組互補(bǔ)的探針序列.后期的分析,通過(guò)探針組的結(jié)果,去分析對(duì)應(yīng)的基因.一一對(duì)應(yīng).
②.總共多少個(gè)基因，基因最多對(duì)應(yīng)多少個(gè)探針，是哪些基因
attach(ids) #這樣后面訪問(wèn)的時(shí)候,就不用加$了
length(unique(symbol)) #查看長(zhǎng)度,也即是個(gè)數(shù)
[1] 8584 #一共有8584個(gè)基因symbol
table(symbol)[table(symbol)==max(table(symbol))] 查找眾數(shù)，也就是出現(xiàn)最多的元素。table可以用來(lái)對(duì)數(shù)據(jù)統(tǒng)計(jì)出現(xiàn)的次數(shù)。
7.第二步提取到的表達(dá)矩陣是12625個(gè)探針在22個(gè)樣本的表達(dá)量矩陣，找到那些不在 hgu95av2.db 包收錄的對(duì)應(yīng)著SYMBOL的探針。
noin_exprSet <- exprSet[!row.names(exprSet)%in%ids$probe_id,] #首先是提取exprSet的行名字,也就是探針,然后用%in%進(jìn)行一個(gè)邏輯判斷,是否相同.probe_id是在hgu95av2SYMBOL的行名.然后一比對(duì)
dim(noin_exprSet)
[1] 1166 22
8.過(guò)濾表達(dá)矩陣，刪除那1165個(gè)沒(méi)有對(duì)應(yīng)基因名字的探針。
exprSet <- exprSet[rownames(exprSet) %in% ids$probe_id, ] ##原理和前面一題一樣。！是取反向操作，這題不用
dim(exprSet)
9.整合表達(dá)矩陣，多個(gè)探針對(duì)應(yīng)一個(gè)基因的情況下，只保留在所有樣本里面平均表達(dá)量最大的那個(gè)探針。
提示，理解 tapply,by,aggregate,split 函數(shù) , 首先對(duì)每個(gè)基因找到最大表達(dá)量的探針。
然后根據(jù)得到探針去過(guò)濾原始表達(dá)矩陣
思路：先把是多個(gè)探針對(duì)應(yīng)一個(gè)基因的那些探針找出來(lái)——然后再比較表達(dá)量——再過(guò)濾。X
參考答案，應(yīng)該用by函數(shù)，進(jìn)行分類，把多個(gè)探針對(duì)應(yīng)一個(gè)基因的這樣的情況分類出來(lái)——再比較表達(dá)量
a <- by(exprSet,ids$symbol, ##exprSet是需要處理的數(shù)據(jù)，ids$symbol是分類的標(biāo)注。分類出每一個(gè)小的data
function(x) ##對(duì)每一個(gè)小的data進(jìn)行函數(shù)處理。x就應(yīng)該是代表著每一個(gè)分類出來(lái)的data
row.names(x)[which.max(rowMeans(x))]) ##先去平均值，再找最大的值，過(guò)濾開來(lái)。
head（a） ##看出a還是有探針和基因?qū)?yīng)的，但是我們只要探針就可以了。
INDICES
AADAC AAK1 AAMP AANAT AARS1 AASDHPPT
"36512_at" "39463_at" "38434_at" "36332_at" "36185_at" "35761_at"
b = as.character(a) ##轉(zhuǎn)換數(shù)據(jù)類型
[1] "36512_at" "39463_at" "38434_at" "36332_at" "36185_at" "35761_at" ##這樣就是有探針了，后面可以用來(lái)對(duì)exprSet進(jìn)行篩選。
exprSet <- exprSet[row.names(exprSet)%in%b,] 從前面選出來(lái)的b，對(duì)exprSet進(jìn)行篩選。
dim(exprSet)
[1] 8584 22
10.把過(guò)濾后的表達(dá)矩陣更改行名為基因的symbol，因?yàn)檫@個(gè)時(shí)候探針和基因是一對(duì)一關(guān)系了。
row.names(exprSet)=ids[rownames(exprSet)%in%ids$probe_id,2] ##邏輯判斷一下，然后更換行名。注意這個(gè)“2”。應(yīng)該是要去第二列的意思。
11、對(duì)第10步得到的表達(dá)矩陣進(jìn)行探索，先畫第一個(gè)樣本的所有基因的表達(dá)量的boxplot,hist,density ， 然后畫所有樣本的 這些圖
參考：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html
理解ggplot2的繪圖語(yǔ)法，數(shù)據(jù)和圖形元素的映射關(guān)系
思路：先畫第一個(gè)樣本的圖：看著參考，要畫的圖有：boxplot、vioplot、boxplot、histogram
density、gpairs、cluster、PCA、heatmap、DEG && volcano plot————再看看能不能做個(gè)循環(huán)，把所有的樣本都畫出來(lái)。
##箱形圖
p=ggplot(long_exprSet,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
大概是這樣的圖。x軸是22個(gè)樣本，y軸是表達(dá)量。不同顏色是每樣本人的疾病進(jìn)展情況（progres和stable）
可以得出的結(jié)論是：
①、看出每個(gè)樣本的異常值的情況。有幾個(gè)樣本的異常值特別多，比如16號(hào)樣本？？有什么用？
②、看出平均值，每個(gè)樣本的表達(dá)量的平均值的？？
###vioplot
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
什么是vioplot？
vioplot叫做小提琴圖。小提琴圖 (Violin Plot) 用于顯示數(shù)據(jù)分布及其概率密度。
得出什么結(jié)論？
①、？？？？
##直方圖
p=ggplot(long_exprSet,aes(value,fill=group))+geom_histogram()+facet_wrap(~sample,nrow = 3)##做個(gè)直方圖，其中facet_wrap是一個(gè)分面函數(shù)。
print(p)
得出的結(jié)論：
①、大部分的樣本，的表達(dá)量分布在3-4左右。
##密度圖
p=ggplot(long_exprSet,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)
可以看出，表達(dá)量value值，在大概3-4的時(shí)候，fre頻率是最大的。
gpairs？？？？看不懂怎么看這圖呀，相關(guān)性分析？
12.理解統(tǒng)計(jì)學(xué)指標(biāo)mean,median,max,min,sd,var,mad并計(jì)算出每個(gè)基因在所有樣本的這些統(tǒng)計(jì)學(xué)指標(biāo)，最后按照mad值排序，取top 50 mad值的基因，得到列表。
注意：這個(gè)題目出的并不合規(guī)，請(qǐng)仔細(xì)看。
思路：可以很簡(jiǎn)單的計(jì)算出一個(gè)基因的這些統(tǒng)計(jì)指數(shù)，如何計(jì)算多個(gè)基因的呢？？
——所以需要用到apply函數(shù)才可以
mean <- apply(exprSet, 1,mean)
mean <- as.data.frame(mean)
median <- apply(exprSet, 1,median)
median <- as.data.frame(median)
max <- apply(exprSet, 1,max)
max <- as.data.frame(max)
min <- as.data.frame(apply(exprSet, 1,min))
var <- as.data.frame(apply(exprSet, 1,var))
mad <- as.data.frame(apply(exprSet, 1,mad))
merge(mean,median,max,min,sd,var,mad,by=rownames)？？想合并起來(lái)，看起來(lái)會(huì)比較好看，不過(guò)好像不行，因?yàn)閐im（）看出來(lái)，只有一列。還是算了，不合并了。
g_mad <- tail(sort(apply(exprSet,1,mad)),50)
names(g_mad)
13.根據(jù)第12步驟得到top 50 mad值的基因列表來(lái)取表達(dá)矩陣的子集，并且熱圖可視化子表達(dá)矩陣。試試看其它5種熱圖的包的不同效果。
思路：先根據(jù)g_mad來(lái)在exprest里面取出部分?jǐn)?shù)據(jù)——然后再對(duì)這份數(shù)據(jù)進(jìn)行做個(gè)熱圖。
TOP50 <- exprSet[g_mad,]
library(pheatmap)
pheatmap(TOP50)
14.取不同統(tǒng)計(jì)學(xué)指標(biāo)mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包來(lái)看他們之間的overlap情況。
UpSetR包：對(duì)集合，進(jìn)行可視化的一個(gè)包。
思路：怎么去做一個(gè)，符合UpSetR包的表格。
g_all <- unique(c(names(g_mean),names(g_median),names(g_max),names(g_min),
names(g_sd),names(g_var),names(g_mad) )) ###行名
dat=data.frame(g_all=g_all,
g_mean=ifelse(g_all %in% names(g_mean) ,1,0), ##ifelse的一個(gè)循環(huán)。
g_median=ifelse(g_all %in% names(g_median) ,1,0),
g_max=ifelse(g_all %in% names(g_max) ,1,0),
g_min=ifelse(g_all %in% names(g_min) ,1,0),
g_sd=ifelse(g_all %in% names(g_sd) ,1,0),
g_var=ifelse(g_all %in% names(g_var) ,1,0),
g_mad=ifelse(g_all %in% names(g_mad) ,1,0)
)
upset(dat,nsets = 7)
點(diǎn)和點(diǎn)相連的就是有交集。交集的個(gè)數(shù)就往上看柱形圖。
15.在第二步的基礎(chǔ)上面提取CLL包里面的data(sCLLex) 數(shù)據(jù)對(duì)象的樣本的表型數(shù)據(jù)。
pdata=pData(sCLLex)
dim(pdata)
Disease <- pdata$Disease
print(Disease)
16.對(duì)所有樣本的表達(dá)矩陣進(jìn)行聚類并且繪圖，然后添加樣本的臨床表型數(shù)據(jù)信息(更改樣本名)
思路：如何聚類？——繪圖？——更改樣本名。
怎么看聚類圖。從大的往里面切，看看自己分幾類
那分類的依據(jù)是什么？不知道用的是什么算法，但分類的依據(jù)應(yīng)該是，基因表達(dá)量。
colnames(exprSet)=paste(group_list,1:22,sep='') ##把CLL11.CEL——換為progres.這樣的病例信息
exprSet <- t(exprSet) ##需要吧行和列轉(zhuǎn)換一下。因?yàn)榻酉聛?lái)要多個(gè)dist，求距離的，是對(duì)dist的是計(jì)算每一行的之間的距離，然后是以每一列的作為變量作為標(biāo)注
hc <- hclust(dist(exprSet)) ##聚類
plot(as.dendrogram(hc), horiz = TRUE) ##horiz = TRUE是吧豎向的換為橫向的。
怎么看這個(gè)圖：從左往右看，看看那些樣本是一大類的。——這樣又能得出什么有用的信息呢？
17.對(duì)所有樣本的表達(dá)矩陣進(jìn)行PCA分析并且繪圖，同樣要添加表型信息。
思路：什么是PCA分析。ggfortify這個(gè)包怎么用。
PCA分析：主成分分析，數(shù)據(jù)的降維？看疾病組和對(duì)照組能否分開
install.packages("ggfortify")
library(ggfortify)
exprSet <- exprs(sCLLex)
df <- as.data.frame(t(exprSet))
df$group <- group_list
autoplot(prcomp(df[,1:(ncol(df)-1)]), data=df, colour = 'group') ？？？
第一成分和第二成分加起來(lái)才30%左右，證明這個(gè)是不太好的？
另外一個(gè)就是progres組和stable組是沒(méi)有明顯差異的，因?yàn)樗麄冇薪患?，距離比較近。是否有差異還得看p值
18.根據(jù)表達(dá)矩陣及樣本分組信息進(jìn)行批量T檢驗(yàn)，得到檢驗(yàn)結(jié)果表格
如何進(jìn)行T檢驗(yàn)？
還是不懂？？雙樣本的檢驗(yàn)？
19.使用limma包對(duì)表達(dá)矩陣及樣本分組信息進(jìn)行差異分析，得到差異分析表格，重點(diǎn)看logFC和P值，畫個(gè)火山圖(就是logFC和-log10(P值)的散點(diǎn)圖)。
？？如何用這個(gè)包？
對(duì)哪些數(shù)據(jù)進(jìn)行分析？
suppressMessages(library(limma)) #使用limma包進(jìn)行差異基因分析時(shí)，做最多的是兩分類的，例如control組和disease組，但也會(huì)碰到按照序列進(jìn)行的分組。這時(shí)，如果逐一使用兩兩比較求差異基因則略顯復(fù)雜。其實(shí)開發(fā)limma包的大神們已經(jīng)替我們考慮到
design <- model.matrix(~0+factor(group_list)) ##構(gòu)建分組矩陣
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design) ##構(gòu)建比較矩陣
contrast.matrix
##這個(gè)矩陣聲明，我們要把progres.組跟stable進(jìn)行差異分析比較
##step1
fit <- lmFit(exprSet,design)
fit
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##這一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 去掉NA值
head(nrDEG )
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
## volcano plot
DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
'\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)
g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab("log2 fold change") + ylab("-log10 p-value") +
ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c('blue','black','red')) ## corresponding to the levels(res$change)
print(g)
