《GEO数据挖掘课程》配套练习题粗浅的答案 / 开普饭

前面我布置了一系列学徒作业，其中有一个虽然是单一作业，但是工作量相当于一个本科毕业设计：《GEO数据挖掘课程》配套练习题，关于这个课程学徒也写了一系列笔记：学徒写的《GEO数据挖掘课程》的配套笔记完结撒花

下面是GEO配套练习题的一个答案

学徒写的代码不一定说有多好，但是胜在努力学习了，同样的值得分享！（代码很多冗余的地方，粗制滥造，甚至会有一些些错误，但确实很真诚！）

基于转录本上的概念

==外显子和内含子==：基因 DNA 分为编码区和非编码区，编码区包含外显子和内含子，一般非编码区具有基因表达的调控功能，如启动子在非编码区。编码区则转录为 mRNA 并最终翻译成蛋白质。外显子和内含子都被转录到 mRNA 前体 hnRNA 中，当 hnRNA 进行剪接变为成熟的 mRNA 时，内含子被切除，而外显子保留。实际上真正编码蛋白质的是外显子，而内含子则无编码功能。内含子存在于DNA 中，在转录的过程中，DNA 上的内含子也会被转录到前体 RNA 中，但前体 RNA 上的内含子会在 RNA 离开细胞核进行翻译前被切除。promoter不属于intron和exon的任何一个，属于noncoding sequence. ==开放读码框====ORF== 开放读码框是从一个起始密码子开始到一个终止密码子结束的一段序列；不是所有读码框都能被表达出蛋白产物，或者能表达出占有优势或者能产生生物学功能的蛋白。内含子和外显子指的就是一个开放阅读框（ORF）内编码的部分和不编码的部分。

基于翻译上的概念

mRNA 包括 UTR 和 CDS ==UTR==：UTR,Untranslated Regions一般指的是一个转录本（transcript）3'和5'不参与编码的区域即非翻译区，是信使RNA（mRNA）分子两端的非编码片段。UTR区不参与编码，但是不是说他们没有功能，只是不被翻译成具有功能的蛋白质。多数基因都有UTR,它们也是外显子拼接的产物。UTR在DNA序列中是算外显子的区域。==CDS==：CDS,Sequencecodingfor aminoacids in protein 蛋白质编码区，CDS 是 Coding sequence的缩写，是编码一段蛋白产物的序列，CDS 必定是一个 ORF 。但也可能包括很多 ORF 。反之，每个 ORF 不一定都是 CDS 。外显子与 CDS 区不是完全一致的，CDS 区一定属于外显子，但是外显子不一定是 CDS 区，也就是说外显子不一定都能翻译成蛋白。

外显子和内含子是转录本上面的概念，是基于转录这个行为的定义。 而5端UTR区域和CDS区域，还有3端UTR区域,是基于翻译这种行为的定义！

背景介绍完毕，正式开始

哦下面正式开始分类，其实还是得自己去了解背景，呃，以前没有这么深入细致的区别这些知识点，真的有点难懂。（总以为自己在学习生物信息学，其实何尝不是在学生物呢？）

主要的目的是要从包里挑出lncRNA 基因的探针....

library(hgu133plus2.db) ids=toTable(hgu133plus2SYMBOL) length(unique(ids$symbol)) tail(sort(table(ids$symbol))) table(sort(table(ids$symbol))) plot(table(sort(table(ids$symbol)))) # 可以看到很多基因都是有着多个对应的探针

从ensembl网站下载了最新的注释文件 Homo_sapiens.GRCh38.101.chr.gtf.gz（就这轻飘飘一句话，对初学者来说，可以跑一天一夜，都不一定能下载到gtf文件。）使用R包读入：

library(rtracklayer) library(dplyr) gtf <- import('Homo_sapiens.GRCh38.101.chr.gtf.gz') gtf <- gtf[gtf$gene_name!=''] gtf <- gtf[!is.na(gtf$gene_name)] save(gtf,file = 'gtf.Rdata') table(gtf$transcript_biotype)

既然已经做到了这一步，顺便也把别的相对应的基因也分类了

首先是蛋白，竟然这么多....为什么不分在一起

###可以参考下面的网页来进行分类 #http://vega.archive.ensembl.org/info/about/gene_and_transcript_types.html #http://asia.ensembl.org/Help/Glossary #https://www.gencodegenes.org/pages/biotypes.html


##protein coding

protein = gtf$gene_name[gtf$transcript_biotype %in% 

 c('nonsense_mediated_decay','nontranslating_CDS','non_stop_decay',

 'polymorphic_pseudogene','protein_coding')]

gtf_df <- as.data.frame(gtf) table(gtf_df$transcript_biotype) t_df <- select(gtf_df, c(gene_id,gene_name,transcript_biotype)) dim(t_df)


## Long non-coding RNA (lncRNA)

lncRNA <- t_df$gene_name[t_df$transcript_biotype %in% 

  c('processed_transcript','lncRNA','Non coding',

 '3prime_overlapping_ncRNA','Antisense','lincRNA','Retained_intron',

 'Sense_intronic','Sense_overlapping','macro_lncRNA')]

最后注释到一起

lncRNA <- as.data.frame(lncRNA) a = lncRNA a = a[!duplicated(a$lncRNA),] a = as.data.frame(a)

ids1 = ids[match(a$a,ids$symbol),] ids2 = na.omit(ids1) dim(ids2)

我没看出来哪里出了问题.....应该得到1950个探针...而我的结果是11488，多了这么多....

什么原因？？？

看了下只包括lncRNA的

lncRNA <- t_df$gene_name[t_df$transcript_biotype %in% c('lncRNA')]

确实最后的结果是只有1780，和需要的1950个探针结果相似，也就是说，并不需要那么多条目，只需要lncRNA？

问题2：GEO数据集下载

下载GEO dataset (GSE4290) 数据集，使用GEOquery包或者其它。然后筛选 77 个 glioblastoma samples and 23 non-tumor controls

首先下载数据

rm(list = ls()) options(stringsAsFactors = F) library(AnnoProbe) library(GEOquery) gse=geoChina('GSE4290') eSet=gse[[1]]

probes_expr <- exprs(eSet);dim(probes_expr) head(probes_expr[,1:4]) boxplot(probes_expr,las=2) probes_expr=log2(probes_expr+1) phenoDat <- pData(eSet) head(phenoDat[,1:4])

可以看到临床信息中

metadata = phenoDat head(metadata)


index1=grep('glioblastoma',metadata$characteristics_ch1)

metadata1=metadata[index1,]
index2=grep('non-tumor',metadata$characteristics_ch1)

metadata2=metadata[index2,]
###这就挑出来了 77 个 glioblastoma samples and 23 non-tumor controls

##临床信息

md = rbind(metadata1,metadata2)

##匹配表达矩阵 probes_expr = probes_expr[,colnames(probes_expr) %in% md$geo_accession] dim(probes_expr) #[1] 54613 100

这个就蛮简单的了，只需要基础的R语言过关即可。（但这些代码写的超级冗余，超级小白，我相信生信技能树的绝大部分读者，都比我写的好！！！）

问题3：差异分析并且绘制热图

在下载的GSE4290表达矩阵里面提取1303 lncRNA基因的热图。（数量级是1000即可，不要求精确数量）

library(limma) exp = normalizeBetweenArrays(probes_expr) boxplot(exp,las = 2)


#group_list----

md$characteristics_ch1

group_list=c(rep("glioblastoma",77),rep("non-tumor",23))

group_list
#(4)提取芯片平台编号

gpl <- gset[[1]]@annotation

gpl 

#"GPL570"

#http://www.bio-info-trainee.com/1399.html

library(AnnoProbe)

##经常无缘无故的抽风~有时能下有时不能下，古怪

tmp <- idmap(gpl,type = 'pipe')
###如果上面不能下载，替代方案如下

if(F){

  if(!require(hgu133plus2.db))BiocManager::install("hgu133plus2.db")

  library(hgu133plus2.db)

  ls("package:hgu133plus2.db")

  ids <- toTable(hgu133plus2SYMBOL)

  head(ids)

}else if(F){

  getGEO(gpl)

  ids = data.table::fread(paste0(gpl,".soft"),header = T,skip = "ID",data.table = F)

  ids = ids[c("ID","Gene Symbol"),]

  colnames(ids) = c("probe_id")

}else if(T){

  ids = idmap(gpl,type = "bioc")

}
tmp = ids

##删除了10个，说明注释中有10个不在表达矩阵中

tmp = tmp[tmp$probe_id %in% rownames(exp),]
tmp2 = annoGene(IDs = tmp$symbol, ID_type = 'SYMBOL', species = "human")

table(tmp2$biotypes)

sort(table(tmp2$biotypes))

##tail就直接看到了最后面最大的往前数

tail(sort(table(tmp2$biotypes)))

##而head排序后从前面往后数，由于没有那么复杂的加上T

head(sort(table(tmp2$biotypes)))

这样就可以得到 tmp2，啊这是新学会的方法，在Jimmy老师的B站的lncRNA视频学的~

可以看到 biotypes 这一列中有我们所需要的lncRNA信息。我做出来直接有1601个而题目写的是1303个，哪里出了什么问题吗，但是这标的清清楚楚明明白白：

不过题目也强调了，只需要数据量是ok的即可，都是一千多，不需要精确的一模一样。

######################挑选lncRNA head(tmp2[tmp2$biotypes == 'lncRNA',]) ##首先挑取tmp2的biotypes这一列为编码的所有行，取第一列symbol ##于是pcg就是19411个symbol,且去除重复探针 lncg = unique(tmp2[tmp2$biotypes == 'lncRNA',1]) ##挑选tmp中对应的pcg中的基因名，就可以得到对应的第一列probe_id lncp = tmp[tmp$symbol %in% lncg,1] ###重构表达矩阵，现在是想要的探针组成 lnc_exp = exp[lncp,]


############################## 匹配注释信息

dat = lnc_exp

ids = tmp

dat[1:5,1:5]  

head(ids)

ids=ids[ids$probe_id %in%  rownames(dat),]

##%in%用于判断是否匹配，然后取匹配的几行，去掉无法匹配的信息。

##去除基因的冗余，一个探针是可以对应多个基因

ids=ids[!duplicated(ids$probe_id),]

#去除探针的冗余：取表达矩阵中可以与探针名匹配的那些，去掉无法匹配的表达数据

dat=dat[ids$probe_id,]
ids$median=apply(dat,1,median) 

#ids新建median这一列，列名为median，同时对dat这个矩阵按行操作，取每一行的中位数，将结果给到median这一列的每一行

ids=ids[order(ids$symbol,ids$median,decreasing = T),]

#对ids$symbol按照ids$median中位数从大到小排列的顺序排序

##即先按symbol排序，相同的symbol再按照中位数从大到小排列，方便后续保留第一个值。

##将对应的行赋值为一个新的ids，这样order()就相当于sort()

ids=ids[!duplicated(ids$symbol),]

#将symbol这一列取取出重复项，'!'为否，即取出不重复的项，去除重复的gene ，保留每个基因最大表达量结果

dat=dat[ids$probe_id,] 

#新的ids取出probe_id这一列，将dat按照取出的这一列中的每一行组成一个新的dat

rownames(dat)=ids$symbol

#把ids的symbol这一列中的每一行给dat作为dat的行名

head(dat)

##至此我们就得到了该数据集的表达矩阵，最后将结果保存。
save(dat,md,group_list,file = 'lnc.Rdata')

接下来就可以直接画热图了，中间还有一个PCA图，我看着觉得不怎么能区分的开，也没见着有批次效应！

#热图 rm(list = ls()) load(file = "lnc.Rdata") library(pheatmap) cg=names(tail(sort(apply(dat,1,sd)),1000)) n=dat[cg,] pheatmap(dat[cg,],show_rownames = F,show_colnames = F) ###scale = "row"，为了重点突出样本与样本之间的差异而缩小基因与基因之间的差异，做了标准化 n2 = t(scale(t(n))) n2[n2 > 2] =2 n2[n2 < -2] = -2


pheatmap(n2,show_rownames = F,show_colnames = F)
#绘制热图

###annotation_col是为了加上样本分组

annotation_col=data.frame(group=group_list)

rownames(annotation_col)=colnames(n2)
pheatmap(n2,

   show_colnames =F,

   show_rownames = F,

   cluster_cols = F,

   annotation_col=annotation_col)

#scale = "row"

pheatmap(n2,

   show_colnames =F,

   show_rownames = F,

   #cluster_cols = F,

   annotation_col=annotation_col,

   filename = 'heatmap')

dev.off()

这个热图太丑了...看不下去

肯定是中间出了什么问题？？？（头发都掉完了，唉！！！！不想学生信了，不想这么辛苦了）

问题4：绘制火山图

仅仅是对1303个 lncRNA基因的表达矩阵进行差异分析，并且根据 (|logFC| >1; P <0.05), 阈值来画火山图

#差异分析——limma rm(list = ls()) ## 魔幻操作，一键清空~ options(stringsAsFactors = F) load(file = "lnc.Rdata") exp = dat dim(exp)


library(limma)

# 做分组矩阵 

group_list=c(rep("glioblastoma",77),rep("nontumor",23))

design <- model.matrix(~0+factor(group_list))

colnames(design)=levels(factor(group_list))

rownames(design)=colnames(exp)

design  #分组矩阵
# 做比较矩阵
# contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = "-"),levels = design)

# contrast.matrix ##这个矩阵声明，我们要把treat组和contorl组进行差异分析比较

# -1和1的意思是contorl是用来被比的，treat是来比的

contrast.matrix<-makeContrasts(paste0(c("glioblastoma","nontumor"),collapse = "-"),levels = design)

contrast.matrix

#到此，做差异分析所需要的三个矩阵就做好了：表达矩阵、分组矩阵、差异比较矩阵

#我们已经制作好了必要的输入数据，下面开始讲如何使用limma这个包来进行差异分析了！
##step1

fit <- lmFit(exp,design)

##step2

fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果

fit2 <- eBayes(fit2)  ## default no trend !!!

##eBayes() with trend=TRUE

##step3

allDiff <- topTable(fit2, adjust="fdr", sort.by="B", number=50000)
nrDEG = na.omit(allDiff) 

#write.csv(nrDEG2,"limma_notrend.results.csv",quote = F)

head(nrDEG)
#####火山图1ok

library(ggplot2)

library(ggpubr)

library(base)
##  新加一列Group

nrDEG$Group = "not-significant"

##  将adj.P.Val小于0.05，logFC大于2的基因设为显著上调基因

##  将adj.P.Val小于0.05，logFC小于2的基因设为显著下调基因

nrDEG$Group[which((nrDEG$P.Value < 0.05)&(nrDEG$logFC > 1))] = "up"

nrDEG$Group[which((nrDEG$P.Value < 0.05)&(nrDEG$logFC < -1))] = "down"

##  查看上调和下调基因数目

table(nrDEG$Group)
##3.可以挑出上调和下调差异基因

deg = nrDEG

colnames(deg)

deg$symbol = rownames(deg)
x1 = head(deg[deg$Group=='up',],2)

x2 = head(deg[deg$Group=='down',],2)
for_label = rbind(x1,x2)

logFC_t = 1 P.Value_t = 0.05 p <- ggplot(data = deg, aes(x = logFC, y = -log10(P.Value))) + geom_point(alpha=0.4, size=3.5, aes(color=Group)) + ylab("-log10(Pvalue)")+ scale_color_manual(values=c("blue", "grey","red"))+ geom_vline(xintercept=c(-logFC_t,logFC_t),lty=4,col="black",lwd=0.8) + geom_hline(yintercept = -log10(P.Value_t),lty=4,col="black",lwd=0.8) + theme_bw() p volcano_plot <- p + geom_point(size = 3, shape = 1, data = for_label) + ggrepel::geom_label_repel( aes(label = symbol), data = for_label, color="black" ) volcano_plot

总共就两个上调2个下调...（我这是做错了什么，我找谁哭去呢？？？）

感觉实在是做不下去了，唉，都是坑啊！！！

问题5：R包和GPL的soft信息差异

比较hgu133plus2.db里面的基因的注释信息和https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GPL570的数据的区别。

主要是更新时间问题，看：hgu133plus2_dbInfo()

1558290_a_at Homo sapiens BG200951 PVT1

该平台对应的bioconductor里面的芯片探针注释包hgu133plus2的信息的提取，而soft就是我们前面说的在GEO数据库里面访问该平台的主页看到的注释信息的提取

###1.通过Jimmy老师的这个包获取注释平台 library(idmap2) id=get_soft_IDs('GPL570')

###2.通过R包注释 library(hgu133plus2.db) ##运行这个包 ls("package:hgu133plus2.db") #对这个包进行探索，看一下有多少元素 ##找到有SYMBOL的那个元素就是我们需要的对应关系 ids=toTable(hgu133plus2SYMBOL)

###3.直接下载平台注释信息 library(GEOquery) #Download GPL file, put it in the current directory, and load it: gpl <- getGEO('GPL570', destdir=".") ##需要修改相应的平台号，把平台信息赋值给gpl colnames(Table(gpl)) head(Table(gpl)[,c(1,11)]) ## you need to check this , which column do you need probe2gene=Table(gpl)[,c(1,11)] head(probe2gene)

通过三种方式获取平台的注释信息，理论上将第一种方式和第三种方式是完全一样的。

主要比较2.3的区别

可以看到，查同一个基因，会出来不一样的探针结果。

###GPL570 hgu133plus2_dbInfo()

library(AnnoProbe) ids1=idmap('GPL570',type = 'bioc') head(ids1) ids2=idmap('GPL570',type = 'soft') head(ids2)


length(unique(ids1[,1]))

length(unique(ids2[,1]))

ids2_filter=ids2[nchar(ids2[,2])>1,] m1=merge(ids1,ids2_filter,by.x = 'probe_id',by.y = 'ID') table(m1[,2]==m1[,3]) #FALSE TRUE # 2495 38953 m1_not=m1[m1[,2] !=m1[,3],]

可以看到，两个注释结果冲突的比例还好，可以接受的误差范围。毕竟38953 探针都是一致的，然后检查那2495个冲突了的探针，在 m1_not 中随便找一个基因看一下，比如说这个

谷歌搜索之后，出来了下面这个...说明了什么，应该是说明了下面这个也就是第二列的注释差不多还行

gencode数据库的注释，发现来源于soft的注释信息（第三列）的基因名，99.48%都是无法去找到记录的，手动搜索一下，发现都是一些被废弃掉的基因名字。

tmp1=annoGene(m1_not[,2],'SYMBOL') tmp2=annoGene(m1_not[,3],'SYMBOL')

所以选择R包就好，平台信息可能有点老，不怎么更新。

问题6: 指定基因（这里是PVT1）画boxplot，多个分组

还是前面的 GSE4290 数据，突然发现这样把每个分组单独提出来，很方便分组耶....如果不希望在分组上纠结一下。其实也差不多..:candy:

rm(list = ls()) options(stringsAsFactors = F) library(AnnoProbe) library(GEOquery) load('GSE4290_eSet.Rdata') #gse=geoChina('GSE4290') eSet=gset[[1]]


probes_expr <- exprs(eSet);dim(probes_expr)

head(probes_expr[,1:4])

#boxplot(probes_expr,las=2)

probes_expr=log2(probes_expr+1)

phenoDat <- pData(eSet)

head(phenoDat[,1:4])
metadata = phenoDat

head(metadata)
index1=grep('glioblastoma',metadata$characteristics_ch1)

metadata1=metadata[index1,]
index2=grep('non-tumor',metadata$characteristics_ch1)

metadata2=metadata[index2,]
index3=grep('astrocytoma',metadata$characteristics_ch1)

metadata3=metadata[index3,]
index4=grep('oligodendroglioma',metadata$characteristics_ch1)

metadata4=metadata[index4,]
length(nrow(metadata1))

md = rbind(metadata1,metadata2,metadata3,metadata4)
##匹配表达矩阵

probes_expr = probes_expr[,colnames(probes_expr) %in% md$geo_accession]

dim(probes_expr)

#[1] 54613   176

group_list=c(rep('nt',nrow(metadata2)), rep('GBM',nrow(metadata1)), rep('astrocytoma',nrow(metadata3)), rep('od',nrow(metadata4))) table(group_list)

给定基因的boxplot，这里是 PVT1

library(ggstatsplot) dat=as.data.frame(t(probes_expr)) dat$group=group_list # 1558290_a_at Homo sapiens BG200951 PVT1 ###这里由于没有进行注释，dat中只有探针信息，直接用探针来画图 ggbetweenstats(data = dat, x = group, y = '1558290_a_at')

函数画多个基因

a =probes_expr boxplot(a[5,]~group_list,annotation_col = group_list) g = function(i){ boxplot(a[rownames(a)[i],]~group_list,annotation_col = group_list) } g(20)

这里缺一张图！

问题7：摸索TCGA的GBM的临床信息，找到classical, neural, proneural and mesenchymal分类方式

建议使用ucsc的xena浏览器探索，下载数据。如果你网速太差，也可以看我这边备份的TCGA数据，就是来源于xena，ucsc的，都在，https://share.weiyun.com/5zLnKmO

需求最大的是tcga数据库的生存分析和表达量差异，看看这两个视频：

https://www.bilibili.com/video/av25643438?p=9
https://www.bilibili.com/video/av49363776?p=6

首先去xena下载表达矩阵和临床信息

表达矩阵

rm(list = ls()) options(stringsAsFactors = F)


## 建议大家先处理表达矩阵，处理过程中会发现表达矩阵中的样本数量少于样本信息中的样本数量，

## 所以，以表达矩阵中的样本数为准

if(!file.exists('TCGA-GBM.HiSeqV2.Rdata')){

  expr<-read.table("HiSeqV2.gz",header = T,row.names = 1,check.names = F)

  expr = as.matrix(expr)

  expr = expr[apply(expr, 1, function(x) sum(x > 1) > 9), ]

  dim(expr)

  expr[1:4,1:4]

  save(expr, file = 'TCGA-GBM.HiSeqV2.Rdata')

}else{

  load('TCGA-GBM.HiSeqV2.Rdata')

}

临床信息

phe = read.table('GBM_clinicalMatrix',header = T,sep = '\t',quote = '') ## 只保留有表达矩阵的样本信息 phe <- phe[phe[,1] %in% colnames(expr),] colnames(phe) phe = phe[phe[,28] != '',] phe = phe[,c(1,28)]

接下来就是分组

group_list = as.character( phe[,2] ) table(group_list) rownames(phe) = phe$sampleID ###看一下表达矩阵的列名是否和临床信息的行名匹配，若不匹配，则改一下 p = identical(rownames(phe),colnames(expr));p exp = expr[,match(rownames(phe),colnames(expr))] p = identical(rownames(phe),colnames(exp));p

group_list = as.character( phe[,2] ) table(group_list)

最后看一下感兴趣的基因在四种分型中的表达。

boxplot(exp['PTEN',]~group_list)

exp['PDGFRA',] exp['EGFR',] exp['PTEN',] exp['CDKN2A',] exp['IDH1',] exp['NF1',]

问题8：摸索TCGA的GBM的mRNA-seq表达矩阵（完成PVT1的boxplot）

下载mRNA-seq表达矩阵，counts矩阵，然后提取PVT1基因的表达量信息，根据临床信息，分类汇总后绘图。

1.下载表达矩阵数据

######step1-getdata rm(list = ls()) expr<-read.csv('TCGA-GBM.htseq_counts.tsv.gz',sep = '\t') colnames(expr) library(stringr) ##a中新增一列为ensembl_id赋值为矩阵第一列的值，经小数点切割后 expr$Ensembl_ID = str_split(expr$Ensembl_ID,'[.]',simplify = T)[,1] rownames(expr)<-expr$Ensembl_ID


expr<-expr[,-1]

colnames(expr)<-gsub('.','-',colnames(expr),fixed = T) 

head(expr)

phe<-read.csv('TCGA-GBM.GDC_phenotype.tsv.gz',sep = '\t') phe[1:4,1:4] colnames(phe)

2.提取lncRNA数据

（刚开始直接提取了mRNA数据，结果发现 PVT1 是 lncRNA ）

##2.lncRNA数据 library(rtracklayer) library(dplyr) gtf <- import('Homo_sapiens.GRCh38.101.chr.gtf.gz') gtf_df <- as.data.frame(gtf) gene_df <- select(gtf_df, c(gene_id,gene_name,gene_biotype)) index <- duplicated(gene_df$gene_id) gene_df = gene_df[!index,] dim(gene_df) lncRNA = gene_df[gene_df$gene_biotype == 'lncRNA',] dim(lncRNA) save(lncRNA,gene_df,file = 'lncRNA_df.Rdata') exprSet = expr[match(lncRNA$gene_id,rownames(expr)),] dim(exprSet) exprSet = na.omit(exprSet) dim(exprSet) save(exprSet,phe,file = 'data.Rdata')

3.ID转换

##3.ID转换，ensemblID转换为symbolID rm(list = ls()) options(stringsAsFactors = F) library(limma) load("data.Rdata") load("lncRNA_df.Rdata") exprSet$names = rownames(exprSet) exprSet$names = lncRNA[match(exprSet$names,lncRNA$gene_id),2] dim(exprSet) table(duplicated(exprSet$names)) # 有9个重复的基因名 # 对重复基因名取平均表达量，然后将基因名作为行名 exprSet = avereps(exprSet[,-ncol(exprSet)],ID = exprSet$names) dim(exprSet) save(exprSet, file = 'exprSet_names_by_symbol.Rdata')

4.数据整理

##4.数据整理 rm(list = ls()) options(stringsAsFactors = F) load('data.Rdata') load("exprSet_names_by_symbol.Rdata")


# 4.1 去除低表达量的基因

exprSet = as.data.frame(exprSet)

pick_row <- apply(exprSet, 1, function(x){

sum(x == 0) < 40

  })

exprSet1 <- exprSet[pick_row,]

exprSet1['PVT1',]
# 4.2 分组(癌症组织和癌旁组织)

library(stringr)

tumor <- colnames(exprSet1)[as.integer(substr(colnames(exprSet1),14,15)) < 10]

normal <- colnames(exprSet1)[as.integer(substr(colnames(exprSet1),14,15)) >= 10]

  

tumor_sample <- exprSet1[,tumor]

normal_sample <- exprSet1[,normal]

  

exprSet_by_group <- cbind(tumor_sample,normal_sample)

group_list <- c(rep('tumor',ncol(tumor_sample)),rep('normal',ncol(normal_sample)))

  

save(exprSet_by_group, group_list, file = 'exprSet_by_group_list.Rdata')

colnames(exprSet_by_group)

5.boxplot

##4.数据整理 rm(list = ls()) options(stringsAsFactors = F) load('data.Rdata') load("exprSet_by_group_list.Rdata")


pData = phe[match(colnames(exprSet_by_group),phe$submitter_id.samples),]

dim(pData)

colnames(pData)

phe_M<-pData[,c(1,5,33,39)]

table(phe_M$initial_pathologic_diagnosis_method)
dat = exprSet_by_group

class(dat)

###哦就是这里，我一直画不出来boxplot图，结果发现应该是matrix

dat = as.matrix(dat)

class(dat)

head(dat)

rownames(dat)

g = as.character(phe_M$initial_pathologic_diagnosis_method)

table(g)

boxplot(dat['PVT1',] ~ g,annotation_col = g)

##应该是个list吗 boxplot(dat) boxplot(a) boxplot(dat['PVT1',]~group_list)

表达矩阵改为matrix就行

根据临床信息中分组画图

根据临床信息中 pathologic_diagnosis_method

完美！（信心又回来了一点，感谢各位走过路过的哥哥姐姐点赞，我会继续努力的！！！）

文末友情推荐

。如果大家没有时间自行慢慢摸索着学习，可以考虑我们生信技能树官方举办的学习班。

《GEO数据挖掘课程》配套练习题粗浅的答案