首先友情宣傳生信技能樹
- 全國巡講:R基礎(chǔ)，Linux基礎(chǔ)和RNA-seq實(shí)戰(zhàn)演練 : 預(yù)告：12月28-30長(zhǎng)沙站
- 廣州珠江新城GEO數(shù)據(jù)挖掘滾動(dòng)開班
題目來源：http://www.bio-info-trainee.com/3409.html
- 安裝一些R包：
數(shù)據(jù)包： ALL, CLL, pasilla, airway
軟件包：limma，DESeq2停巷，clusterProfiler
工具包：reshape2
繪圖包：ggplot2
不同領(lǐng)域的R包使用頻率不一樣瓮顽，在生物信息學(xué)領(lǐng)域缅刽，尤其需要掌握bioconductor系列包藤滥。
rm(list=ls())
options(stringsAsFactors = F)
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")}
if (!require("ggplot2")) {BiocManager::install("ggplot2")}
- 了解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
help("ExpressionSet")
library("CLL")
data(sCLLex)
exprSet <- exprs(sCLLex)
dim(exprSet)
#[1] 12625 22
head(rownames(exprSet))
#[1] "1000_at" "1001_at" "1002_f_at" "1003_s_at" "1004_at"
#[6] "1005_at"
colnames(exprSet)
# [1] "CLL11.CEL" "CLL12.CEL" "CLL13.CEL" "CLL14.CEL" "CLL15.CEL"
# [6] "CLL16.CEL" "CLL17.CEL" "CLL18.CEL" "CLL19.CEL" "CLL20.CEL"
#[11] "CLL21.CEL" "CLL22.CEL" "CLL23.CEL" "CLL24.CEL" "CLL2.CEL"
#[16] "CLL3.CEL" "CLL4.CEL" "CLL5.CEL" "CLL6.CEL" "CLL7.CEL"
#[21] "CLL8.CEL" "CLL9.CEL"
sCLLex
#ExpressionSet (storageMode: lockedEnvironment)
#assayData: 12625 features, 22 samples
# element names: exprs
#protocolData: none
#phenoData
# sampleNames: CLL11.CEL CLL12.CEL ... CLL9.CEL (22 total)
# varLabels: SampleID Disease
# varMetadata: labelDescription
#featureData: none
#experimentData: use 'experimentData(object)'
#Annotation: hgu95av2
- 了解 str,head,help函數(shù)枢希，作用于第2題提取到的表達(dá)矩陣
str(exprSet)
#num [1:12625, 1:22] 5.74 2.29 3.31 1.09 7.54 ...
# - attr(*, "dimnames")=List of 2
# ..$ : chr [1:12625] "1000_at" "1001_at" "1002_f_at" "1003_s_at" ...
# ..$ : chr [1:22] "CLL11.CEL" "CLL12.CEL" "CLL13.CEL" "CLL14.CEL" ...
head(exprSet)
help(exprSet)
help
- 安裝并了解
hgu95av2.db包,看看ls("package:hgu95av2.db")后顯示的那些變量
if (!require("hgu95av2.db")){BiocManager::install("hgu95av2.db")}
ls("package:hgu95av2.db")
# [1] "hgu95av2" "hgu95av2.db"
# [3] "hgu95av2_dbconn" "hgu95av2_dbfile"
# [5] "hgu95av2_dbInfo" "hgu95av2_dbschema"
# [7] "hgu95av2ACCNUM" "hgu95av2ALIAS2PROBE"
# [9] "hgu95av2CHR" "hgu95av2CHRLENGTHS"
#[11] "hgu95av2CHRLOC" "hgu95av2CHRLOCEND"
#[13] "hgu95av2ENSEMBL" "hgu95av2ENSEMBL2PROBE"
#[15] "hgu95av2ENTREZID" "hgu95av2ENZYME"
#[17] "hgu95av2ENZYME2PROBE" "hgu95av2GENENAME"
#[19] "hgu95av2GO" "hgu95av2GO2ALLPROBES"
#[21] "hgu95av2GO2PROBE" "hgu95av2MAP"
#[23] "hgu95av2MAPCOUNTS" "hgu95av2OMIM"
#[25] "hgu95av2ORGANISM" "hgu95av2ORGPKG"
#[27] "hgu95av2PATH" "hgu95av2PATH2PROBE"
#[29] "hgu95av2PFAM" "hgu95av2PMID"
#[31] "hgu95av2PMID2PROBE" "hgu95av2PROSITE"
#[33] "hgu95av2REFSEQ" "hgu95av2SYMBOL"
#[35] "hgu95av2UNIGENE" "hgu95av2UNIPROT"
- 理解
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[IDs$symbol=='TP53',]
# probe_id symbol
#966 1939_at TP53
#997 1974_s_at TP53
#1420 31618_at TP53
- 理解探針與基因的對(duì)應(yīng)關(guān)系，總共多少個(gè)基因苞轿，基因最多對(duì)應(yīng)多少個(gè)探針茅诱，是哪些基因，是不是因?yàn)檫@些基因很長(zhǎng)搬卒，所以在其上面設(shè)計(jì)多個(gè)探針呢瑟俭？
length(unique(IDs$symbol))
#[1] 8584
#總共有8584個(gè)基因
max(table(IDs$symbol))
#[1] 8
#一個(gè)基因最多對(duì)應(yīng)8個(gè)探針
table(table(IDs$symbol)==8)
#FALSE TRUE
# 8579 5
#有5個(gè)同時(shí)對(duì)應(yīng)8個(gè)探針的基因
a <- table(IDs$symbol)
for (i in 1:length(a)){
if (a[i]==8) {print(names(a[i]))}
}
[1] "GAPDH"
[1] "INPP4A"
[1] "MYB"
[1] "PTGER3"
[1] "STAT1"
- 第2題提取到的表達(dá)矩陣是12625個(gè)探針在22個(gè)樣本的表達(dá)量矩陣，找到那些不在
hgu95av2.db包收錄的對(duì)應(yīng)著SYMBOL的探針契邀。提示：有1165個(gè)探針是沒有對(duì)應(yīng)基因名字的摆寄。
rownames(exprSet)[!(rownames(exprSet) %in% IDs$probe_id)]
library(dplyr)
rownames(exprSet)[!(rownames(exprSet) %in% IDs$probe_id)] %>% length
#[1] 1166
rownames(exprSet)[!(rownames(exprSet) %in% IDs$probe_id)] %>% head
#[1] "1007_s_at" "1047_s_at" "1089_i_at" "108_g_at" "1090_f_at"
#[6] "1099_s_at"
- 過濾表達(dá)矩陣，刪除那1165個(gè)沒有對(duì)應(yīng)基因名字的探針坯门。
tmp <-rownames(exprSet) %in% IDs$probe_id
head(tmp)
#[1] TRUE TRUE TRUE TRUE TRUE TRUE
new_eSet <- exprSet[tmp,]
nrow(exprSet)
#[1] 12625
nrow(new_eSet)
#[1] 11459
- 整合表達(dá)矩陣微饥，多個(gè)探針對(duì)應(yīng)一個(gè)基因的情況下，只保留在所有樣本里面平均表達(dá)量最大的那個(gè)探針古戴。
提示：
1.理解tapply欠橘，by，aggregate允瞧，split函數(shù) , 首先對(duì)每個(gè)基因找到最大表達(dá)量的探針简软。
2.然后根據(jù)得到探針去過濾原始表達(dá)矩陣
IDs$median <- apply(new_eSet,1,median)
nrow(IDs)
#[1] 11459
IDs <- IDs[order(IDs$symbol,IDs$median,decreasing = T),]
IDs <- IDs[!duplicated(IDs$symbol),]
nrow(IDs)
#[1] 8584
genes_expr <- new_eSet[as.character(IDs$probe_id),]
dim(genes_expr)
#[1] 8584 22
genes_expr[1:4,1:4]
# CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL
#36456_at 6.645791 7.350613 6.333290 6.603640
#38626_at 5.289264 6.677600 4.447104 7.008260
#36958_at 3.949769 5.423343 3.540189 5.234420
#35995_at 4.316881 2.705329 3.131087 2.821306
- 把過濾后的表達(dá)矩陣更改行名為基因的symbol，因?yàn)檫@個(gè)時(shí)候探針和基因是一對(duì)一關(guān)系了述暂。
rownames(genes_expr) <- IDs$symbol
genes_expr <- genes_expr[order(rownames(genes_expr)),]
genes_expr[1:4,1:4]
# CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL
#AADAC 2.837647 2.893664 2.848301 2.868344
#AAK1 3.755016 3.807378 3.716134 3.759093
#AAMP 1.953082 2.010410 2.244479 1.498089
#AANAT 1.179319 1.202284 1.180669 1.211086
- 對(duì)第10題得到的表達(dá)矩陣進(jìn)行探索痹升，先畫第一個(gè)樣本的所有基因的表達(dá)量的boxplot,hist,density ， 然后畫所有樣本的這些圖
1.參考：http://bio-info-trainee.com/tmp/basic_visualization_for_expression_matrix.html
2.理解ggplot2的繪圖語法畦韭，數(shù)據(jù)和圖形元素的映射關(guān)系
library(ggplot2)
library(reshape2)
library(patchwork)
#讀入表型信息
phenoData <- pData(sCLLex)
colnames(genes_expr) <- phenoData[colnames(genes_expr),1]
#表達(dá)矩陣轉(zhuǎn)置為長(zhǎng)數(shù)據(jù)框
expr_L <- melt(genes_expr)
colnames(expr_L) <- c("gene","sample","exp")
expr_L$pheno <- as.character(phenoData[expr_L$sample,2]) # 長(zhǎng)數(shù)據(jù)框添加表型信息
expr_L[,1] <- as.character(expr_L[,1])
expr_L[,2] <- as.character(expr_L[,2])
#把第一個(gè)樣本的所有數(shù)據(jù)提取出來作為一個(gè)新的數(shù)據(jù)框expr_L_1
expr_L_1 <- expr_L[expr_L$sample==expr_L$sample[1],]
#箱線圖
(boxplot_1 <- ggplot(data = expr_L_1,aes(x=sample,y=exp)) +
geom_boxplot())
(boxplot_all <- ggplot(data = expr_L,aes(x=sample,y=exp,fill=pheno)) +
geom_boxplot() +
theme(axis.text.x = element_text(angle = 90, hjust = 1)))
draw_boxplot <- boxplot_1 / boxplot_all
ggsave('boxplot.PNG',draw_boxplot)
#頻數(shù)圖
(histogram_1 <- ggplot(data = expr_L_1,aes(exp,fill=pheno)) +
geom_histogram(bins=200))
(histogram_all <- ggplot(data = expr_L,aes(exp,fill=pheno)) +
geom_histogram() +
facet_wrap(~sample,nrow = 4))
ggsave('histogram_all.png',histogram_all)
#密度圖
(density_1 <- ggplot(data = expr_L_1,aes(exp,fill=pheno)) +
geom_density())
(density_all <- ggplot(data = expr_L,aes(exp,fill=pheno)) +
geom_density() +
facet_wrap(~sample,nrow = 4))
ggsave("density_all.png",density_all)
箱線圖
頻數(shù)圖
密度圖
- 理解統(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ì)看砖织。
e_mean <- apply(genes_expr,1,mean)
e_median <- apply(genes_expr,1,median)
e_max <- apply(genes_expr,1,max)
e_min <- apply(genes_expr,1,min)
e_sd <- apply(genes_expr,1,sd)
e_var <- apply(genes_expr,1,var)
e_mad <- apply(genes_expr,1,mad)
top50_mad <- head(sort(e_mad,decreasing = T),50)
head(top50_mad)
# FAM30A IGF2BP3 DMD TCF7 SLAMF1 FOS
#3.982191 3.234011 3.071541 2.993352 2.944105 2.938078
names(top50_mad)
# [1] "FAM30A" "IGF2BP3" "DMD" "TCF7" "SLAMF1" "FOS"
# [7] "LGALS1" "IGLC1" "ZAP70" "FCN1" "LHFPL2" "HBB"
#[13] "S100A8" "GUSBP11" "COBLL1" "VIPR1" "PCDH9" "IGH"
#[19] "ZNF804A" "TRIB2" "OAS1" "CCL3" "GNLY" "CYBB"
#[25] "VAMP5" "RNASE6" "RGS2" "PLXNC1" "CAPG" "RBM38"
#[31] "VCAN" "APBB2" "ARF6" "TGFBI" "NR4A2" "S100A9"
#[37] "ZNF266" "TSPYL2" "CLEC2B" "FLNA" "H1-10" "DUSP5"
#[43] "DUSP6" "ANXA4" "LPL" "THEMIS2" "P2RY14" "ARHGAP44"
#[49] "TNFSF9" "PFN2"
- 根據(jù)第12題得到top 50 mad值的基因列表來取表達(dá)矩陣的子集款侵，并且熱圖可視化子表達(dá)矩陣。試試看其它5種熱圖的包的不同效果侧纯。
library(pheatmap)
top50_mad_matrix <- genes_expr[names(top50_mad),]
top50_mad_matrix <- t(scale(t(top50_mad_matrix)))
phenoData <- pData(sCLLex)
rownames(phenoData) <- phenoData$SampleID
annotation_col <- data.frame(
pheno <- factor(phenoData[colnames(top50_mad_matrix),2])
)
colnames(annotation_col) <- "phenoData"
rownames(annotation_col) <- colnames(top50_mad_matrix)
pheatmap(top50_mad_matrix,
cellwidth = 35, cellheight = 12, fontsize = 8,
annotation_col = annotation_col,
filename = "heatmap.png")
heatmap.png
- 取不同統(tǒng)計(jì)學(xué)指標(biāo)
mean新锈，median，max眶熬，min妹笆，sd，var聋涨，mad的各top50基因列表晾浴，使用UpSetR包來看他們之間的overlap情況。
- 在第2題的基礎(chǔ)上面提取
CLL包里面的data(sCLLex)數(shù)據(jù)對(duì)象的樣本的表型數(shù)據(jù)牍白。
phenoData <- pData(sCLLex)
phenoData
# SampleID Disease
#CLL11.CEL CLL11 progres.
#CLL12.CEL CLL12 stable
#CLL13.CEL CLL13 progres.
#CLL14.CEL CLL14 progres.
#CLL15.CEL CLL15 progres.
#CLL16.CEL CLL16 progres.
#CLL17.CEL CLL17 stable
#CLL18.CEL CLL18 stable
#CLL19.CEL CLL19 progres.
#CLL20.CEL CLL20 stable
#CLL21.CEL CLL21 progres.
#CLL22.CEL CLL22 stable
#CLL23.CEL CLL23 progres.
#CLL24.CEL CLL24 stable
#CLL2.CEL CLL2 stable
#CLL3.CEL CLL3 progres.
#CLL4.CEL CLL4 progres.
#CLL5.CEL CLL5 progres.
#CLL6.CEL CLL6 progres.
#CLL7.CEL CLL7 progres.
#CLL8.CEL CLL8 progres.
#CLL9.CEL CLL9 stable
- 對(duì)所有樣本的表達(dá)矩陣進(jìn)行聚類并且繪圖脊凰，然后添加樣本的臨床表型數(shù)據(jù)信息(更改樣本名)
- 對(duì)所有樣本的表達(dá)矩陣進(jìn)行PCA分析并且繪圖，同樣要添加表型信息茂腥。
- 根據(jù)表達(dá)矩陣及樣本分組信息進(jìn)行批量T檢驗(yàn)狸涌，得到檢驗(yàn)結(jié)果表格
- 使用
limma包對(duì)表達(dá)矩陣及樣本分組信息進(jìn)行差異分析，得到差異分析表格最岗，重點(diǎn)看logFC和P值帕胆，畫個(gè)火山圖(就是logFC和-log10(P值)的散點(diǎn)圖。)般渡。
- 對(duì)T檢驗(yàn)結(jié)果的P值和limma包差異分析的P值畫散點(diǎn)圖懒豹，看看哪些基因相差很大芙盘。
