Day15-R學(xué)習(xí)-生信人應(yīng)該這樣學(xué)R語(yǔ)言-筆記(轉(zhuǎn)載)

01. 介紹r語(yǔ)言及rstudio編輯器

安裝r及rstudio
打開rstudio編輯器:4塊
1.新建腳本,markdown, 寫source代碼,點(diǎn)任何地方會(huì)進(jìn)行下一步。
2.命令在console運(yùn)行,初學(xué)一步一步運(yùn)行代碼,上下左右鍵看歷史命令記錄。
3.history為歷史記錄,可以to source 或 to console進(jìn)行重新運(yùn)行,environment環(huán)境變量,賦值的變量會(huì)出現(xiàn)在這個(gè)位置。
4.files查看有哪些文件或?yàn)g覽文件
plot顯示畫出的圖,dev.off()關(guān)閉畫板
packages顯示安裝的所有包,.libPaths()顯示包安裝路徑
Help 幫助文檔查看運(yùn)行函數(shù)example,體會(huì)如何使用該函數(shù)

定位當(dāng)前文件位置:getwd()

02. R語(yǔ)言基礎(chǔ)變量講解

變量類型
**1.向量vector **
a=c(1,2,3)
b=c(1,'2',3)
class() 查看變量類型(a為numeric,b為character)

可以用函數(shù)創(chuàng)建向量,或直接使用內(nèi)置變量,左邊變量,右邊值。
數(shù)字不加引號(hào),字符串加引號(hào)(*單引號(hào)或雙引號(hào)一般可通用,當(dāng)值內(nèi)本身含有單引號(hào)時(shí),需用雙引號(hào)。例如:a='Hello World!'a="Hello World!"輸出內(nèi)容相同.b="This's a nice example."b='This's a nice example.'不同。)

2.矩陣matix
向量加維度為矩陣 dim()加維度,矩陣中任何一個(gè)元素類型發(fā)生變化,其他元素的類型也會(huì)變化。
取元素的方式:1.通過下標(biāo)來取,#逗號(hào)左邊是行右邊是列 a[,]; 2.邏輯符取
判斷元素類型: class(), str()

3.數(shù)組array

4.data.frame()數(shù)據(jù)框
轉(zhuǎn)變其中一個(gè)元素的類型,其他元素的類型不變

is.系列函數(shù),as.系列函數(shù)
是什么函數(shù)和轉(zhuǎn)變?yōu)槭裁春瘮?shù)(用Table鍵補(bǔ)全時(shí)可看到很多函數(shù))

提問方式
新建文件夾,放入需要提問/報(bào)錯(cuò)的“.Rproj”以及代碼".R"標(biāo)注報(bào)錯(cuò)的位置,發(fā)壓縮包。
如果其中有變量,報(bào)錯(cuò)之前一步寫
save(filename,file='input.Rdata')
load(file='input.Rdata')
然后注釋掉load之前的代碼(加上#)

5.列表list
用'$'從列表中取出的是元素,數(shù)據(jù)框中取出的是一列
length()有多少元素,lapply()對(duì)每個(gè)元素操作,返回list,用unlist()讓返回值為character,用as.numeric()將其轉(zhuǎn)變?yōu)閿?shù)值型。

取下標(biāo)

  1. 用下標(biāo)取索引
    b[,]逗號(hào)左行右列。若沒有該元素,返回NA;若有,返回對(duì)應(yīng)元素
  2. 用判斷方式取下標(biāo)
    b[,c(F,F,T,F,F,)]取b第三列的元素

read.table()函數(shù)讀取table
grep()#搜索函數(shù),例grep('RNA-Seq',a$Assay_Type)搜索在a的Assay_Type一列中,搜索含有RNA-Seq的下標(biāo)。
grepl()#取匹配到的所有的行,返回TRUE/FALSE,例grepl('RNA-Seq',a$Assay_Type)

取list
list[]取元素,不能取元素里的內(nèi)容,拿出的元素可能還是一個(gè)list,如果用[[*]][*]取出的是元素里面的內(nèi)容

03. 外部數(shù)據(jù)導(dǎo)入導(dǎo)出

讀入
1.直接用import dataset
2.代碼讀取
read.table('filename',sep='\t',header = T)#參數(shù)與參數(shù)之間用逗號(hào)分開,sep指定制表符,header加上表頭。

write.csv()#讀出為csv文件

取b第一列作為b的行名,然后去掉第一列
row.names(b)=b[,1]
b=[,-1]

保存為R讀取的格式,避免讀入讀出時(shí)格式的更改
save(b,file = 'b_input.Rdata')
load(file = 'b_input.Rdata')

04. 中級(jí)變量操作

讀入數(shù)據(jù)
read.table(filename,header = T,sep = '\t')
read.table(tablename,comment.char = "!",header = T,sep = '\t')#手動(dòng)指定注釋符,以!開頭的不讀入comment.char = "!"

讀出csv文件,并且去掉行名兩種方式
write.csv(b,'tmp.csv',row.names=F)
write.table(b,'tmp.csv',sep = ',')

英文單詞為函數(shù),函數(shù)有參數(shù),互通
sort()#分類
max()#最大值
min()#最小值
fivenum()#包括最小值, 25%分位數(shù), 中位數(shù), 75%分位數(shù), 最大值

table()#查看多少個(gè)元素
boxplot(y~x)#分組取域值,可用boxplot查看峰位數(shù)
先分組,分別賦值,"="是賦值,"=="是判斷,
rna=a[a$Assay_type=='Rna-Seq',]#a中Assay_type這一列中Rna-Seq對(duì)應(yīng)的行
wxs=a[a$Assay_type=='WXS',]#a中Assay_type這一列中WXS對(duì)應(yīng)的行

表達(dá)矩陣,不同的表達(dá)量
str(mean(b[1,]))#取b的第一行平均值報(bào)錯(cuò),查看其類型
mean(as.numeric(b[1,]))#轉(zhuǎn)換為數(shù)值型取均值
head(rowMeans(b))#對(duì)所有的行取均值,查看前10個(gè)

取b中行的均值
rowMeans(b)#直接用函數(shù)
for (i in 1:nrow(b)){ print(mean(as.numeric(b[i,])))}#用循環(huán)
apply(b,1,function(x){mean(x)})#apply函數(shù)
apply(b,1,mean)#與上一個(gè)相同

有函數(shù)可直接用函數(shù),沒有函數(shù)可以自定義后進(jìn)行操作
rowMax=function(x){apply(x,1,max)}#定義一個(gè)函數(shù)
rowMax(b)#用定義的函數(shù)進(jìn)行操作

向量取元素直接寫下標(biāo),data.frame()list()取元素[,]左邊取行右邊取列

top50的方差對(duì)應(yīng)的基因名做熱圖
cg=names(sort(apply(b,1,sd),decreasing = T)[1:50])
library(pheatmap)
pheatmap(b[1:50,])
pheatmap(sample(1:nrow(b),50))#隨機(jī)選取50個(gè)方差

05. 熱圖

不同大小的數(shù)據(jù)達(dá)標(biāo)不同顏色的深淺
rnorm(n, mean = 0, sd = 1)#隨機(jī)正態(tài)分布

畫個(gè)熱圖

a1=rnorm(100)
#dim加維度
dim(a1)=c(5,20)
a2=rnorm(100)+2
dim(a2)=c(5,20)
#as.data.frame將向量變?yōu)榫仃?b=as.data.frame(cbind(a1,a2))
#paste函數(shù)給a1和a2取名并以此作為b的名字看熱圖的變化,
names(b)=c(paste('a1',1:20,sep = '_'),paste('a2',1:20,sep = '_'))
#加載pheatmap包
library(pheatmap)
#加注釋信息
tmp=data.frame(group=c(rep('a1',20),rep('a2',20)))
rownames(tmp)=colnames(b)
pheatmap(b,cluster_cols = F) #cluster_cols = F 不進(jìn)行排序

pheatmap學(xué)習(xí),查看example
?pheatmap

06. 選取差異名明顯的基因的表達(dá)矩陣?yán)L熱圖

#均一化(normalization)調(diào)整差異過大的數(shù),轉(zhuǎn)置-scale-轉(zhuǎn)置。
n=t(scale(t(dat[cg,])))
n[n>2]=2#大于二的值均等于二
n[n<-2]=-2#小于負(fù)二的值均等于負(fù)二

?scale()#查看scale函數(shù)的用法

若代碼過長(zhǎng),可以定義一個(gè)函數(shù)將代碼包裹起來,后面如果有需要可以直接用不重復(fù)寫代碼
d_h <- function(dat,group_list){*}
d_h(dat,group_list)

07. id轉(zhuǎn)換

#矩陣a的V1列元素含有以'.'分割的內(nèi)容,想要去除‘.’后的內(nèi)容
library(stringr)#載入stringr
str_split(a$V1,'[.]',simplify = T)[,1]#用str_split()函數(shù),simplify=T返回值為字符型矩陣,F(xiàn)返回為字符型向量,[,1]取第一列

library(org.Hs.eg.db)#加載org.Hs.eg.db包
toTable()#org.Hs.eg.db包規(guī)定讀入數(shù)據(jù)用toTable

#兩個(gè)數(shù)據(jù)框做關(guān)聯(lián),兩個(gè)數(shù)據(jù)框共有相同的列/行名,保留沒有關(guān)聯(lián)上的元素 
merge(a,b,by='*',all.x=T)
#出現(xiàn)的頻率,對(duì)大于1的做統(tǒng)計(jì)
table(table(*)>1)
table(*)[table(*)>1]
#dataframe(d)的行去重
d=d[!duplicated(),]
#a所在的順序放在b這邊來
match(a,b)

08. 任意基因任意癌癥表達(dá)量分組的生存分析

一個(gè)基因在TCGA數(shù)據(jù)庫(kù)中的各個(gè)癌癥的生存分析
網(wǎng)頁(yè)中輸入感興趣的基因,分組(高低表達(dá)量),下載數(shù)據(jù),
Rstudio中讀入數(shù)據(jù),ggstatsplot包中的代碼進(jìn)行分析

a=read.table('*.csv',header = T,sep = ',',fill = T)
ggbetweenstats(data=a,x='Group',y='Expression')

添加Group和Expression最好復(fù)制加'',避免打錯(cuò)和不識(shí)別。
學(xué)會(huì)基礎(chǔ)變量及函數(shù)操作,然后對(duì)包的說明書進(jìn)行學(xué)習(xí),數(shù)據(jù)應(yīng)用

09. 任意基因任意癌癥表達(dá)量和臨床性狀關(guān)聯(lián)

網(wǎng)頁(yè)中輸入感興趣的基因,下載數(shù)據(jù),導(dǎo)入Rstudio
a=read.table('*.txt',header = T,sep = '\t',fill = T)
列名賦予簡(jiǎn)單的值,dat與之前寫過的代碼中一致,不需改動(dòng)代碼可直接運(yùn)行(個(gè)人習(xí)慣)

colnames(a)=c('id','stage','gene','mut')
dat=a

10.表達(dá)矩陣的樣本的相關(guān)性

cor() 關(guān)聯(lián)函數(shù)

#先安裝Bioconductor的代碼,需要R版本3.5.0
if (!requireNamespace("BiocManager"))
    install.packages("BiocManager")
BiocManager::install()

Bioconductor三個(gè)包:airway(數(shù)據(jù)包),Annotatiobdbi(注釋包)和GenomicFeatures(功能函數(shù)包)

#加載數(shù)據(jù),獲取表達(dá)矩陣,查看表達(dá),用dim()維度函數(shù)獲得矩陣行列情況
data(airway)
exprSet=assay(airway)
colnames(exprSet)
dim(exprSet)
#樣本與樣本之間的相關(guān)性,
cor(exprSet[,1],exprSet[,2])

相關(guān)性很高警惕兩種情況
1.同一個(gè)樣本技術(shù)重復(fù)
2.大多數(shù)值為低表達(dá)值或零,掩蓋真實(shí)相關(guān)性
cor(exprSet)#整個(gè)樣本做相關(guān)性

#篩選數(shù)據(jù),例如一個(gè)基因表達(dá)量大于一樣本量少于5,用dim()維度函數(shù)獲得矩陣行列情況,可以與之前的dim()做比較
x= exprSet[1,]
table(x>1)
true=1,false=0
sum(x>1)>5
exprSet= exprSet[apply(exprSet,1,sum(x>1)>5),1]
#繼續(xù)篩選,edgeR::cpm()函數(shù)去除文庫(kù)大小差異,mad的前50個(gè)變化最大的基因
#最后cor做分析,畫熱圖
M=cor(log2(exprSet+1))

#畫heatmap 加注釋信息,先分組(group_list)并轉(zhuǎn)為一個(gè)數(shù)據(jù)框,加注釋信息(多少列不重要,行名等于矩陣的列名很重要)tmp=data.frame(g=group_list)
rownames(tmp)=colnames(M)
library(pheatmap)
pheatmap(cor(exprSet),annotation_col = tmp)

11.芯片表達(dá)矩陣下游分析

構(gòu)建表達(dá)矩陣,分組信息

#單個(gè)做差異分析
t.test(exprSet[1,]~ group_list)
boxplot(exprSet[1,]~ group_list)
#或用limma()包,構(gòu)造一個(gè)design矩陣,矩陣中顯示分組信息屬于(值為1)不屬于(值為0)

12.RNA-seq表達(dá)矩陣差異分析

airway()
assay()
表達(dá)矩陣,分組準(zhǔn)備,用DEseq2()
統(tǒng)計(jì)學(xué)相關(guān)的知識(shí):statquest

原文

最后編輯于
?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請(qǐng)聯(lián)系作者
【社區(qū)內(nèi)容提示】社區(qū)部分內(nèi)容疑似由AI輔助生成,瀏覽時(shí)請(qǐng)結(jié)合常識(shí)與多方信息審慎甄別。
平臺(tái)聲明:文章內(nèi)容(如有圖片或視頻亦包括在內(nèi))由作者上傳并發(fā)布,文章內(nèi)容僅代表作者本人觀點(diǎn),簡(jiǎn)書系信息發(fā)布平臺(tái),僅提供信息存儲(chǔ)服務(wù)。

友情鏈接更多精彩內(nèi)容