R语言作业—初级

教程对应B站：【生信技能树】生信人应该这样学R语言
配套资料：B站的11套生物信息学公益视频配套讲义、练习题及思维导图
先仔细观看视频，理解代码含义

初级题目

题目链接：http://www.bio-info-trainee.com/3793.html

R语言基础知识

1、找到工作目录
getwd()

2、数据类型：数值、字符串、逻辑、复数

num <- c(1,2,3,4,5,6)
cha <- c("hello world")
log <- c(TRUE,FALSE)
com <- c(1 + 0i)

3、数据结构类型：向量，列表，矩阵，数组，因子，数据帧

c <- c(1,2,3)
list <- list(c(1,2,3),'3,14')
M = matrix( c('a','a','b','c','b','a'), nrow = 2, ncol = 3, byrow = TRUE)
a <- array(c('a','b'),dim = c(3,3,2))
f <- factor(c)
df <- data.frame(gene=paste0("gene",1:15),
                 s1=rnorm(n=15),s2=rnorm(n=15),s3=rnorm(n=15),s4=rnorm(n=15),s5=rnorm(n=15))
# 取df的第1，3行，取第4，6列
df[c(1,3),]
df[,c(4,6)]

4、使用data函数来加载R内置数据集 rivers，其他数据集：https://mp.weixin.qq.com/s/dZPbCXccTzuj0KkOL7R31g

data()
#Lengths of Major North American Rivers
rivers  
head(rivers)
tail(rivers)

# 数据集的长
length(rivers)
# structure 显示对象内部结构
str(rivers) 
# 获取描述性统计量（最小值/最大值/四分位数/数值型变量/因子向量/逻辑型向量）
summary(rivers) 
rm(list = ls())

5、下载 https://www.ncbi.nlm.nih.gov/sra?term=SRP133642 里面的RunInfo Table文件读入到R里面，了解这个数据框，多少列，每一列都是什么属性的元素。 （参考B站生信小技巧获取runinfo table: https://www.bilibili.com/video/av25131640/?p=7）

初级-RunInfo Table图示-1.png 初级-Runinfo Table图示-2.png

这是一个单细胞转录组项目的数据，共768个单细胞转录组的测序数据，如果你找不到RunInfo Table文件，可以点击下载，然后读入你的R里面也可以。

# 读取RunInfo Table,防止默认读取为Factors
options(stringsAsFactors = F)
# 有报错信息，加上了一些参数就可以了
Table <-read.table(file="SraRunTable.txt",header = T,sep = '\t',fileEncoding = "utf-8")
str(Table)
## 'data.frame':    768 obs. of  31 variables:
## ......
colSums(Table=="yes")

colnames(Table)
##  [1] "BioSample"          "Experiment"         "MBases"            
##  [4] "MBytes"             "Run"                "SRA_Sample"        
##  [7] "Sample_Name"        "Assay_Type"         "AssemblyName"      
## [10] "AvgSpotLen"         "BioProject"         "Center_Name"       
## [13] "Consent"            "DATASTORE_filetype" "DATASTORE_provider"
## [16] "InsertSize"         "Instrument"         "LibraryLayout"     
## [19] "LibrarySelection"   "LibrarySource"      "LoadDate"          
## [22] "Organism"           "Platform"           "ReleaseDate"       
## [25] "SRA_Study"          "age"                "cell_type"         
## [28] "marker_genes"       "source_name"        "strain"            
## [31] "tissue"

6、下载 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，也可以点击下载。

初级-samples图示-1.png 初级-samples图示-2.png 初级-samples图示-3.png

# 读取样本信息
sample <-read.csv("sample.csv")
colnames(sample)
##  [1] "Accession"           "Title"               "Sample.Type"        
##  [4] "Taxonomy"            "Channels"            "Platform"           
##  [7] "Series"              "Supplementary.Types" "Supplementary.Links"
## [10] "SRA.Accession"       "Contact"             "Release.Date"

7、两个表格关联(RunInfo Table文件和样本信息文件sample.csv)

dim(Table)
## [1] 768  31
dim(sample)
## [1] 768  12
d = merge(Table,sample,by.x = "Sample_Name",by.y = "Accession")
# merge() 函数 此时Table的Sample_Name列和sample的Accession列内容相同，合并这一列
dim(d)
## [1] 768  42

统计可视化基础

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

# Mbases样本的碱基数
# 举例：
boxplot(Table$MBases, main = "boxplot of MBases")
plot(fivenum(Table$MBases), main = "fivenum of MBases")
hist(Table$MBases, main = "hist of MBases")
plot(density(Table$MBases,na.rm=T), main = "density of MBases")

初级-boxplot of MBases.png 初级-hist of MBases.png 初级-density of MBases.png 初级-fivenum of MBases.png

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

# lapply 'l'代表list 将指定操作应用于列表中的所有元素
# 字符串的分割函数，指定分隔符，生成list
title = sample$Title
class(title)
## [1] "character"

plate = unlist(lapply(title,function(x){
  x
  strsplit(x,'_')[[1]][3]
}))
# strsplit(title[1],'_')[[1]][3]  结果为"0048"
# 这里写了一个函数，对每个title值进行上一行的操作，获得768个值

# table() 函数各个值出现的频率
table(plate) 
## 0048 0049 
##  384  384

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

# plate 指384孔PCR板，编号分别是48号和49号

t.test(d$MBases~plate)
##  Welch Two Sample t-test
## 
## data:  d$MBases by plate
## t = 2.3019, df = 728.18, p-value = 0.02162
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  0.1574805 1.9831445
## sample estimates:
## mean in group 0048 mean in group 0049 
##           13.08854           12.01823
# 显著性检验相关知识还需要学习...

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

boxplot(d$MBases~plate)
typeof(plate)
e = d[,c("MBases","Title")]
e$plate = plate
hist(e$MBases,freq = F, breaks = "sturges")
plot(density(e$MBases,na.rm=T))
# 比较简陋，可以尝试用ggplot2和ggpubr画图包

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

suppressMessages(library(ggplot2))
e$plate = plate
e$num=c(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(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()

初级-箱线图.png
初级-频数图-1.png 初级-频数图-2.png
初级-密度图-1.png 初级-密度图-2.png

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

suppressMessages(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"))

初级-ggpubr-箱型图.png 初级-ggpubr-频数图.png 初级-ggpubr-密度图.png

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

# sample() 函数 随机抽样
data <- e[sample(nrow(e),384),][,c(3,1,2)]
str(data)
## 'data.frame':    384 obs. of  3 variables:
##  $ plate : chr  "0049" "0048" "0049" "0049" ...
##  $ MBases: int  3 16 5 2 8 14 25 11 16 7 ...
##  $ Title : chr  "SS2_15_0049_J6" "SS2_15_0048_N2" "SS2_15_0049_M5" "SS2_15_0049_N24" ...

更多学习资源：
生信技能树公益视频合辑
生信技能树简书账号
生信工程师入门最佳指南
...
欢迎关注生信菜鸟团、生信技能树！