[ZHUHAI_Biotrainee] R语言初级作业（下）

基于下午的统计可视化

1、对前面读取的 RunInfo Table 文件在R里面探索其MBases列，包括 箱线图(boxplot)和五分位数(fivenum)，还有频数图(hist)，以及密度图(density) 。

e在a中探索MBases列，重新载入e：

rm(list = ls())

a <- read.csv('SraRunTable.txt',header = T,sep = "\t")

b <- read.csv('sample.csv',header = T,sep = ",")d <- merge(a,b,by.x = 'Sample_Name',by.y = 'Accession')

e=d[,c("MBases","Title")]

e[,2]

class(e[,2])

e[,2] <- as.character(e[,2])

class(e[,2])

plate=unlist(lapply(e[,2],function(x){ # x=e[1,2] x strsplit(x,'_')[[1]][3] }))

箱线图、五分位数、频数图、密度图

boxplot(e[,1])
fivenum(e[,1])
[1]  0  8 12 16 74
hist(e[,1])
plot(density(e[,1]))

2、把前面读取的样本信息表格的样本名字根据下划线分割看第3列元素的统计情况。第三列代表该样本所在的plate

class(e[,2])
plate=unlist(lapply(e[,2],function(x){
  # x=e[1,2]
  x
  strsplit(x,'_')[[1]][3]
  
}))

Result: 未完待续... 部分

当我们看到这样的结果时，不禁思考：
如何把plate中的“0048”、“0049”全部取出来呢？这是需要解决的。对于前几步的densty的作图中我们也可以将0048与0049分开。不过，对于初学者我们今天这道题先到这里，后续的0048与0049分开的部分留到之后（本页面的第4题）再进行解决。

3、根据plate把关联到的 RunInfo Table 信息的MBases列分组检验是否有统计学显著的差异。

t.test(e[,1]~plate)

4、分组绘制箱线图(boxplot)，频数图(hist)，以及密度图(density) 。

分组箱图

boxplot(e[,1]~plate)

Resul

分组频数图（hist）

e$plate=plate
plate0 <- e[,c("MBases","plate")]
plate1 <- plate0[1:384,1]
hist(plate1)          #此步为0048的频数图


plate2<- plate0[385:768,1]
hist(plate2)   #此步为0049的频数图

Result:

0048_hist 0049_hist

分组密度图

density(plate1)
plot(density(plate1))    #此步为0048的密度图
plot(density(plate2))   #此步为0049的密度图

5、使用ggplot2把上面的图进行重新绘制

library(ggplot2)
colnames(e)
ggplot(e,aes(x=plate,y=MBases))+geom_boxplot()
ggplot(e,aes(x=plate,y=MBases))+geom_histogram()   #geom_histogram() 不完整，会报错，需要把（）内的部分填充
plate1 <- as.data.frame(plate1)
plate2 <- as.data.frame(plate2)            #必须转化成数据框才可以作图，先以plate1为例，colname为plate1,其实这里是全部0048的MBases
colnames(plate1)="MBases"   
colnames(plate2)="MBases" 
ggplot(plate1,aes(x=MBases))+geom_histogram(bins = 60,binwidth=0.01,color="blue")
ggplot(plate1,aes(x=MBases))+geom_histogram(bins = 30,color="blue")
#可以前后对比，发现前后差别。  
ggplot(plate2,aes(x=MBases))+geom_histogram(bins = 60,binwidth=0.01,color="blue")
ggplot(plate2,aes(x=MBases))+geom_histogram(bins = 30,color="blue")

#密度图
ggplot(plate1,aes(x=MBases))+geom_density()
ggplot(plate2,aes(x=MBases))+geom_density()

Result:
结果请参考：“E:\生信技能树\Qi_study(1 month)\7.3初级作业\ggplot_hist_0048-0049、hist_0048、hist_0049、density_0048、density_0049

值得注意有以下两点：
①我们在作图时利用了“ggplot2"
②在作图一时，使用boxplot()函数；作图二时叠加了一个geom_boxplot函数，如果我们在以后的作图中还想叠加其他的函数或者添加相应的操作可以在上一个函数结束后用“+”进行得加。>
Eg:
ggplot(e,aes(x=plate,y=MBases))+geom_boxplot()

6、使用ggpubr把上面的图进行重新绘制。

library(ggpubr)
p <- ggboxplot(e, x = "plate", y = "MBases",
               color = "plate", palette = "jco",
               add = "jitter")
plot(p)
p + stat_compare_means(method = 't.test')
#下面代码为hist_0048
ggplot(plate0,aes(x=MBases))+geom_histogram(bins = 40,binwidth=0.01,color="blue")

结果见file.name: E:\生信技能树\Qi_study(1 month)\7.3初级作业\“ggboxplot","p+1"

show_ggboxplot\p+1结果在file

7、随机取384个MBases信息，跟前面的两个plate的信息组合成新的数据框，第一列是分组，第二列是MBases，总共是384*3行数据。

plate3 <- sample(1:nrow(e),384)

colnames(plate1) <- "MBases_0048"
colnames(plate2) <- "MBases_0049"

data1 <- data.frame(plate2,plate1)
data2 <- data.frame(plate3,data1)

心得：
在此处合并时，建议采用data.frame的形式。而不是merge 可以试试merge的格式，但是merge一般都是都是有相同的列或者行才会用的。这里是没有的。要学会用灵活思考问题、解决问题！