R语言初级作业10题
2019-07-05 本文已影响0人
DrKu
打开 Rstudio 告诉我它的工作目录
getwd()
新建6个向量，基于不同的原子类型。（重点是字符串，数值，逻辑值）
a <- c('apple','pear','hello')
b <- c(3,5,7,9)
c <- c(F,T)
告诉我在你打开的rstudio里面 getwd() 代码运行后返回的是什么？
getwd()
新建一些数据结构，比如矩阵，数组，数据框，列表等重点是数据框，矩阵）
mat <- matrix(1:20,4,5)
array1 <- array(1:24,c(2,3,4))
patientid <- c(1,2,3,4)
age <- c(23,44,56,78)
etilogy <- c('CAD','MI','HF','CAD')
sur <- c(T,F,T,F)
EF <- c(25,45,22,40)
HF <- data.frame(patientid,age,etilogy,sur,EF)
myhead <-'My <- list'
g <- c(23,4,5,6)
h <- matrix(1:10,2)
k <- c("allen",'pear','watermelon')
mylist <- list(myhead=myhead,age_sg=g,h,k)
在你新建的数据框进行切片操作，比如首先取第1，3行， 然后取第4，6列
HF[1:3,]
HF[,4:5]
使用data函数来加载R内置数据集 rivers 描述它。并且可以查看更多的R语言内置的数据集https://mp.weixin.qq.com/s/dZPbCXccTzuj0KkOL7R31g
data(rivers)
rivers
下载 https://www.ncbi.nlm.nih.gov/sra?term=SRP133642 里面的 RunInfo Table 文件读入到R里面，了解这个数据框，多少列，每一列都是什么属性的元素。（参考B站生信小技巧获取runinfo table） 这是一个单细胞转录组项目的数据，共768个细胞，如果你找不到RunInfo Table 文件，可以点击下载，然后读入你的R里面也可以。
SRP <- read.table("SRP133642.txt",header = T,sep='\t')
dim(SRP)
class(SRP$BioSample)
mode(SRP$BioSample)
ncol(SRP)
1:ncol(SRP)
class(SRP[,1])
mode(SRP[,2])
for (i in 1:ncol(SRP)) {
x=class(SRP[,i])
print(x)
}
下载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，也可以点击下载。
load("sample.Rdata")
for (i in 1:ncol(pdata)) {
x=class(pdata[,i])
print(x)
}
把前面两个步骤的两个表（RunInfo Table 文件，样本信息sample.csv）关联起来， 使用merge函数。
m1 <- merge(pdata,a,by.y = 'Sample_Name',by.x = 'geo_accession')
基于下午的统计可视化。对前面读取的 RunInfo Table(SRP) 文件在R里面探索其MBases列，包括 箱线图(boxplot)和五分位数(fivenum)，还有频数图(hist)，以及密度图(density) 。
SRP$MBases
boxplot(SRP$MBases,main='boxplot')
plot(fivenum(SRP$MBases),main='fivenum')
hist(SRP$MBases,main = 'hist')
plot(density(SRP$MBases),main='density')
把前面读取的样本信息表格的样本名字根据下划线分割看第3列元素的统计情况。第三列代表该样本所在的plate
load(file="sample.Rdata")
Title=pdata$title
Title=as.character(pdata$title)
class(Title)
strsplit(Title[1],'_')[[1]][3]
length(Title)
plate <- unlist(lapply(Title,function(x){
x
strsplit(x,'_')[[1]][3]
}))
plate
根据plate把关联到的 RunInfo Table 信息的MBases列分组检验是否有统计学显著的差异。
a <- read.table('SRP133642.txt',header = T,sep='\t')
save(a,file='RunInfo.Rdata')
m1 <- merge(pdata,a,by.y = 'Sample_Name',by.x = 'geo_accession')
t.test(m1$MBases~plate)
分组绘制箱线图(boxplot)，频数图(hist)，以及密度图(density)
boxplot(m1$MBases~plate,main='boxplot')
e <- m1[,c("title",'MBases')]
e$plate <- plate
hist(e$MBases,main='boxplot',freq = F, breaks = "sturges")
plot(density(e$MBases))
使用ggplot2把上面的图进行重新绘制。
library(ggplot2)
e$num <- 1:768
colnames(e)
ggplot(e,aes(x=plate,y=MBases)) + geom_boxplot()
# 频数图
ggplot(e,aes(x=MBases)) + geom_histogram(fill="lightblue",colour="grey") + facet_grid(plate ~ .) # 频数图
ggplot(e,aes(x=MBases,fill=plate))+geom_histogram()
ggplot(e,aes(x=MBases,fill=plate))+geom_histogram()
# 密度图
ggplot(e,aes(y=MBases,x=num)) + geom_point() + stat_density2d(aes(alpha=..density..),geom = "raster",contour = F)+ facet_grid(plate ~ .)
ggplot(e,aes(x=MBases,fill=plate))+geom_density()
使用ggpubr把上面的图进行重新绘制。
# install.packages("ggpubr")
library(ggpubr)
ggboxplot(e, x="plate", y="MBases", color = "plate", palette = "aaas",add = "jitter") + stat_compare_means(method = "t.test")
gghistogram(e, x="MBases", fill = "plate",palette = c("#f4424e", "#41a6f4"))
ggdensity(e, x="MBases", fill = "plate", color = "plate", add = "mean",palette = c("#f4424e", "#41a6f4"))
随机取384个MBases信息，跟前面的两个plate的信息组合成新的数据框，第一列是分组，第二列是MBases，总共是384*3行数据。
sample(e$MBases,384)
nrow(e)
data <- e[sample(nrow(e),384),]
data <- data[,c(3,2,4,1)]
str(data)
