? ? ? ? ? ? ? ? ? ? ? ?? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? ? 學習者：駱栢維

題目來源：生信基石之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不等于基因的表達量與其系數的乘積的加權呢?

總結

以上是生活随笔為你收集整理的语言 提取列名_学习健明老师发布的R语言练习题的学习笔记（二）的全部內容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網站內容還不錯，歡迎將生活随笔推薦給好友。