学习者：骆栢维

题目来源：生信基石之R语言

中级10 个题目：http://www.bio-info-trainee.com/3750.html

备注:本文为笔者学习健明老师GitHub答案代码的学习收获！

作业1：根据R包org.Hs.eg.db找到下面ensembl 基因ID 对应的基因名(symbol)

library(org.Hs.eg.db)symbol1=toTable(org.Hs.egSYMBOL)ensembl1=toTable(org.Hs.egENSEMBL)head(ensembl1)#查看前6行head(symbol1)#查看前6行exprset1 'gene_id')#orexprset2 'gene_id')c 'ENSG00000000003.13',       'ENSG00000000005.5',       'ENSG00000000419.11',       'ENSG00000000460.15',       'ENSG00000000938.11')c substr(c,c#看一下吧answer #这里按照逻辑值提取元素##来探讨一下吧exprset1$ensembl_id %in%c#%in%表示前者是否在后者里面，返回逻辑值table(exprset1$ensembl_id %in%c)#ORc answer1 "ensembl_id",by.y=View(answer)

学习心得：数据的提取

在上一篇分享文章的Q6，笔者提及过：提取数据的方法可以根据“名称”、“逻辑值” 和 “位置”。在上面的题目中，笔者是通过逻辑值提取数据，在下面的“读入数据” 中笔者也是采用逻辑值的方式提取数据。

作业2：根据R包hgu133a.db找到下面探针对应的基因名(symbol)

install.packages("BiocManager")BiocManager::install('hgu133a.db')library(hgu133a.db)symbol2=toTable(hgu133aSYMBOL)#获取探针表格symbol2[1:4,]#看一下前4行c '1053_at',       '1320_at', '1405_i_at', '1431_at', '1438_at' ,'1487_at' ,'1494_f_at',       '1598_g_at', '160020_at', '1729_at','177_at' )answer2 in% c,]table(symbol2$probe_id %in% c)#看一下结果#另一种方法同Q1

作业3：找到R包CLL内置的数据集的表达矩阵里面的TP53基因的表达量，并且绘制在 progres.-stable分组的boxplot图

BiocManager::install('CLL')BiocManager::install('hgu95av2.db')library(CLL)data(package="CLL")#查看R包里面具体数据集data("sCLLex")#加载R包内置数据,是列表biaoxing #提取表型数据exprSet #提取表达矩阵 library(hgu95av2.db)symbol3=toTable(hgu95av2SYMBOL)# 在symbol3中搜索TP53得到三个探针"1939_at"；'31618_at'；1974_s_at'probe in%probeexprSet boxplot(exprSet["1939_at",] ~ biaoxing$Disease,data = exprSet,col="red") #其他的同上，换基因名就可以了

作业3(补充)：ggpubr

#ggpubr(以"1939_at"为例)library(ggpubr)#首先在绘图之前，我们要处理数据，把分组信息和表达信息合并起来#都是笔者发现他们的sampleID和表达数据的行名并不一致，怎么办呢？visual1 as.data.frame(exprSet)[biaoxing$SampleID2 function(x){  x  paste(c(x,"CEL"),collapse = ".")#paste函数主要是为了将两个字符串合并起来}))##或者#############################################biaoxing$SampleID1 function(x){  x  paste(c(x,"CEL"),collapse = ".")})##或者#######################################################上面是改分组信息名，也可以改表达数据的列名exprSet as.data.frame(exprSet)colnames(exprSet) ".CEL",colnames(exprSet)############################################################identical(colnames(exprSet),biaoxing$SampleID)#判断表达数据列名和分组信息是否对应，结果为Truebiaoxing$Disease as.character(biaoxing$Disease)b as.data.frame(t(biaoxing$Disease as.character(biaoxing$Disease)visual1 #按行合并visual1 %  as.matrix()%>%  t()%>%  as.data.frame()visual1$`1939_at` as.numeric(ggboxplot(visual1, "Disease", "`1939_at`",          color = "Disease", palette =c("#00AFBB", "#E7B800"),          add = "jitter", shape = "Disease")

实战一下吧——lncRNA数据处理(模拟)

##########################实战一下吧############################a 1:b 'tumor', c "male",e ":")d "_")shizhan #创造一个数据框#1  191:normal_female#2  131:tumor_male#3  13:normal_male#4  266:normal_female#5  257:normal_female#6  246:normal_female###如果要把他们分三列怎么做呢library(tidyr)shizhan$d ":",shizhan "a",#完美############################################另一种方法shizhan1 "a",###发现宝藏！！！！！#      a status gender#1   191 normal female#2   131  tumor   male#3    13 normal   male#4   266 normal female#5   257 normal female#6   246 normal female

学习心得

1.字符串的处理

*截取：在上一篇分享文章中，笔者介绍了unlist与strsplist截取字符串，还有笔者也拓展了substr和substring的用法

*合并：使用paste函数，在这里要注意有时候collapse参数得不到理想的结果，可以把参数变为sep，二者都可以作为连接参数(具体为什么collapse有时候得不到理想的结果，笔者也不确定)

*替换：使用gsub函数，可以做批量处理。

*拆分：使用separate函数，这里要注意的是，它第一个参数为数据框，第二个参数为列名，不要觉得“数据框$列名” 就可以了。笔者就犯了这么愚蠢的错误。

2.apply ; lapply 和sapply的区别

*输入的数据：apply为数据框和矩阵；lapply为列表，向量和数据框；sapply同lapply。

*输出结果：apply为向量，列表和数组；lapply为列表；sapply为向量和矩阵。

*回到上面的合并字符串的前两种方法，笔者在lapply前面加了unlist，而sapply却没有。可以理解吗？

温馨提示：由于笔者的的网络不好，没有办法下载所有数据(郁闷！！！)，所以笔者按照作业相应的考点进行操作。接下来的内容可能与原作业存在差异，还请大家海涵。

读入数据(对应作业8)

exprSet "GSE24673_series_matrix.txt",comment.char = library(tibble)rownames(exprSet) 1]exprSet -1]platformMap "platformMap.txt")platformMap[platformMap$gpl=="GPL6244",3]#找到对应的注释包"hugene10sttranscriptcluster"if (!requireNamespace("BiocManager", quietly = TRUE))  install.packages("BiocManager")BiocManager::install("hugene10sttranscriptcluster.db")library(hugene10sttranscriptcluster.db)prode head(prode)

基因注释(对应作业4和9)

exprSet$probe_id #将行名转变为第一列，笔者真是手残fexprset "probe_id")fexprset "probe_id",fexprset -1]fexprset 2:head(fexprset)#发现在第一列的是"RN7SK"基因maxgene "RN7SK",maxgenefexprset #去重rownames(fexprset) 1]final -1]save(final,file = "final_exprset.Rdata")#保存数据

备注：很感激果子老师为我们这些生信小白这里的平台文件的注释R包！详情请看：https://mp.weixin.qq.com/s/nWbMO4mULgN__nPjooRDlg

差异分析(对应作业5、7和10一部分)

library(limma)group "tumor",group group,levels = c(design group)colnames(design) group)design#  tumor normal#1      1      0#2      1      0#3      1      0#4      1      0#5      1      0#6      1      0#7      1      0#8      1      0#9      1      0#10     1      1#11     1      1fit #线性模型拟合fit2 #贝叶斯检验allDiff=topTable(fit2,adjust='fdr',coef=2,number=Inf) #输出结果save(allDiff,file = "allDiff.Rda")#保存dif_subset 1 & adj.P.Val < 0.05)#找出差异两倍以上，pvalue小于0.05#其他方法##################################################3library(tibble)dif_filter %   rownames_to_column() %>%   filter(adj.P.Val < 0.05 & abs(logFC) >1) %>%   column_to_rownames()save(diffLab,group,file = "diffLab.Rda")##########################################################

热图绘制(对应作业6、7和10)

library(pheatmap)heatdata #筛选出差异表达明显的基因表达谱pheatmap(heatdata)#可以简单画一下，都是没有分组信息group "tumor",annotation_col group)rownames(annotation_col) #设置分组信息pheatmap(heatdata, #热图的数据         annotation_col =annotation_col, #标注样本分类         annotation_legend=TRUE, # 显示注释         show_rownames = F,# 显示行名         scale = "row", #以行来标准化，使结果更加明显         color = colorRampPalette(c("green", "white", "red"))(100),         cellwidth = 20, cellheight = 0.2,# 设置单元格大小         fontsize = 10)

学习思考

——在实战中的问题中(不同分割符号的分割问题)

这是笔者的师兄在做lncRNA数据分析时分享的一个问题，当时师兄向我分享了先把分割符变成相同的分割符，然后直接切割这个方法。由于笔者对于separate函数的用法产生兴趣，help()了一下，发现原来不同的分割符也可以"一步到位"，只需要设置sep的参数即可。

所以学习R语言是一个循序渐进的过程，不可能一蹴而就，不断学习，不断提升(共勉！！！)

请多多指教！！！

编辑：骆栢维

校审：梁晓杰

相关阅读：

学习健明老师发布的R语言练习题的学习笔记(一)

利用ROC曲线寻找最佳cutoff值(连续型变量组成的riskscore)

如何使用X-tile软件寻找最佳cutoff值

m6A评分可作为胃癌的预后标志物并预测免疫治疗的效果，IF=10.679，实时IF=15.83

为何predict()函数计算的Riskscore不等于基因的表达量与其系数的乘积的加权呢?