盤一盤 生信技能樹R語言小作業(yè)（高級）

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

1.安裝一些R包：

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

# 設(shè)置鏡像
options(repos<- c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/") )
source("https://bioconductor.org/biocLite.R")
options("BioC_mirror"<- "https://mirrors.ustc.edu.cn/bioc/")
# 安裝必須的包
install.packages("RSQLite")
#哎呀具滴，出現(xiàn)錯誤蚂维，試一下這個
install.packages('RSQLite', dependencies=TRUE, repos='http://cran.rstudio.com/')
library(DESeq2)
library(reshape2)

2.了解ExpressionSet對象，比如CLL包里面就有data(sCLLex) 力图，找到它包含的元素，提取其表達矩陣(使用exprs函數(shù))室抽，查看其大小
參考：http://www.bio-info-trainee.com/bioconductor_China/software/limma.html
參考：https://github.com/bioconductor-china/basic/blob/master/ExpressionSet.md

# 以下在小作業(yè)中級均提到過
# 參考 http://www.reibang.com/p/e15ee2cd3174
library(CLL)
data(sCLLex)
sCLLex
exprSet <- exprs(sCLLex)
dim(exprSet)

##sCLLex是依賴于CLL這個package的一個對象
samples=sampleNames(sCLLex)
pdata=pData(sCLLex)
# 生成分組信息
group_list=as.character(pdata[,2])
dim(exprSet)
exprSet[1:5,1:5]

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

# 看一下exprSet的結(jié)構(gòu)
str(exprSet)
# 顯示前六行（默認）
head(exprSet)
# 如何想看幾行就看幾行呢 10為例
head(exprSet,n=10)
# 萬物皆可"help"&"?"
help()

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

# BiocManager::install('hgu95av2.db')
library(hgu95av2.db)
# 如前
ls("package:hgu95av2.db")

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

ids <- toTable(hgu95av2SYMBOL)
head(ids)
# 找到 TP53 基因?qū)?yīng)的探針I(yè)D坪圾，也可直接去數(shù)據(jù)框搜
ID <- ids[ids$symbol%in%"TP53",][,1]
# 保存下數(shù)據(jù)晓折，方便下次繼續(xù)，避免從頭開始兽泄。
save(ids,exprSet,pdata,file = 'input.Rdata')

6.理解探針與基因的對應(yīng)關(guān)系漓概，總共多少個基因，基因最多對應(yīng)多少個探針病梢，是哪些基因胃珍，是不是因為這些基因很長，所以在其上面設(shè)計多個探針呢蜓陌？

# dim一下觅彰，發(fā)現(xiàn)11460行，2列
dim(ids)
# unique可見钮热，總共8585個基因
length(unique(ids$symbol))
# 給多個基因排序填抬，顯示最靠后的6個，同樣可以n=10
tail(sort(table(ids$symbol)))
# 查看各出現(xiàn)頻數(shù)的分布隧期，其中6555個基因只出現(xiàn)了一次
table(sort(table(ids$symbol)))
# 可視化
plot(table(sort(table(ids$symbol))))

7.第二步提取到的表達矩陣是12625個探針在22個樣本的表達量矩陣飒责，找到那些不在 hgu95av2.db 包收錄的對應(yīng)著SYMBOL的探針。
提示：有1165個探針是沒有對應(yīng)基因名字的仆潮。

# 查看多少%in%宏蛉，多少!%in%,即分布⌒灾茫可知拾并，11460個TRUE，1165個FALSE
table(rownames(exprSet) %in% ids$probe_id)
dim(exprSet)

8.過濾表達矩陣蚌讼，刪除那1165個沒有對應(yīng)基因名字的探針辟灰。

# 篩選
exprSet=exprSet[rownames(exprSet) %in% ids$probe_id,]
# 11460行
dim(exprSet)

9.整合表達矩陣，多個探針對應(yīng)一個基因的情況下篡石，只保留在所有樣本里面平均表達量最大的那個探針芥喇。
提示，理解 tapply,by,aggregate,split 函數(shù) , 首先對每個基因找到最大表達量的探針凰萨。
然后根據(jù)得到探針去過濾原始表達矩陣

ids=ids[match(rownames(exprSet),ids$probe_id),]
head(ids)
exprSet[1:5,1:5]
if(F){
    # 定義最大平均表達量的...
  tmp = by(exprSet,ids$symbol,
           function(x) rownames(x)[which.max(rowMeans(x))] )
  probes = as.character(tmp)
  
  dim(exprSet)
  #篩選表達矩陣
  exprSet=exprSet[rownames(exprSet) %in% probes ,]
  # 看下篩選之后行情況继控，8585
  dim(exprSet)
  # 表達矩陣探針換為基因名
  rownames(exprSet)=ids[match(rownames(exprSet),ids$probe_id),2]
  exprSet[1:5,1:5]
}         

# 兩者一樣嗎，F(xiàn)ALSE
identical(ids$probe_id,rownames(exprSet))
# match一下胖眷，再鑒定一遍武通，TRUE
ids=ids[match(rownames(exprSet),ids$probe_id),]
identical(ids$probe_id,rownames(exprSet))
#新建dat
dat=exprSet
#ids新建median這一列，列名為median珊搀，同時對dat這個矩陣按行操作冶忱，取每一行的中位數(shù)，將結(jié)果給到median這一列的每一行
ids$median=apply(dat,1,median)
#對ids$symbol按照ids$median中位數(shù)從大到小排列的順序排序境析，將對應(yīng)的行賦值為一個新的ids（跟隨改變）
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
#將symbol這一列取取出重復(fù)項囚枪，'!'為否派诬，即取出不重復(fù)的項，去除重復(fù)的gene 链沼，保留每個基因最大表達量結(jié)果
ids=ids[!duplicated(ids$symbol),]
dim(ids)
#新的ids取出probe_id這一列默赂，將dat按照取出的這一列中的每一行組成一個新的dat
dat=dat[ids$probe_id,] 
#把ids的symbol這一列中的每一行給dat作為dat的行名
rownames(dat)=ids$symbol
#保留每個基因ID第一次出現(xiàn)的信息
dat[1:4,1:4] 
dim(dat)

10.把過濾后的表達矩陣更改行名為基因的symbol，因為這個時候探針和基因是一對一關(guān)系了括勺。見上

#把ids的symbol這一列中的每一行給dat作為dat的行名
rownames(dat)=ids$symbol
#保留每個基因ID第一次出現(xiàn)的信息
dat[1:4,1:4] 
dim(dat)

11.對第10步得到的表達矩陣進行探索缆八，先畫第一個樣本的所有基因的表達量的boxplot,hist,density ， 然后畫所有樣本的 這些圖
參考：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html
理解ggplot2的繪圖語法疾捍，數(shù)據(jù)和圖形元素的映射關(guān)系

# 管道回來
exprSet=dat
# 
exprSet['GAPDH',]
boxplot(exprSet[,1])
boxplot(exprSet['GAPDH',])
exprSet['ACTB',]
# 用reshape2包畫
library(reshape2)
# 整理分組矩陣
exprSet_L=melt(exprSet)
colnames(exprSet_L)=c('probe','sample','value')
# 獲得分組信息
group_list=as.character(pdata[,2])
exprSet_L$group=rep(group_list,each=nrow(exprSet))
head(exprSet_L)
### ggplot2畫圖 
library(ggplot2)
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)
p=ggplot(exprSet_L,aes(value,col=group))+geom_density() 
print(p)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun.y="mean",geom="point",shape=23,size=3,fill="red")
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face='bold'),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print(p)

放最后一張圖吧


image.png

12.理解統(tǒng)計學(xué)指標(biāo)mean,median,max,min,sd,var,mad并計算出每個基因在所有樣本的這些統(tǒng)計學(xué)指標(biāo)奈辰，最后按照mad值排序，取top 50 mad值的基因乱豆，得到列表冯挎。

# 異曲同工
g_mean <- tail(sort(apply(exprSet,1,mean)),50)
g_median <- tail(sort(apply(exprSet,1,median)),50)
g_max <- tail(sort(apply(exprSet,1,max)),50)
g_min <- tail(sort(apply(exprSet,1,min)),50)
g_sd <- tail(sort(apply(exprSet,1,sd)),50)
g_var <- tail(sort(apply(exprSet,1,var)),50)
g_mad <- tail(sort(apply(exprSet,1,mad)),50)
g_mad
names(g_mad)
image.png

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

# 熱一下試試
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.png

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

## UpSetR
# 看下說明書 https://cran.r-project.org/web/packages/UpSetR/README.html
library(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),
               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)
image.png

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

pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
group_list
dim(exprSet)
exprSet[1:5,1:5]

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

## hclust
# 更改表達矩陣列名
colnames(exprSet)=paste(group_list,1:22,sep='')
# 定義nodePar
nodePar <- list(lab.cex = 0.6, pch = c(NA, 19), 
                cex = 0.7, col = "blue")
# 聚類
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10)) 
# 繪圖
plot(as.dendrogram(hc), nodePar = nodePar,  horiz = TRUE)
image.png

17.對所有樣本的表達矩陣進行PCA分析并且繪圖疲酌，同樣要添加表型信息蜡峰。

library(ggfortify)
# 互換行和列，dim一下
df=as.data.frame(t(exprSet))
dim(df)
# 不要view df朗恳，列太多湿颅，軟件會崩掉；
df$group=group_list 
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = 'group')
#或者
#畫主成分分析圖需要加載這兩個包
library("FactoMineR")
library("factoextra") 
df=as.data.frame(t(exprSet))
dat.pca <- PCA(df, graph = FALSE)#現(xiàn)在dat最后一列是group_list粥诫，需要重新賦值給一個dat.pca,這個矩陣是不含有分組信息的
fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind = group_list, # color by groups
             # palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
)
image.png

image.png

18.根據(jù)表達矩陣及樣本分組信息進行批量T檢驗油航，得到檢驗結(jié)果表格

dat = exprSet
group_list=as.factor(group_list)
group1 = which(group_list == levels(group_list)[1])
group2 = which(group_list == levels(group_list)[2])
dat1 = dat[, group1]
dat2 = dat[, group2]
dat = cbind(dat1, dat2)
pvals = apply(exprSet, 1, function(x){
  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檢驗結(jié)果表格，包含log2FC怀浆、pvals和p.adj等谊囚，通常認為t<0.05即有統(tǒng)計學(xué)意義
head(DEG_t.test)
image.png

19.使用limma包對表達矩陣及樣本分組信息進行差異分析，得到差異分析表格执赡，重點看logFC和P值镰踏，畫個火山圖(就是logFC和-log10(P值)的散點圖)。

中級作業(yè)已經(jīng)涉及沙合。
參考 http://www.reibang.com/p/e15ee2cd3174

# DEG by limma 
suppressMessages(library(limma)) 
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design

## 下面的 contrast.matrix 矩陣非常重要奠伪，制定了誰比誰這個規(guī)則
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix 
##這個矩陣聲明，我們要把progres.組跟stable進行差異分析比較
##step1
fit <- lmFit(exprSet,design)
##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) 
#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)
head(nrDEG)
image.png
## 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',])
)
this_tile
head(DEG)
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)
image.png

20.對T檢驗結(jié)果的P值和limma包差異分析的P值畫散點圖绊率，看看哪些基因相差很大含末。

### different P values 
head(nrDEG)
head(DEG_t.test)
# 將limma生成的nrDEG與t檢驗合并
DEG_t.test=DEG_t.test[rownames(nrDEG),]
## 可以看到logFC是相反的
plot(DEG_t.test[,3],nrDEG[,1]) 
# 可以看到使用limma包和t.test本身的p值差異尚可接受
plot(DEG_t.test[,4],nrDEG[,4]) 
plot(-log10(DEG_t.test[,4]),-log10(nrDEG[,4]))
image.png

image.png

image.png
# 找3個基因看一下
exprSet['GAPDH',]
exprSet['ACTB',]
exprSet['DLEU1',]

library(ggplot2)
library(ggpubr)
my_comparisons <- list(
  c("stable", "progres.")
)
dat=data.frame(group=group_list,
               sampleID= names(exprSet['DLEU1',]),
               values= as.numeric(exprSet['DLEU1',]))
ggboxplot(
  dat, x = "group", y = "values",
  color = "group",
  add = "jitter"
)+
  stat_compare_means(comparisons = my_comparisons, method = "t.test")
image.png
## heatmap 
library(pheatmap)
choose_gene=head(rownames(nrDEG),25)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)
image.png

后面幾個題目會用函數(shù)就行，搞清楚根據(jù)不同的數(shù)據(jù)如何修改即舌。

參考來源：生信技能樹

友情鏈接：

課程分享
生信技能樹全球公益巡講
（https://mp.weixin.qq.com/s/E9ykuIbc-2Ja9HOY0bn_6g）
B站公益74小時生信工程師教學(xué)視頻合輯
（https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw）
招學(xué)徒：
（https://mp.weixin.qq.com/s/KgbilzXnFjbKKunuw7NVfw）

歡迎關(guān)注公眾號：青島生信菜鳥團

最后編輯于 ：
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者
  • 人面猴
    序言：七十年代末，一起剝皮案震驚了整個濱河市挎袜，隨后出現(xiàn)的幾起案子顽聂，更是在濱河造成了極大的恐慌，老刑警劉巖盯仪，帶你破解...
    沈念sama閱讀 221,273評論 6贊 515
  • 死咒
    序言：濱河連續(xù)發(fā)生了三起死亡事件紊搪，死亡現(xiàn)場離奇詭異，居然都是意外死亡全景，警方通過查閱死者的電腦和手機耀石，發(fā)現(xiàn)死者居然都...
    沈念sama閱讀 94,349評論 3贊 398
  • 救了他兩次的神仙讓他今天三更去死
    文/潘曉璐 我一進店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來爸黄，“玉大人滞伟，你說我怎么就攤上這事】还螅” “怎么了梆奈？”我有些...
    開封第一講書人閱讀 167,709評論 0贊 360
  • 道士緝兇錄：失蹤的賣姜人
    文/不壞的土叔 我叫張陵，是天一觀的道長称开。 經(jīng)常有香客問我亩钟，道長，這世上最難降的妖魔是什么鳖轰？ 我笑而不...
    開封第一講書人閱讀 59,520評論 1贊 296
  • ?港島之戀（遺憾婚禮）
    正文 為了忘掉前任清酥，我火速辦了婚禮，結(jié)果婚禮上蕴侣，老公的妹妹穿的比我還像新娘焰轻。我一直安慰自己，他們只是感情好昆雀，可當(dāng)我...
    茶點故事閱讀 68,515評論 6贊 397
  • 惡毒庶女頂嫁案：這布局不是一般人想出來的
    文/花漫 我一把揭開白布鹦马。 她就那樣靜靜地躺著，像睡著了一般忆肾。 火紅的嫁衣襯著肌膚如雪审轮。 梳的紋絲不亂的頭發(fā)上，一...
    開封第一講書人閱讀 52,158評論 1贊 308
  • 城市分裂傳說
    那天扒最，我揣著相機與錄音赏殃，去河邊找鬼。 笑死，一個胖子當(dāng)著我的面吹牛和悦，可吹牛的內(nèi)容都是我干的退疫。 我是一名探鬼主播，決...
    沈念sama閱讀 40,755評論 3贊 421
  • 雙鴛鴦連環(huán)套：你想象不到人心有多黑
    文/蒼蘭香墨 我猛地睜開眼鸽素，長吁一口氣：“原來是場噩夢啊……” “哼褒繁！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起馍忽，我...
    開封第一講書人閱讀 39,660評論 0贊 276
  • 萬榮殺人案實錄
    序言：老撾萬榮一對情侶失蹤棒坏，失蹤者是張志新（化名）和其女友劉穎，沒想到半個月后遭笋，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體坝冕，經(jīng)...
    沈念sama閱讀 46,203評論 1贊 319
  • ?護林員之死
    正文 獨居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛 以下內(nèi)容為張勛視角 年9月15日...
    茶點故事閱讀 38,287評論 3贊 340
  • ?白月光啟示錄
    正文 我和宋清朗相戀三年瓦呼，在試婚紗的時候發(fā)現(xiàn)自己被綠了喂窟。 大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
    茶點故事閱讀 40,427評論 1贊 352
  • 活死人
    序言：一個原本活蹦亂跳的男人離奇死亡央串，死狀恐怖磨澡，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情质和，我是刑警寧澤钱贯，帶...
    沈念sama閱讀 36,122評論 5贊 349
  • ?日本核電站爆炸內(nèi)幕
    正文 年R本政府宣布，位于F島的核電站侦另，受9級特大地震影響秩命，放射性物質(zhì)發(fā)生泄漏。R本人自食惡果不足惜褒傅，卻給世界環(huán)境...
    茶點故事閱讀 41,801評論 3贊 333
  • 男人毒藥：我在死后第九天來索命
    文/蒙蒙 一弃锐、第九天 我趴在偏房一處隱蔽的房頂上張望。 院中可真熱鬧殿托，春花似錦霹菊、人聲如沸。這莊子的主人今日做“春日...
    開封第一講書人閱讀 32,272評論 0贊 23
  • 一樁弒父案旋廷，背后竟有這般陰謀
    文/蒼蘭香墨 我抬頭看了看天上的太陽。三九已至礼搁，卻和暖如春饶碘，著一層夾襖步出監(jiān)牢的瞬間，已是汗流浹背馒吴。 一陣腳步聲響...
    開封第一講書人閱讀 33,393評論 1贊 272
  • 情欲美人皮
    我被黑心中介騙來泰國打工扎运， 沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留瑟曲，地道東北人。 一個月前我還...
    沈念sama閱讀 48,808評論 3贊 376
  • 代替公主和親
    正文 我出身青樓豪治，卻偏偏與公主長得像洞拨，于是被迫代替她去往敵國和親。 傳聞我的和親對象是個殘疾皇子负拟，可洞房花燭夜當(dāng)晚...
    茶點故事閱讀 45,440評論 2贊 359