原帖鏈接: 生信菜鳥團
R語言基礎(chǔ)知識的一些檢驗兼砖，最好是對照幾本R基礎(chǔ)語法書籍來理解语婴。
全部答案及視頻在：https://github.com/jmzeng1314/R_bilibili
首先做完了周末班工作儒飒， 包括軟件安裝以及R包安裝：http://www.bio-info-trainee.com/3727.html
打開 Rstudio 告訴我它的工作目錄颈将。
getwd()
新建6個向量锻霎，基于不同的原子類型嫉称。（重點是字符串茬暇，數(shù)值霹期，邏輯值）
a <- c(1,2,3,4)
b <- c("hello","world","!")
d <- c(FALSE, TRUE, T, F) #注意大小寫
e <- c(1+0i, 2+4i) #complex
f <- c(1:4) # integer; or f <- c(1L,2L,3L,4L)
print(class(f))
g <- charToRaw("Hello"); # raw(?)
print(class(g))
> [1] "raw"
告訴我在你打開的rstudio里面 getwd() 代碼運行后返回的是什么？
當前工作目錄
新建一些數(shù)據(jù)結(jié)構(gòu)快压，比如矩陣圆仔，數(shù)組，數(shù)據(jù)框蔫劣，列表等坪郭。(重點是數(shù)據(jù)框，矩陣）
myvector <- c(1,2,3,4) #包括數(shù)據(jù)向量脉幢、邏輯向量截粗、字符串向量
myfactor <- c('green','green','yellow','red','red','red','green')
mymatrix <- matrix( c('a','a','b','c','b','a'), nrow = 2, ncol = 3, byrow = TRUE)
myarray <- array(c('green','yellow'),dim = c(3,3,2))
mydataframe <- data.frame(
gender = c("Male", "Male","Female"),
height = c(152, 172, 165),
weight = c(81,93, 78),
Age = c(42,38,26)
)
# list
a <- "My List"
b <- c(25, 26, 18, 39)
c <- matrix(1:10, nrow=5)
d <- c("one", "two", "three")
mylist <- list(title=a ,b,c,d)
在你新建的數(shù)據(jù)框進行切片操作信姓，比如首先取第1，3行绸罗， 然后取第2，4列
mydf1 <- mydataframe[c(1,3),]
mydf2 <- mydf1[,c(2,4)]
# method 2
mydf2 <- mydataframe[c(1,3),c(2,4)]
使用data函數(shù)來加載R內(nèi)置數(shù)據(jù)集 rivers 并描述它豆瘫。并且可以查看更多的R語言內(nèi)置的數(shù)據(jù)集：https://mp.weixin.qq.com/s/dZPbCXccTzuj0KkOL7R31g
data("rivers")
class(rivers)
珊蟀？rivers
下載 https://www.ncbi.nlm.nih.gov/sra?term=SRP133642 里面的 RunInfo Table 文件讀入到R里面，了解這個數(shù)據(jù)框外驱，多少列育灸，每一列都是什么屬性的元素（參考B站生信小技巧獲取runinfo table 這是一個單細胞轉(zhuǎn)錄組項目的數(shù)據(jù)，共768個細胞）昵宇。如果你找不到RunInfo Table 文件磅崭，可以點擊下載，然后讀入你的R里面也可以瓦哎。
a <- read.csv("sample.csv")
dim(a)
str(a)
下載 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111229 里面的樣本信息sample.csv讀入到R里面砸喻，了解這個數(shù)據(jù)框，多少列蒋譬，每一列都是什么屬性的元素割岛。（參考 https://mp.weixin.qq.com/s/fbHMNXOdwiQX5BAlci8brA 獲取樣本信息sample.csv）如果你實在是找不到樣本信息文件sample.csv，也可以點擊下載犯助。
b <- read.table("SraRunTable.txt", header = T,sep = "\t")
dim(b)
str(b)
把前面兩個步驟的兩個表（RunInfo Table 文件癣漆，樣本信息sample.csv）關(guān)聯(lián)起來，使用merge函數(shù)剂买。
Merged_ab <- merge(a, b, by.x = "Accession",by.y = "Sample_Name")
基于下午的統(tǒng)計可視化
對前面讀取的 RunInfo Table 文件在R里面探索其MBases列惠爽，包括: 箱線圖(boxplot)和五分位數(shù)(fivenum)，還有頻數(shù)圖(hist)瞬哼，以及密度圖(density) 婚肆。
fivenum(b$MBases)
par(mfrow=c(1,3)) #在一張圖上顯示
boxplot(MBases ~ MBytes, data = b)
plot(density(b$MBases))
hist(b$MBases)
顯示結(jié)果：
> fivenum(b$MBases)
[1] 0 8 12 16 74
把前面讀取的樣本信息表格的樣本名字 根據(jù)下劃線分割 看第3列元素的統(tǒng)計情況。第3列代表該樣本所在的plate倒槐。
d <- Merged_ab[,c("MBases", "Title")]
save(d, file = "input.Rdata")
rm(list = ls())
options(stringsAsFactors = F)
load(file = "input.Rdata")
e <- lapply(d[,2], function(x){
x
strsplit(x,"_")[[1]][3]
})
plate <- unlist(e)
#method 2
plate=unlist(lapply(as.character(d[,2]),function(x){
x
strsplit(x,'_')[[1]][3]
}))
boxplot(d[,1]~plate)
d$plate <= plate
根據(jù)plate把關(guān)聯(lián)到的 RunInfo Table 信息的MBases列分組檢驗是否有顯著的統(tǒng)計學差異旬痹。
t.test(d[,1]~plate)$p.value
分組繪制箱線圖(boxplot)，頻數(shù)圖(hist)讨越，以及密度圖(density) 两残。
# 分組畫出boxplot
d1 <- d[d$plate == '0048',]
d2 <- d[d$plate == '0049',]
par(mfrow=c(1,2))
boxplot(d1$MBases~d1$plate)
boxplot(d2$MBases~d2$plate)
# 分組畫出density
plot(density(d[d$plate == '0048',]$MBases))
plot(density(d[d$plate == '0049',]$MBases))
# 分組畫出histogram
hist((d[d$plate == '0048',]$MBases))
hist((d[d$plate == '0049',]$MBases))
使用ggplot2把上面的圖進行重新繪制。
library(ggplot2)
colnames(d)
ggplot(d,aes(x=plate,y=MBases))+geom_boxplot()
ggplot(d,aes(x=MBases))+geom_density()
ggplot(d,aes(x=MBases))+geom_histogram(bins = 50)
使用ggpubr把上面的圖進行重新繪制把跨。
library(ggpubr)
p <- ggboxplot(d, x = "plate", y = "MBases",
color = "plate", palette = "jco",
add = "jitter")
# Add p-value
p + stat_compare_means(method = 't.test')
ggdensity(d, x = "MBases",
color = "plate", palette = "jco", rug = T,
add = "mean")
#color 可以顯示分組信息人弓；rug是下面的柵欄效果；add是添加none,mean或者median線
gghistogram(d, x = "MBases",
color = "plate", palette = "jco",bins = 50,
add = "mean")
References：
簡書作者[Forest_Lee] 盤一盤 生信技能樹R語言小作業(yè)（初級）
