R语言初级作业10题

打开 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)