原帖链接: 生信菜鸟团
R语言基础知识的一些检验，最好是对照几本R基础语法书籍来理解。
全部答案及视频在：https://github.com/jmzeng1314/R_bilibili
首先做完了周末班工作， 包括软件安装以及R包安装：http://www.bio-info-trainee.com/3727.html
打开 Rstudio 告诉我它的工作目录。
getwd()
新建6个向量，基于不同的原子类型。（重点是字符串，数值，逻辑值）
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() 代码运行后返回的是什么？
当前工作目录
新建一些数据结构，比如矩阵，数组，数据框，列表等。(重点是数据框，矩阵）
myvector <- c(1,2,3,4) #包括数据向量、逻辑向量、字符串向量
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)
在你新建的数据框进行切片操作，比如首先取第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函数来加载R内置数据集 rivers 并描述它。并且可以查看更多的R语言内置的数据集：https://mp.weixin.qq.com/s/dZPbCXccTzuj0KkOL7R31g
data("rivers")
class(rivers)
？rivers
下载 https://www.ncbi.nlm.nih.gov/sra?term=SRP133642 里面的 RunInfo Table 文件读入到R里面，了解这个数据框，多少列，每一列都是什么属性的元素（参考B站生信小技巧获取runinfo table 这是一个单细胞转录组项目的数据，共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里面，了解这个数据框，多少列，每一列都是什么属性的元素。（参考 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）关联起来，使用merge函数。
Merged_ab <- merge(a, b, by.x = "Accession",by.y = "Sample_Name")
基于下午的统计可视化
对前面读取的 RunInfo Table 文件在R里面探索其MBases列，包括: 箱线图(boxplot)和五分位数(fivenum)，还有频数图(hist)，以及密度图(density) 。
fivenum(b$MBases)
par(mfrow=c(1,3)) #在一张图上显示
boxplot(MBases ~ MBytes, data = b)
plot(density(b$MBases))
hist(b$MBases)
显示结果：
> fivenum(b$MBases)
[1] 0 8 12 16 74
Rplot.png
把前面读取的样本信息表格的样本名字 根据下划线分割 看第3列元素的统计情况。第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
根据plate把关联到的 RunInfo Table 信息的MBases列分组检验是否有显著的统计学差异。
t.test(d[,1]~plate)$p.value
分组绘制箱线图(boxplot)，频数图(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))
Boxplot-grouped.png
density plot.png
Histograms.png
使用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")
ggpubr_boxplot.png
ggpubr_density.png
ggpubr_hist.png
References：
简书作者[Forest_Lee] 盘一盘 生信技能树R语言小作业（初级）
网友评论