1 介紹R及R Studio
右上角history，To Console選擇到左下角，To Source選擇到左上編輯器內(nèi)。
畫圖沒出圖：dev.off
包的安裝路徑：.libPaths()
善用幫助文檔：?函數(shù)
2 R語言基礎(chǔ)變量
變量結(jié)構(gòu)：向量、矩陣、數(shù)組、數(shù)據(jù)框、列表
當創(chuàng)建向量元素出現(xiàn)字符，那么全部元素都是字符型。
重要函數(shù)：class( ) str( ) 報錯的時候查看元素性質(zhì)
is. 函數(shù)：判斷
as. 函數(shù)：轉(zhuǎn)換成
索引：取數(shù)據(jù) b[,3] 等價于 b[,F,F,T,F]
直接搜索：grep('RNA-Seq', a$Assay_type)、grepl() 結(jié)果判斷為TRUE FALSE、
table() 計數(shù)
3 外部數(shù)據(jù)導入導出
把txt文件導入R
read.table('SraRunTable.txt')
R無法識別正確的列數(shù)
增加參數(shù)：
a <- read.table('SraRunTable.txt', sep = '\t')
繪制表頭：
a <- read.table('SraRunTable.txt', header = TRUE,
sep = '\t')
image.png
更復雜的文件導入方式：!：有感嘆號開頭的內(nèi)容不要
b <- read.table('R_bilibili-master/R語言小作業(yè)/GSE17215_series_matrix.txt.gz',
comment.char = "!",
header = TRUE,
sep = '\t')
把剛才文件導出成csv文件
write.csv(b,'GSE17215_series_matrix.csv')
把第一列名去掉
參數(shù)：row.name=FALSE
讀取文件：
write.table(b,'tmp.csv',sep=',')
把第一列添加到列名，且刪除原本第一列
rownames(b) <- b[,1]
b <- b[,-1]
在配置好R包后，繪制圖像：
pheatmap::pheatmap(b[1:10,])
image.png
b數(shù)值太大了，取對數(shù)
b <- log2(b)
pheatmap::pheatmap(b[1:10,])
image.png
導出數(shù)據(jù)最好保存成Rdata文件
save(b, file = 'b_input.Rdata')
讀取
load(file='b_input.Rdata')
4 中級變量操作
a <- read.table('SraRunTable.txt', header = TRUE,
sep = '\t')
以a為例子。
選取最大值和最小值
> sort(a$MBases)[1]
[1] 2128
> sort(a$MBases,decreasing = TRUE)[1]
[1] 19506
也可以直接用函數(shù)
max(a$MBases)
min(a$MBases)
也可以用統(tǒng)計學五分位數(shù)：
> fivenum(a$MBases)
[1] 2128.0 5513.0 6873.5 8415.0 19506.0
統(tǒng)計該變量小于5000的個數(shù)
> table(a$MBases < 5000)
FALSE TRUE
151 31
把小于5000的數(shù)據(jù)新建一個文件(列數(shù)據(jù))
d <- a[a$MBases < 5000,]
發(fā)現(xiàn)吧絕大部分RNA-Seq的數(shù)據(jù)取出來了，我們看一下外顯子和RNA-Seq的數(shù)據(jù)分布
boxplot(a$MBases~a$Assay_Type)
image.png
差異很大，所以要分組
wes <- a[a$Assay_Type=='Wxs',]
rna <- a[a$Assay_Type=='RNA-Seq',]
然后再根據(jù)每個組的數(shù)據(jù)分布來篩選
b <- read.table('R_bilibili-master/R語言小作業(yè)/GSE17215_series_matrix.txt.gz',
comment.char = "!",
header = TRUE,
sep = '\t')
rownames(b) <- b[,1]
b <- b[,-1]
b <- log2(b)
以b為例子。
計算b第一行的數(shù)據(jù)，但是查看格式
> str(b[1,])
'data.frame': 1 obs. of 6 variables:
$ GSM431121: num 8.91
$ GSM431122: num 9.22
$ GSM431123: num 11.4
$ GSM431124: num 11.3
$ GSM431125: num 11.4
$ GSM431126: num 11.4
是數(shù)據(jù)框格式，要轉(zhuǎn)變?yōu)閿?shù)字
> as.numeric(b[1,])
[1] 8.911691 9.221081 11.410364 11.325483 11.418782 11.360438
就可以進行取平均值操作
> mean(as.numeric(b[1,]))
[1] 10.60797
想要批量計算均值？采用函數(shù)（相對最簡單）
> head(rowMeans(b))
1007_s_at 1053_at 117_at 121_at 1255_g_at 1294_at
10.607973 7.925899 5.193894 7.168633 4.275652 5.686036
也可以建立循環(huán)，類似C語言的for循環(huán)
for (i in 1:nrow(b)) {
print(mean(as.numeric(b[1,])))
}
或者
apply(b, 1, function(x){
x <- as.numeric(b[1,])
mean(x)
})
自定義函數(shù)：選取每一列最大值
rowMax=function(x){
apply(x,1,max)
}
計算每一行方差值
apply(b,1,sd)
按照遞減順序，取前50，并拿到基因名字
cg=names(sort(apply(b,1,sd),decreasing=TRUE)[1:50])
可以隨機取50個畫熱圖
pheatmap::pheatmap(b[sample(1:nrow(b),50),])
用cg來畫熱圖
pheatmap::pheatmap(b[cg,])
image.png
5 熱圖
生成向量，定義成矩陣進行畫圖
a1 <- rnorm(100)
dim(a1) <- c(5,20)
a2 <- rnorm(100)+2
dim(a2) <- c(5,20)
library(pheatmap)
pheatmap(a1,cluster_rows = FALSE,show_colnames = FALSE)
pheatmap(cbind(a1,a2))
pheatmap(cbind(a1,a2),show_rownames = FALSE, show_colnames = FALSE)
不排序
pheatmap(cbind(a1,a2),cluster_cols = F)
可以看到a2整體比a1大（左側(cè)藍a1偏小，右側(cè)紅a2偏大）
image.png
學習函數(shù)paste
> paste('a1',1:20,sep = '_')
[1] "a1_1" "a1_2" "a1_3" "a1_4" "a1_5" "a1_6" "a1_7" "a1_8" "a1_9"
[10] "a1_10" "a1_11" "a1_12" "a1_13" "a1_14" "a1_15" "a1_16" "a1_17" "a1_18"
[19] "a1_19" "a1_20"
b <- cbind(a1,a2)
b <- as.data.frame(b)
names(b) <- c(paste('a1',1:20,sep = '_'),paste('a2',1:20,sep = '_'))
pheatmap(b,cluster_cols = F)
image.png
圖就有橫坐標了。
增加group名
b <- cbind(a1,a2)
b <- as.data.frame(b)
pheatmap(cbind(a1,a2),cluster_cols = F)
names(b) <- c(paste('a1',1:20,sep = '_'),paste('a2',1:20,sep = '_'))
tmp <- data.frame(group=c(rep('a1',20), rep('a2',20)))
rownames(tmp) <- colnames(b)
pheatmap(b,annotation_col = tmp)
image.png
根據(jù)help文件里面的example，可以不斷豐富熱圖。
6 選取差異明顯的基因的表達量矩陣繪制熱圖
1轉(zhuǎn)置
sort選取極值，用scale標準化。拉平大值和小值
n[n>2] <- 2
n[n<-2] <- -2
2再繪制
names(tail((sort(apply(dat, 1, sd))),1000)),t(轉(zhuǎn)置)
7 id轉(zhuǎn)換
首先導入ENSG文本
a <- read.table('e1.txt')
ENSG00000000003.13
分割點號：
> strsplit('ENSG00000000003.13', '[.]')
[[1]]
[1] "ENSG00000000003" "13"
取這個向量第一個元素，再取第一個元素
> strsplit('ENSG00000000003.13', '[.]')[[1]][1]
[1] "ENSG00000000003"
合并成功
引入循環(huán)，用stringr包
library(stringr)
a$ensembl_id=str_split(a$V1, '[.]', simplify = T)[,1]
我們就在a中加入了名為ensembl_id，合并數(shù)據(jù)的一列
image.png
載入包
library(org.Hs.eg.db)
ls("package:org.Hs.eg.db")
g2s=toTable(org.Hs.egSYMBOL);head(g2s)
g2e=toTable(org.Hs.egENSEMBL);head(g2e)
把a數(shù)據(jù)框與g2e關(guān)聯(lián)
b <- merge(a,g2e,by='ensembl_id', all.x=T)
d <- merge(b,g2s,by='gene_id', all.x=T)
最后在d中恢復與a相同的排序
d <- d[order(d$V1),]
image.png
去除d中的重復序號
d <- d[!duplicated(d$V1)，]
合并d與a的順序
d <- d[match(a$V1,d$V1),]
最后導出成csv文件
write.csv(d,'geneID2symbol.csv')
image.png
8 任意基因任意癌癥表達量分組的生存分析
從生存分析網(wǎng)頁工具下載數(shù)據(jù)
http://www.oncolnc.org/kaplan/?lower=50&upper=50&cancer=LGG&gene_id=93663&raw=ARHGAP18&species=mRNA
導入R
options(stringsAsFactors = F)
a <- read.table('R_bilibili-master/R語言小作業(yè)/LGG_93663_50_50.csv',
header = T,
sep = ',',
fill = T)
colnames(a)
dat <- a
繪制第一個圖
library(ggstatsplot)
ggbetweenstats(data = dat, x= Group, y= Expression)
image.png
繪制生存曲線
table(dat$Status)
dat$Status <- ifelse(dat$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=dat)
summary(sfit)
ggsurvplot(sfit, conf.int = F, pval = T)
image.png
更進一步：
ggsurvplot(sfit, palette = c("#E7B800", "#2E9FDF"),
risk.table = T, pval = T,
conf.int = T,
xlab="Time in months",
ggtheme = theme_light(),
ncensor.plot=T)
保存圖片
ggsave('survial_ARHGAP18_in_LGG.png')
image.png
9 任意基因任意癌癥表達量和臨床性關(guān)聯(lián)
下載plot文件。
options(stringsAsFactors = F)
a=read.table('R_bilibili-master/R語言小作業(yè)/plot.txt',
sep = '\t',fill = T,header = T)
a <- a[,-5]
colnames(a)=c('id','subtype','expression','mut')
dat=a
library(ggstatsplot)
ggbetweenstats(data = dat, x = subtype, y = expression)
image.png
中間出現(xiàn)了一段報錯，原因是Package 'PMCMRplus'沒有安裝，報錯顯示安裝即可畫出。
10 表達矩陣的樣本相關(guān)性
獲取airway數(shù)據(jù)：
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("airway")
library(airway)
options(stringsAsFactors = F)
library(airway)
data("airway")
exprSet <- assay(airway)
colnames(exprSet)
查看樣本：
> dim(exprSet)
[1] 64102 8
8列，6萬多個值。
第一列和第二列相關(guān)性
> cor(exprSet[,1],exprSet[,2])
[1] 0.9632268
group_list <- colData(airway)[,3]
tmp <- data.frame(g=group_list)
rownames(tmp) <- colnames(exprSet)
pheatmap::pheatmap(cor(exprSet),annotation_col = tmp)
image.png
刪除數(shù)據(jù)里的空白數(shù)值，判斷每一行是否都有數(shù)據(jù)。
exprSet <- exprSet[apply(exprSet,1, function(x) sum(x>1)>5),]
篩選后從原先
> dim(exprSet)
[1] 64102 8
6萬多個基因變成
> dim(exprSet)
[1] 19481 8
接近2萬個。
進一步整理數(shù)據(jù)，需要安裝edgeR包：
> BiocManager::install("edgeR")
取mad數(shù)值最大前500個基因畫圖
exprSet <- log(edgeR::cpm(exprSet)+1)
exprSet <- exprSet[names(sort(apply(exprSet,1,mad),decreasing = T)[1:500]),]
M <- cor(log2(exprSet+1))
tmp <- data.frame(g=group_list)
rownames(tmp) <- colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)
調(diào)整后的圖形：
image.png
11 芯片表達矩陣下游分析
獲取CLL包的數(shù)據(jù)：
if (!require("BiocManager", quietly = TRUE))
install.packages("BiocManager")
BiocManager::install("CLL")
正式開始：
library(CLL)
data(sCLLex)
sCLLex
exprSet <- exprs(sCLLex)
samples <- sampleNames(sCLLex)
pdata <- pData(sCLLex)
group_list <- as.character(pdata[,2])
dim(exprSet)
[1] 12625 22
exprSet[1:5,1:5]
CLL11.CEL CLL12.CEL CLL13.CEL CLL14.CEL CLL15.CEL
1000_at 5.743132 6.219412 5.523328 5.340477 5.229904
1001_at 2.285143 2.291229 2.287986 2.295313 2.662170
1002_f_at 3.309294 3.318466 3.354423 3.327130 3.365113
1003_s_at 1.085264 1.117288 1.084010 1.103217 1.074243
1004_at 7.544884 7.671801 7.474025 7.152482 6.902932
獲得差異矩陣
boxplot(exprSet)
構(gòu)造比較矩陣
#DEG by limma
suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design) <- levels(factor(group_list))
rownames(design) <- colnames(exprSet)
design
progres. stable
CLL11.CEL 1 0
CLL12.CEL 0 1
CLL13.CEL 1 0
CLL14.CEL 1 0
CLL15.CEL 1 0
略
把progress組與stable組進行差異分析比較
contrast.matrix <- makeContrasts(paste0(unique(group_list),collapse="-"),levels=design)
contrast.matrix
Contrasts
Levels progres.-stable
progres. 1
stable -1
差異分析第一步：
fit <- lmFit(exprSet,design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)
tempOutput <- topTable(fit2, coef=1, n=Inf)
nrDEG <- na.omit(tempOutput)
最后得到nrDEG，畫火山圖，富集分析等。
image.png
12 RNA-Seq表達矩陣差異分析
options(stringsAsFactors = F)
library(airway)
data(airway)
exprSet <- assay(airway)
colnames(exprSet)
group_list <- colData(airway)[,3]
exprSet <- exprSet[apply(exprSet,1, function(x) sum(x>1) >5),]
table(group_list)
方法：#DESeq2
if(T){
library(DESeq2)
(colData <- data.frame(row.names = colnames(exprSet),
group_list = group_list))
dss <- DESeqDataSetFromMatrix(countData = exprSet,
colData = colData,
design = ~group_list)
tmp_f <- 'airway_DESeq2-dss.Rdata'
if(!file.exists(tmp_f)){
dss <- DESeq(dss)
save(dss, file = tmp_f)
}
load(file=tmp_f)
res <- results(dds,
contrasts=c("group_list","trt","untrt"))
resOrdered <- res[order(res$padj),]
head(resOrdered)
DEG <- as.data.frame(resOrdered)
DESeq2_DEG <- na.omit(DEG)
nrDEG <- DESeq2_DEG[,c(2,6)]
colnames(nrDEG) <- c('log2FoldChange','pvalue')
}
colnames(nrDEG) <- c('logFC','P.Value')
attach(nrDEG)
plot(logFC,-log10(P.Value))
library(ggpubr)
df <- nrDEG
df$v <- -log10(P.Value)
ggscatter(df, x="logFC", y="v", size=0.5)
df$g=ifelse(df$P.Value>0.01, 'stable',
ifelse(df$logFC>1.5,'up',
ifelse(df$logFC < -1.5,'down','stable')))
table(df$g)
df$name <- rownames(df)
ggscatter(df, x="logFC", y='v', size=0.5, color='g')
ggscatter(df, x="logFC", y='v',color='g', size=0.5,
palette=c("#00AFBB","#E7B800","#FC4E07"))
13 R語言習題
作業(yè) 1
根據(jù)R包org.Hs.eg.db找到下面ensembl 基因ID 對應(yīng)的基因名(symbol)
ENSG00000000003.13
ENSG00000000005.5
ENSG00000000419.11
ENSG00000000457.12
ENSG00000000460.15
ENSG00000000938.11
提示：
library(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
rm(list = ls())
options(stringsAsFactors = F)
a <- read.table('R_bilibili-master/R語言小作業(yè)/e1.txt')
library(org.Hs.eg.db)
g2s=toTable(org.Hs.egSYMBOL)
g2e=toTable(org.Hs.egENSEMBL)
head(g2e)
library(stringr)
a$ensembl_id <- unlist(lapply(a$V1,function(x){
strsplit(x,'[.]')[[1]][1]
})
)
tmp = merge(a,g2e,by='ensembl_id')
tmp = merge(tmp,g2s,by='gene_id')
作業(yè) 2
根據(jù)R包hgu133a.db找到下面探針對應(yīng)的基因名(symbol)。
1053_at
117_at
121_at
1255_g_at
1316_at
1320_at
1405_i_at
1431_at
1438_at
1487_at
1494_f_at
1598_g_at
160020_at
1729_at
177_at
提示：
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
兩種方法，tmp1修改a列名，tmp2不修改分別設(shè)置
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('R_bilibili-master/R語言小作業(yè)/e2.txt')
colnames(a)='probe_id'
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
tmp1=merge(ids,a,by='probe_id')
tmp2=ids[match(a$probe_id,ids$probe_id),]
比較tmp1tmp2：tmp1==tmp2
作業(yè) 3
找到R包CLL內(nèi)置的數(shù)據(jù)集的表達矩陣里面的TP53基因的表達量，并且繪制在 progres.-stable分組的boxplot圖。
提示：
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex)
library(hgu95av2.db)
suppressPackageStartupMessages(library(CLL))
data(sCLLex)
sCLLex
exprSet=exprs(sCLLex)
library(hgu95av2.db)
pd=pData(sCLLex)
ids=toTable(hgu95av2SYMBOL)
head(ids)
boxplot(exprSet['1939_at',] ~ pd$Disease)
boxplot(exprSet['1974_s_at',] ~ pd$Disease)
boxplot(exprSet['31618_at',] ~ pd$Disease)
作業(yè) 4
找到BRCA1基因在TCGA數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集(Breast Invasive Carcinoma (TCGA, PanCancer Atlas))的表達情況。
提示：使用http://www.cbioportal.org/index.do 定位數(shù)據(jù)集：http://www.cbioportal.org/datasets
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('R_bilibili-master/R語言小作業(yè)/e4-plot.txt',
sep = '\t',
fill = T,
header = T)
colnames(a)=c('id','subtype','expression','mut')
dat=a
dat=dat[,-4]
library(ggstatsplot)
ggbetweenstats(data = dat, x = subtype, y = expression)
lastlibrary(ggplot2)
ggsave('plot-again-BRCA1-TCGA-BRCA-cbioportal.png')
tips:如果出現(xiàn)報錯，Names must be unique.，那么將expression列名改為exp即可成功繪圖。（可能是expression與某一函數(shù)同名會起沖突）
作業(yè) 5
找到TP53基因在TCGA數(shù)據(jù)庫的乳腺癌數(shù)據(jù)集的表達量分組看其是否影響生存
提示使用：http://www.oncolnc.org/
rm(list = ls())
options(stringsAsFactors = F)
a=read.table('BRCA_7157_50_50.csv',sep = ',',fill = T,header = T)
dat=a
library(ggplot2)
library(survival)
library(survminer)
table(dat$Status)
dat$Status=ifelse(dat$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=dat)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
ggsave('survival_TP53_in_BRCA_TCGA.png')
### 分割線
head(a)
b=read.table('e4-plot.txt',sep = '\t',fill = T,header = T)
colnames(b)=c('Patient','subtype','expression','mut')
head(b)
b$Patient=substring(b$Patient,1,12)
tmp=merge(a,b,by='Patient')
table(tmp$subtype)
type='BRCA_LumB'
x=tmp[tmp$subtype==type,]
library(ggplot2)
library(survival)
library(survminer)
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
table(tmp$subtype)
type='BRCA_Normal'
x=tmp[tmp$subtype==type,]
library(ggplot2)
library(survival)
library(survminer)
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
table(tmp$subtype)
type='BRCA_Basal'
x=tmp[tmp$subtype==type,]
library(ggplot2)
library(survival)
library(survminer)
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
table(tmp$subtype)
type='BRCA_Her2'
x=tmp[tmp$subtype==type,]
library(ggplot2)
library(survival)
library(survminer)
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
table(tmp$subtype)
type='BRCA_LumA'
x=tmp[tmp$subtype==type,]
library(ggplot2)
library(survival)
library(survminer)
#table(x$Status)
x$Status=ifelse(x$Status=='Dead',1,0)
sfit <- survfit(Surv(Days, Status)~Group, data=x)
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)
作業(yè)6
下載數(shù)據(jù)集GSE17215的表達矩陣并且提取下面的基因畫熱圖
ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T
提示：根據(jù)基因名拿到探針I(yè)D，縮小表達矩陣繪制熱圖，沒有檢查到的基因直接忽略即可。
rm(list = ls())
options(stringsAsFactors = F)
# 注意查看下載文件的大小，檢查數(shù)據(jù)
f='GSE17215_eSet.Rdata'
library(GEOquery)
# 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE17215', destdir=".",
AnnotGPL = F, ## 注釋文件
getGPL = F) ## 平臺文件
save(gset,file=f) ## 保存到本地
}
load('GSE17215_eSet.Rdata') ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因為這個GEO數(shù)據(jù)集只有一個GPL平臺，所以下載到的是一個含有一個元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
library(hgu133a.db)
ids=toTable(hgu133aSYMBOL)
head(ids)
dat=dat[ids$probe_id,]
dat[1:4,1:4]
ids$median=apply(dat,1,median)
ids=ids[order(ids$symbol,ids$median,decreasing = T),]
ids=ids[!duplicated(ids$symbol),]
dat=dat[ids$probe_id,]
rownames(dat)=ids$symbol
dat[1:4,1:4]
dim(dat)
ng='ACTR3B ANLN BAG1 BCL2 BIRC5 BLVRA CCNB1 CCNE1 CDC20 CDC6 CDCA1 CDH3 CENPF CEP55 CXXC5 EGFR ERBB2 ESR1 EXO1 FGFR4 FOXA1 FOXC1 GPR160 GRB7 KIF2C KNTC2 KRT14 KRT17 KRT5 MAPT MDM2 MELK MIA MKI67 MLPH MMP11 MYBL2 MYC NAT1 ORC6L PGR PHGDH PTTG1 RRM2 SFRP1 SLC39A6 TMEM45B TYMS UBE2C UBE2T'
ng=strsplit(ng,' ')[[1]]
table(ng %in% rownames(dat))
ng=ng[ng %in% rownames(dat)]
dat=dat[ng,]
dat=log2(dat)
pheatmap::pheatmap(dat,scale = 'row')
作業(yè)7
下載數(shù)據(jù)集GSE24673的表達矩陣計算樣本的相關(guān)性并且繪制熱圖，需要標記上樣本分組信息
rm(list = ls())
options(stringsAsFactors = F)
# 注意查看下載文件的大小，檢查數(shù)據(jù)
f='GSE24673_eSet.Rdata'
library(GEOquery)
# 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE24673', destdir=".",
AnnotGPL = F, ## 注釋文件
getGPL = F) ## 平臺文件
save(gset,file=f) ## 保存到本地
}
load('GSE24673_eSet.Rdata') ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因為這個GEO數(shù)據(jù)集只有一個GPL平臺，所以下載到的是一個含有一個元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
group_list=c('rbc','rbc','rbc',
'rbn','rbn','rbn',
'rbc','rbc','rbc',
'normal','normal')
dat[1:4,1:4]
M=cor(dat)
pheatmap::pheatmap(M)
tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
pheatmap::pheatmap(M,annotation_col = tmp)
作業(yè)8
找到 GPL6244 platform of Affymetrix Human Gene 1.0 ST Array 對應(yīng)的R的bioconductor注釋包，并且安裝它
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/"))
BiocManager::install("hugene10sttranscriptcluster.db",ask = F,update = F)
options()$repos
options()$BioC_mirror
作業(yè)9
下載數(shù)據(jù)集GSE42872的表達矩陣，并且分別挑選出 所有樣本的(平均表達量/sd/mad/)最大的探針，并且找到它們對應(yīng)的基因。
rm(list = ls())
options(stringsAsFactors = F)
# 注意查看下載文件的大小，檢查數(shù)據(jù)
f='GSE42872_eSet.Rdata'
library(GEOquery)
# 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE42872', destdir=".",
AnnotGPL = F, ## 注釋文件
getGPL = F) ## 平臺文件
save(gset,file=f) ## 保存到本地
}
load('GSE42872_eSet.Rdata') ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因為這個GEO數(shù)據(jù)集只有一個GPL平臺，所以下載到的是一個含有一個元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# (平均表達量/sd/mad/)最大的探針
boxplot(dat)
sort(apply(dat,1,mean),decreasing = T)[1]
sort(apply(dat,1,sd),decreasing = T)[1]
sort(apply(dat,1,mad),decreasing = T)[1]
作業(yè)10
下載數(shù)據(jù)集GSE42872的表達矩陣，并且根據(jù)分組使用limma做差異分析，得到差異結(jié)果矩陣
rm(list = ls())
options(stringsAsFactors = F)
# 注意查看下載文件的大小，檢查數(shù)據(jù)
f='GSE42872_eSet.Rdata'
library(GEOquery)
# 這個包需要注意兩個配置，一般來說自動化的配置是足夠的。
#Setting options('download.file.method.GEOquery'='auto')
#Setting options('GEOquery.inmemory.gpl'=FALSE)
if(!file.exists(f)){
gset <- getGEO('GSE42872', destdir=".",
AnnotGPL = F, ## 注釋文件
getGPL = F) ## 平臺文件
save(gset,file=f) ## 保存到本地
}
load('GSE42872_eSet.Rdata') ## 載入數(shù)據(jù)
class(gset)
length(gset)
class(gset[[1]])
# 因為這個GEO數(shù)據(jù)集只有一個GPL平臺，所以下載到的是一個含有一個元素的list
a=gset[[1]]
dat=exprs(a)
dim(dat)
pd=pData(a)
# (平均表達量/sd/mad/)最大的探針
boxplot(dat)
group_list=unlist(lapply(pd$title,function(x){
strsplit(x,' ')[[1]][4]
}))
exprSet=dat
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
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)
contrast.matrix<-makeContrasts("progres.-stable",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)
