生信人的20個(gè)R語言習(xí)題

1 安裝R包

數(shù)據(jù)包： ALL, CLL, pasilla, airway 
軟件包：limma，DESeq2，clusterProfiler 
工具包：reshape2
繪圖包：ggplot2

options()$repos 
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install(c("ALL","CLL","pasilla", "airway") ,ask = F,update = F)
BiocManager::install(c("limma","DESeq2", "clusterProfiler") ,ask = F,update = F)
install.packages("reshape2") #工具包
install.packages("ggplot2") #繪圖包

2 了解ExpressionSet對(duì)象

CLL包里面就有data(sCLLex) 计寇，找到它包含的元素元莫，提取其表達(dá)矩陣(使用exprs函數(shù))踱蠢，查看其大小

1. 參考：http://www.bio-info-trainee.com/bioconductor_China/software/limma.html
2. 參考：https://github.com/bioconductor-china/basic/blob/master/ExpressionSet.md

suppressPackageStartupMessages(library(CLL)) #隱藏具體加載信息
data(sCLLex)
exprSet=exprs(sCLLex)   ##sCLLex是依賴于CLL這個(gè)package的一個(gè)對(duì)象
samples=sampleNames(sCLLex)
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
dim(exprSet)

3 了解 str,head,help函數(shù)，作用于第二步提取到的表達(dá)矩陣

str(exprSet)
head(exprSet) 
help("exprSet")

4 安裝并了解hgu95av2.db包

安裝并了解 hgu95av2.db 包,看看 ls("package:hgu95av2.db") 后顯示的那些變量

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

image

5.理解 head(toTable(hgu95av2SYMBOL)) 的用法

理解 head(toTable(hgu95av2SYMBOL)) 的用法,找到 TP53 基因?qū)?yīng)的探針I(yè)D

head(toTable(hgu95av2SYMBOL)) 
   probe_id  symbol
1   1000_at   MAPK3
2   1001_at    TIE1
3   1002_f_at CYP2C19
4   1003_s_at   CXCR5
5   1004_at   CXCR5
6   1005_at   DUSP1
ids=toTable(hgu95av2SYMBOL)
#在ids中搜索TP53
> ids[which(ids[,2]=="TP53"),]#記住
      probe_id symbol
966    1939_at   TP53
997  1974_s_at   TP53
1420  31618_at   TP53
> ids[grep("^TP53$",ids$symbol),]#記住
      probe_id symbol
966    1939_at   TP53
997  1974_s_at   TP53
1420  31618_at   TP53

基因TP53對(duì)應(yīng)的探針

6.理解探針與基因的對(duì)應(yīng)關(guān)系

理解探針與基因的對(duì)應(yīng)關(guān)系播聪，總共多少個(gè)基因朽基，基因最多對(duì)應(yīng)多少個(gè)探針，是哪些基因离陶，是不是因?yàn)檫@些基因很長，所以在其上面設(shè)計(jì)多個(gè)探針呢招刨？

summary(hgu95av2SYMBOL)
SYMBOL map for chip hgu95av2 (object of class "ProbeAnnDbBimap")
|
| Lkeyname: probe_id (Ltablename: probes)
|    Lkeys: "1000_at", "1001_at", ... (total=12625/mapped=11460)
|
| Rkeyname: symbol (Rtablename: gene_info)
|    Rkeys: "A1BG", "A2M", ... (total=61050/mapped=8585)
|
| direction: L --> R
#左邊: 探針總數(shù)是12625霎俩，能匹配上的是11460個(gè)；
#右邊: 基因名的總數(shù)是30071沉眶，實(shí)際上只有8585種
#或者
mapped_probes <- mappedkeys(hgu95av2SYMBOL) #手動(dòng)過濾打却，然后進(jìn)行統(tǒng)計(jì)
count.mappedkeys(hgu95av2SYMBOL) #統(tǒng)計(jì)：與summary結(jié)果相符
[1] 11460
#有1165個(gè)探針沒有匹配上基因。（12625-11460=1165）

ids <- toTable(hgu95av2SYMBOL)
1.dim(ids)
[1] 11460     2
2.length(unique(ids$symbol))
[1] 8585
3.tail(sort(table(ids$symbol)))
YME1L1  GAPDH INPP4A    MYB PTGER3  STAT1 
     7      8      8      8      8      8 
4.table(sort(table(ids$symbol)))#sort(table(ids$symbol))統(tǒng)計(jì)的是每個(gè)基因?qū)?yīng)的探針數(shù)目谎倔。
   1    2    3    4    5    6    7    8 
6555 1428  451  102   22   16    6    5 
##一個(gè)探針對(duì)應(yīng)唯一一個(gè)基因的為6555個(gè)探針柳击，也就是6555個(gè)基因是對(duì)應(yīng)唯一一個(gè)探針的。還有一個(gè)基因?qū)?yīng)多個(gè)探針的片习，比如最后一列捌肴，就是一個(gè)基因?qū)?yīng)8個(gè)探針，這樣的基因共有5個(gè)藕咏。
#所以其實(shí)6555+1428+451+102+22+16+6+5=8585状知，1428、451等是基因數(shù)目孽查，只不過對(duì)應(yīng)的探針數(shù)目不同而已饥悴，但也是unique的，因此同樣應(yīng)該是算在unique(ids$symbol)里面的盲再。

7.找到那些不在 hgu95av2.db 包收錄的對(duì)應(yīng)著SYMBOL的探針西设。

第二步提取到的表達(dá)矩陣是12625個(gè)探針在22個(gè)樣本的表達(dá)量矩陣，找到那些不在 hgu95av2.db 包收錄的對(duì)應(yīng)著SYMBOL的探針答朋。

從上一步了解到有1165個(gè)探針是不能匹配上基因名字的贷揽，現(xiàn)在就是要找出這些探針。

#第二步獲得的exprSet這個(gè)表達(dá)矩陣
dim(exprSet)
[1] 12625    22
#以下為e1和e2在 hgu95av2.db 包收錄的對(duì)應(yīng)著SYMBOL的探針绿映。
e1<-exprSet[rownames(exprSet) %in% ids$probe_id,]
dim(e1)
[1] 11460    22
e2<-exprSet[match(rownames(exprSet),ids$probe_id, nomatch = 0),]
dim(e2)
[1] 11460    22
#但要找1165個(gè)不能匹配上基因名字的探針
e3<-exprSet[!rownames(exprSet) %in% ids$probe_id,]
dim(e3)
[1] 1165   22
chpro_id<-as.data.frame(row.names(e3))
> dim(chpro_id)
[1] 1165    1

不在hgu95av2.db包收錄的對(duì)應(yīng)著SYMBOL的探針

8.過濾表達(dá)矩陣，刪除那1165個(gè)沒有對(duì)應(yīng)基因名字的探針。

#思路是不刪除叉弦，直接取有基因名字的探針為一個(gè)表達(dá)矩陣丐一。同第7題中的e1和e2。
e1<-exprSet[rownames(exprSet) %in% ids$probe_id,]
e2<-exprSet[match(rownames(exprSet),ids$probe_id, nomatch = 0),]

9.整合表達(dá)矩陣

整合表達(dá)矩陣淹冰，多個(gè)探針對(duì)應(yīng)一個(gè)基因的情況下库车，只保留在所有樣本里面平均表達(dá)量最大的那個(gè)探針。
提示樱拴，理解 tapply,by,aggregate,split 函數(shù) , 首先對(duì)每個(gè)基因找到最大表達(dá)量的探針柠衍，然后根據(jù)得到探針去過濾原始表達(dá)矩陣。

#新建dat
dat=exprSet[rownames(exprSet) %in% ids$probe_id,]
ids$median=apply(dat,1,median)#對(duì)dat每一行求median值
ids=ids[order(ids$symbol,ids$median,decreasing = T),]#order是按從小到達(dá)排序晶乔，加了decreasing = T是按從大到小排序珍坊。ids$symbol是將symbol基因名按照字母順序從最后即"z"排序的，然后每個(gè)z開頭的再按ids$median排序
ids=ids[!duplicated(ids$symbol),]#將ids$symbol去重后作為新的行正罢，ids取這些行
dat=dat[ids$probe_id,]
dim(dat)
[1] 8585   22 #此時(shí)這個(gè)8585就是前面第6題中unique(ids$symbol)的基因的數(shù)目

10.把過濾后的表達(dá)矩陣更改行名為基因的symbol阵漏，因?yàn)檫@個(gè)時(shí)候探針和基因是一對(duì)一關(guān)系了。

rownames(dat)<-ids$symbol

dat

11.對(duì)第10步得到的表達(dá)矩陣進(jìn)行探索

對(duì)第10步得到的表達(dá)矩陣進(jìn)行探索翻具，先畫第一個(gè)樣本的所有基因的表達(dá)量的boxplot,hist,density履怯，然后畫所有樣本的這些圖。
參考：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html 理解ggplot2的繪圖語法裆泳，數(shù)據(jù)和圖形元素的映射關(guān)系叹洲。

##需要加載以下包
library(CLL)
library(ggplot2)
library(reshape2)
library(gpairs)
library(corrplot)
##把我們的寬矩陣用reshape2包變成長矩陣,用melt這個(gè)函數(shù)
library(reshape2)
exprSet<-dat #后面的代碼是exprSet，因此又轉(zhuǎn)換成exprSet
exprSet_L=melt(exprSet)#后面的矩陣都為exprSet_L
colnames(exprSet_L)=c('probe','sample','value')
exprSet_L$group=rep(group_list,each=nrow(dat))
head(exprSet_L)
 probe    sample    value    group
1   ZZZ3 CLL11.CEL 6.645791 progres.
2  ZZEF1 CLL11.CEL 5.289264 progres.
3    ZYX CLL11.CEL 3.949769 progres.
4  ZWINT CLL11.CEL 4.316881 progres.
5   ZW10 CLL11.CEL 4.382004 progres.
6 ZSWIM8 CLL11.CEL 4.091876 progres.

p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print(p)
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print(p)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print(p)

image

image

image

image

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值的基因，得到列表帜篇。

e_mean = tail(sort(apply(exprSet,1,mean)),30)
e_median = tail(sort(apply(exprSet,1,median)), 30)
e_max <- tail(sort(apply(exprSet,1,max)),30)
e_min <- tail(sort(apply(exprSet,1,min)),30)
e_sd <- tail(sort(apply(exprSet,1,sd)),30)
e_var <- tail(sort(apply(exprSet,1,var)),30)
e_mad <- tail(sort(apply(exprSet,1,mad)),30) #絕對(duì)中位差來估計(jì)方差,先計(jì)算出數(shù)據(jù)與它們的中位數(shù)之間的偏差糙捺，然后這些偏差的絕對(duì)值的中位數(shù)就是mad

13.根據(jù)第12步驟得到top 50 mad值的基因列表來取表達(dá)矩陣的子集，并且熱圖可視化子表達(dá)矩陣笙隙。試試看其它5種熱圖的包的不同效果洪灯。

dat_top50<-dat[tail(sort(apply(exprSet,1,mad)),50),] #圖左
choose_gene=names(sort(apply(exprSet, 1, mad),decreasing = T)[1:50])
choose_matrix=exprSet[choose_gene,]
choose_matrix=scale(choose_matrix)
heatmap(choose_matrix)

library(gplots) #圖右上
heatmap.2(choose_matrix)

library(pheatmap) #圖右下
choose_gene=names(tail(sort(apply(exprSet,1,mad)),50))
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)

image

14.取不同統(tǒng)計(jì)學(xué)指標(biāo)mean,median,max,mean,sd,var,mad的各top50基因列表，使用UpSetR包來看他們之間的overlap情況竟痰。

library(UpSetR)
g_all <- unique(c(names(e_mean),names(e_median),names(e_max),names(e_min),
                  names(e_sd),names(e_var),names(e_mad) ))
dat=data.frame(e_all=g_all,
               e_mean=ifelse(g_all %in%  names(e_mean) ,1,0),
               e_median=ifelse(g_all %in%  names(e_median) ,1,0),
               e_max=ifelse(g_all %in%  names(e_max) ,1,0),
               e_min=ifelse(g_all %in%  names(e_min) ,1,0),
               e_sd=ifelse(g_all %in%  names(e_sd) ,1,0),
               e_var=ifelse(g_all %in%  names(e_var) ,1,0),
               e_mad=ifelse(g_all %in%  names(e_mad) ,1,0)
)
upset(dat,nsets = 7)

overlap

15.在第二步的基礎(chǔ)上面提取CLL包里面的data(sCLLex) 數(shù)據(jù)對(duì)象的樣本的表型數(shù)據(jù)签钩。

pd = pData(sCLLex)
group_list = as.character(pd[,2])
table(group_list)
group_list
progres.   stable 
      14        8 

16.對(duì)所有樣本的表達(dá)矩陣進(jìn)行聚類并且繪圖，然后添加樣本的臨床表型數(shù)據(jù)信息(更改樣本名)

out.dist=dist(t(exprSet),method='euclidean') #用dist這個(gè)函數(shù)計(jì)算樣本間的距離坏快，這里的t(exprSet)是因?yàn)闃颖疽谛星﹂荩壳暗膃xprSet是樣本在列，基因名字在行莽鸿，所以t(exprSet)將行和列換位置
out.hclust=hclust(out.dist,method='complete')
plot(out.hclust)

'complete'

out.dist=dist(t(exprSet),method='euclidean')
out.hclust=hclust(out.dist,method='ward.D2')
plot(out.hclust)

'ward.D2'

可以在hclust()里面調(diào)整參數(shù)method昧旨，選擇不同的方法：

hclust(E.dist, method="single") # 最近法
hclust(E.dist, method="complete") # 最遠(yuǎn)法
hclust(E.dist, method="average") # 平均法
hclust(E.dist, method="centroid") # 中心法
hclust(E.dist, method="ward.D2") # 華德法

17.對(duì)所有樣本的表達(dá)矩陣進(jìn)行PCA分析并且繪圖拾给，同樣要添加表型信息。

pc <- prcomp(t(exprSet),scale=TRUE) #exprSet本是行為基因名兔沃，列為樣本名蒋得，但t(exprSet)后就應(yīng)該反過來了
pcx=data.frame(pc$x)
pcr=cbind(samples=rownames(pcx),group_list, pcx) 
p=ggplot(pcr, aes(PC1, PC2))+geom_point(size=5, aes(color=group_list)) +#能最大反映樣本差異性的兩個(gè)成分（PC1、PC2）
  geom_text(aes(label=samples),hjust=-0.1, vjust=-0.3)#label=samples可以加上樣本名稱
print(p)

PCA

心得：由于前面有一步將exprSet給轉(zhuǎn)換成行為樣本名乒疏，列為基因名了额衙，所以這時(shí)再t(exprSet)就是錯(cuò)誤的了，用的原始的exprSet都是行為基因名怕吴，列為樣本窍侧。在這個(gè)地方搗鼓了好久。

18.根據(jù)表達(dá)矩陣及樣本分組信息進(jìn)行批量T檢驗(yàn)转绷，得到檢驗(yàn)結(jié)果表格

exprSet=dat #從以上步驟得到的是 先是dat獲得的行為基因名伟件，列為正常順序的樣本的這樣一個(gè)矩陣，給exprSet后暇咆，它倆現(xiàn)在一樣了
group_list=as.factor(group_list)
group1 = which(group_list == levels(group_list)[1])#見下圖中progress.中所示
group2 = which(group_list == levels(group_list)[2])#見下圖中stable中所示
dat1 = dat[, group1]#取出分組信息為progress.的表達(dá)矩陣
dat2 = dat[, group2]##取出分組信息為stable的表達(dá)矩陣
dat = cbind(dat1, dat2)#cbind：列數(shù)形同時(shí)锋爪，橫向追加。rbind：行數(shù)相同時(shí)爸业，縱向追加
pvals = apply(exprSet, 1, function(x){ #后面導(dǎo)致19題出現(xiàn)的問題是這個(gè)地方的exprSet在之前被我給換成rbind的dat了其骄，也就是順序是14個(gè)progress.（11、13扯旷、14拯爽、15、16钧忽、19...）毯炮，然乎再全是stable。而此時(shí)分組信息group_list就是按照最初group_list中的12耸黑、13桃煎、14、15大刊、16为迈、16...。由于矩陣和分組信息不是匹配的缺菌，所以不該該exprSet為cbind后的dat葫辐。那么問題是，既然用exprSet伴郁，那么為什么要有dat = cbind(dat1, dat2)這步呢耿战？這時(shí)目前唯一不理解的啦，但是如果不用這個(gè)exprSet（行為基因名焊傅，列為正常順序的樣本）剂陡，后面第20題就會(huì)做出錯(cuò)誤的圖狈涮。
  t.test(as.numeric(x)~group_list)$p.value
})
p.adj = p.adjust(pvals, method = "BH")
avg_1 = rowMeans(dat1)
avg_2 = rowMeans(dat2)
log2FC = avg_2-avg_1
DEG_t.test = cbind(avg_1, avg_2, log2FC, pvals, p.adj)
DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),]
DEG_t.test=as.data.frame(DEG_t.test)
# 查看t檢驗(yàn)結(jié)果表格，包含log2FC鸭栖、pvals和p.adj等薯嗤，通常認(rèn)為t<0.05即有統(tǒng)計(jì)學(xué)意義
head(DEG_t.test)

• 結(jié)果為

DEG_t.test

找到分組信息為progres.的下角標(biāo)（即在dat矩陣中所處于的位置（哪一列））

progres.

找到分組信息為stable的下角標(biāo)（即在dat矩陣中所處于的位置（哪一列））

stable

• 簡單羅列一下which 和 %in% 的用法，如下舉例

a=c(1,3,4,5,3,2,5,6,3,2,5,6,7,5,8)  
which.max(a)  #which列出符合函數(shù)結(jié)果的下角標(biāo)纤泵，即列出默認(rèn)為TRUE的下角標(biāo)
[1] 15
10:1 %in% a
[1] FALSE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE
b <- which(10:1 %in% a)  
b
[1]  3  4  5  6  7  8  9 10

• cbind函數(shù)用法，如下舉例

> a <- matrix(1:6, 2, 3)
> print(a)
     [,1] [,2] [,3] [,4]
[1,]    1    4    7   10
[2,]    2    5    8   11
[3,]    3    6    9   12
> 
> b <- matrix(-1:-12, 3, 4)
> print(b)
     [,1] [,2] [,3] [,4]
[1,]   -1   -4   -7  -10
[2,]   -2   -5   -8  -11
[3,]   -3   -6   -9  -12
> 
> x=cbind(a,b)
> print(x)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    1    4    7   10   -1   -4   -7  -10
[2,]    2    5    8   11   -2   -5   -8  -11
[3,]    3    6    9   12   -3   -6   -9  -12

19.使用limma包對(duì)表達(dá)矩陣及樣本分組信息進(jìn)行差異分析镜粤，得到差異分析表格捏题，重點(diǎn)看logFC和P值，畫個(gè)火山圖(就是logFC和-log10(P值)的散點(diǎn)圖肉渴。

下面代碼塊中的dat是18題得到的dat = cbind(dat1, dat2)的dat公荧，把它賦值給exprSet。這就是后面為什么出現(xiàn)的火山圖的差異表達(dá)基因和其他人不同的原因了

#exprSet=dat #這一步也不該加同规，因?yàn)榇藭r(shí)的dat是18題中cbind后的dat循狰，而我之前加了，加了后矩陣和分組信息又是不對(duì)應(yīng)的了券勺，見下圖中錯(cuò)誤的
exprSet[1:4,1:4]
# DEG by limma 
suppressMessages(library(limma)) 
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design #注意此步的design圖片與下張design圖片是不同的绪钥，下面的exprSet才是正確的
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix<-makeContrasts("progres.-stable",levels = design)#關(guān)于makeContrasts下面有個(gè)示例
contrast.matrix 
##這個(gè)矩陣聲明，我們要把progres.組跟stable進(jìn)行差異分析比較
##step1
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) 
fit2 <- eBayes(fit2)  
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)

#下面這些代碼是調(diào)整好的关炼，基本可以不用替換的呢
DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs(logFC)) + 2*sd(abs(logFC)) ) 
DEG$result = 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), #round保留小數(shù)位數(shù)
                    '\nThe number of up gene is ',nrow(DEG[DEG$result =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$result =='DOWN',])
)
ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=result)) +
  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'))

行名是按照順序的

錯(cuò)誤的design程腹，因?yàn)榉纸M信息的1、0都是一直沒變動(dòng)的group_list

錯(cuò)誤的nrDEG,因?yàn)榇藭r(shí)的用的是cbind后的dat

錯(cuò)誤火山圖

• 插播插曲

接下來下面3張圖片的代碼與上面的代碼相同儒拂，只不過第一步中exprSet=dat的dat是未經(jīng)過第18題中的dat = cbind(dat1, dat2)的dat寸潦，而是第10題中的dat，就是說社痛，用的exprSet就是原始從第2題中提取出的表達(dá)矩陣见转，經(jīng)過同ids對(duì)應(yīng)好以后的exprSet。因此樣本信息也是按照沒有修改過的順序11蒜哀、12斩箫、13、14凡怎、15來的校焦，同時(shí)design信息即樣本和分組信息的對(duì)應(yīng)是正確的，因?yàn)檫@個(gè)里面的exprSet是最初的统倒，其實(shí)group_list一直都沒有變過寨典，變的是矩陣，即上面的dat是dat1和dat2組合在一起的房匆，它的樣本順序是改變了的耸成，是dat1全部是progres.,而dat2全部是stable报亩。但是我們并沒有操作group_list,group_list是最開始第2題中從pData中提取出來的。所以在上面圖中design的樣本和分組信息是錯(cuò)誤對(duì)應(yīng)的井氢。因此上面最后畫出的火山圖也是錯(cuò)誤的弦追。總結(jié)來說花竞，這里面的分組信息和矩陣exprSet是獨(dú)立來的劲件，可以class(exprSet)可以得到matrix，而不是data frame约急，所以在把dat分為progres.和stable后在cbind這兩組后零远，并不能讓group_list跟著一起分組。

#下面生成圖的代碼同上面的代碼一樣厌蔽，只不過第一步中exprSet=dat的dat是未經(jīng)過第18題中的dat = cbind(dat1, dat2)的dat牵辣，而是第10題中的dat，即行為基因名字   奴饮，列為同exprSet中一樣的順序纬向。

exprSet

正確的design

正確nrDEG

正確火山圖

• 代碼中設(shè)計(jì)的函數(shù)的一些簡單示例，如下
• makeContrasts函數(shù)用法戴卜，如下舉例

x <- c("B-A","C-B","C-A")
makeContrasts(contrasts=x,levels=c("A","B","C"))
  Contrasts
Levels B - A C - B C - A
     A    -1     0    -1
     B     1    -1     0
     C     0     1     1

#相當(dāng)于下面這個(gè)逾条，輸出結(jié)果是一樣的
makeContrasts(B-A,C-B,C-A,levels=c("A","B","C"))

   

• unique、duplicated函數(shù)用法投剥，如下舉例

#unique主要是返回一個(gè)把重復(fù)元素或行給刪除的向量膳帕、數(shù)據(jù)框或數(shù)組
> rt
    年 月 公司名 利率
1 2000  1      A    a
2 2000  1      A    a
3 2001  2      A    b
4 2001  3      A    c
5 2000  1      B    d
6 2000  2      B    e
7 2000  2      B    e
> unique(rt)
    年 月 公司名 利率
1 2000  1      A    a
3 2001  2      A    b
4 2001  3      A    c
5 2000  1      B    d
6 2000  2      B    e
> unique(rt,fromLast=TRUE)
    年 月 公司名 利率
2 2000  1      A    a
3 2001  2      A    b
4 2001  3      A    c
5 2000  1      B    d
7 2000  2      B    e
以上是根據(jù)你的數(shù)據(jù)得到的，R中默認(rèn)的是fromLast=FALSE,即若樣本點(diǎn)重復(fù)出現(xiàn)薇缅，則取首次出現(xiàn)的危彩；
否則取最后一次出現(xiàn)的。
#duplicated主要是判定向量或數(shù)據(jù)框中的元素是否重復(fù)泳桦，它返回一個(gè)元素（行）是不是重復(fù)的邏輯向量
> data.set
   Ensembl.Gene.ID Gene.Biotype Chromosome.Name Gene.Start..bp. Gene.End..bp.
1  ENSG00000236666    antisense              22        16274560      16278602
2  ENSG00000236666    antisense              22        16274560      16278602
3  ENSG00000234381   pseudogene              22        16333633      16342783
4  ENSG00000234381   pseudogene              22        16333633      16342783
5  ENSG00000234381   pseudogene              22        16333633      16342783
6  ENSG00000234381   pseudogene              22        16333633      16342783
7  ENSG00000234381   pseudogene              22        16333633      16342783
8  ENSG00000234381   pseudogene              22        16333633      16342783
9  ENSG00000234381   pseudogene              22        16333633      16342783
10 ENSG00000224435   pseudogene              22        16345912      16355362

#構(gòu)建一個(gè)布爾向量汤徽，索引
> index<-duplicated(data.set$Ensembl.Gene.ID)
> index
 [1] FALSE  TRUE FALSE  TRUE  TRUE  TRUE  TRUE  TRUE  TRUE FALSE

#篩選數(shù)據(jù)
> data.set2<-data.set[!index,]  #選中了非重復(fù)的數(shù)據(jù)

• paste和paste0函數(shù)舉例

#R連接函數(shù)paste和paste0
#paste()與paste0()不僅可以連接多個(gè)字符串，還可以將對(duì)象自動(dòng)轉(zhuǎn)換為字符串再相連灸撰，另外還能處理向量
paste("fitbit", month, ".jpg", sep="")
#這個(gè)函數(shù)的特殊地方在于默認(rèn)的分隔符是空格谒府，所以必須指定sep=""，這樣如果month=10時(shí)浮毯，就會(huì)生成fitbit10.jpg這樣的字符串
#要生成12個(gè)月的fitbit文件名
paste("fitbit", 1:12, ".jpg", sep = "")
[1] "fitbit1.jpg"  "fitbit2.jpg"  "fitbit3.jpg"  "fitbit4.jpg"  "fitbit5.jpg"  "fitbit6.jpg"  "fitbit7.jpg"
[8] "fitbit8.jpg"  "fitbit9.jpg"  "fitbit10.jpg" "fitbit11.jpg" "fitbit12.jpg"
#可以看出參數(shù)里面有向量時(shí)的捉對(duì)拼接的效果,如果某個(gè)向量較短,就自動(dòng)補(bǔ)齊
a <- c("甲","乙","丙","丁","戊","己","庚","辛","壬","癸")
b <- c("子","丑","寅","卯","辰","巳","午","未","申","酉","戌","亥")
paste0(a, b)
[1] "甲子" "乙丑" "丙寅" "丁卯" "戊辰" "己巳" "庚午" "辛未" "壬申" "癸酉" "甲戌" "乙亥"

#collapse是合并成一個(gè)字符串時(shí)的分隔符
paste("fitbit", 1:3, ".jpg", sep = "", collapse = "; ")
[1] "fitbit1.jpg; fitbit2.jpg; fitbit3.jpg"

20.對(duì)T檢驗(yàn)結(jié)果的P值和limma包差異分析的P值畫散點(diǎn)圖完疫，看看哪些基因相差很大?

head(nrDEG)
head(DEG_t.test)
DEG_t.test=DEG_t.test[rownames(nrDEG),]
plot(DEG_t.test[,3],nrDEG[,1]) ## 可以看到logFC是相反的
plot(DEG_t.test[,4],nrDEG[,4]) # 可以看到使用limma包和t.test本身的p值差異尚可接受
plot(-log10(DEG_t.test[,4]),-log10(nrDEG[,4]))

重溫下把我完全弄迷糊的正確的nrDEG和正確的DEG_t.test

• 最后畫出的三圖是這樣的

image

• 不然，按照我原來的思路债蓝，第二張圖的畫風(fēng)是這樣的壳鹤，所以說老大一直說的，你要時(shí)刻知道你現(xiàn)在操作的對(duì)象是什么饰迹，就是赤裸裸的真理啊芳誓，不然會(huì)激起密集恐懼癥的余舶，哈哈。

image

最后锹淌，還有個(gè)疑問匿值，那就是第18題中的疑問，當(dāng)時(shí)看別人的代碼赂摆，把exprSet給改成了cbind后的代碼挟憔，后來向老大要來老大代碼，其實(shí)就在R-20題里呢烟号，沒有找到曲楚。發(fā)現(xiàn)那個(gè)地方是和老大一樣的，所以那18題里的第8行代碼中用一開始的未經(jīng)cbind的正常順序樣本的矩陣褥符，那為什么要有dat = cbind(dat1, dat2)這步呢？

最后再仔細(xì)想想老大的代碼抚垃，截圖如下：

image

我想這個(gè)cbind后的dat就是合起來看一下喷楣，了解一下cbind的用法，而我確實(shí)也是又去查了cbind的用法（想知道它是怎么用的鹤树，和merge的區(qū)別铣焊，merge是用by的，通過相同的一列來疊加（如by probe_id）而cbind是行數(shù)相同橫向追加）罕伯，從圖中粉色圈出的地方可以看的其實(shí)后面用的dat1和dat2曲伊。所以我覺得是就可以這樣解釋了！疑問解決了追他。開心坟募。

總結(jié)，對(duì)象很重要邑狸，你要時(shí)刻知道你現(xiàn)在面對(duì)的對(duì)象是誰懈糯！