【生信技能树】R语言初级练习题

作者: 函数IsLazy | 来源:发表于2020-03-27 22:12 被阅读0次

2019-04-19
生信人的20个R语言习题--by me
【生信技能树】R语言初级练习题
R语言视频笔记
R语言作业—中级
R语言作业-20题
R语言作业—初级
【生信技能树::作业&习题】R语言初级练习题
学习小组Day4笔记--镰羲
【生信技能树】R语言练习题 - 高级

getwd()

新建6个向量，基于不同的原子类型。（重点是字符串，数值，逻辑值）

a = c(1,2,3,4,5,6)
b = c("a","b","c","d","e","f","g")
c = c(T,F,T,F)

新建5个其它数据结构，矩阵，数组，数据框，列表，因子（重点是数据框，矩阵）

 ## c一定要小写，大写C会报错：object not interpretable as a factor
mymatrix = matrix(1:20, nrow = 5, ncol = 4)
myarray = array(1:24, c(2,3,4))
mydatef = data.frame(PatientID = c("pa01","pa02","pa03"),
                   age = c(22,33,44),
                   exp = c(12,23,34))
mylist = list(mymatrix, mydatef, a, b, c)
status = c("good", "moderate", "poor")
myfactor = factor(status, levels = c("poor", "moderate", "good"))
class(mymatrix)
class(myarray)
class(mydatef)
class(mylist)
class(myfactor)

在你新建的数据框进行切片操作，比如首先取第1，3行，然后取第4，6列

mydatef[c(1,3),]
mydatef[,c(2,3)] #列数不够就随便取了

使用data函数来加载R内置数据集 rivers 描述它。并且可以查看更多的R语言内置的数据集：

data("rivers")
class(rivers) #数据类型"numeric"
str(rivers)  #紧凑地显示对象内部结构，即对象里有什么：num [1:141] 735 320 325 392 524 ...
length(rivers)  # 数据的宽度
summary(rivers)  # 可以提供最小值、最大值、四分位数和数值型变量的均值，以及因子向量和逻辑型向量的频数统计
head(rivers); tail(rivers) # 数据的前6个  后6个数据

下载 https://www.ncbi.nlm.nih.gov/sra?term=SRP133642 里面的 RunInfo Table 文件读入到R里面，了解这个数据框，多少列，每一列都是什么属性的元素。

SRA（Sequence ReadArchive，https://www.ncbi.nlm.nih.gov/sra/ ）数据库是专门用于存储二代测序的原始数据，包括 454, IonTorrent, Illumina, SOLiD, Helicos and CompleteGenomics等。除了原始序列数据外，SRA现在也存在raw reads在参考基因的aligment information。

数据下载参考解读SRA数据库规律一文就够
 https://www.ncbi.nlm.nih.gov/Traces/study/?acc=SRP133642 下载样本ID表格--新版本的已经换了名字，如下图

image.png

读入数据并检查

rm(list = ls())
options(stringsAsFactors = F)
sra = read.table("SraRunTable.txt",header = T,sep = ',')
dim(sra) # 查看data.frame维度  [1] 768  30
ncol(sra) #查看data.frame列数View(sra) #用表格形式观察data.frame
#查看每一列都是什么属性的元素
install.packages("magrittr") 
install.packages("dplyr")    # alternative installation of the %>%
library(magrittr) # need to run every time you start R and want to use %>%
library(dplyr)    # alternative, this also loads %>%
for (i in colnames(sra)) paste(i,class(sra[,i])) %>% print()

%>%是管道函数，就是把左件的值发送给右件的表达式，并作为右件表达式函数的第一个参数，不加载那两个包会报错-----抄的猫叽先森的做法。

下载 https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE111229 里面的样本信息读入到R里面，了解这个数据框，多少列，每一列都是什么属性的元素
方法1

library(GEOquery)
GSE_name = 'GSE111229'
options( 'download.file.method.GEOquery' = 'libcurl' ) 
gset = getGEO( GSE_name, getGPL = F )
save( gset, file = 'gset.Rdata' )
load("gset.Rdata")
Gset = gset[[1]]      
pdata = pData(Gset)
class(pdata)  #[1] "data.frame"
View(pdata) 
dim(pdata)  #[1] 768  47
colnames(pdata)

方法2

image.png

image.png
一定要选All search results不然只会下载你页面看到的20个左右

pd = read.csv('sample.csv')
class(pdata)   # [1] "data.frame"
dim(pd)    # [1] 768  12

两种方法得到的观测数是一样的，变量数第一种方法更多，相应的信息量更大。

##查看每列的数据类型
library(magrittr) 
library(dplyr) 
for (i in colnames(sra)) paste(i,class(sra[,i])) %>% print()

把前面两个步骤的两个表（RunInfo Table 文件，样本信息sample.csv）关联起来，使用merge函数。

#查看两个data.frame里相关的列
rm(list = ls())
options(stringsAsFactors = F)
sra <- read.table(file = "SraRunTable.txt",header = T,sep = ',')
geo<- read.csv(file = "sample.csv")
library(magrittr)   #需要用到 %>%

for (i in colnames(sra)) {
    if (i %in% colnames(geo)) print(i)
}
#这里得到的是相同列名的数据名 [1] "Platform"，查看一下platform的数据
head(sra$Platform)   #[1] "ILLUMINA" "ILLUMINA" "ILLUMINA" "ILLUMINA" "ILLUMINA"
[6] "ILLUMINA"
head(geo$Platform)   # [1] "GPL13112" "GPL13112" "GPL13112" "GPL13112" "GPL13112"
[6] "GPL13112"
#列名相同但是内容完全不同。

#然后考虑在数据里查看相同内容：
for (i in colnames(sra)) {
   for (j in colnames(geo)) {
    if (sra[,i] %in% geo[,j]) {paste(i,j) %>% print()}
      }
  }
# [1] "Experiment SRA.Accession"
  [1] "GEO_Accession..exp. Accession"
  [1] "Organism Taxonomy"
  [1] "Sample.Name Accession"
#以上四组，就是两个data.frame内容一致的列，任选一对为标准合并data.frame
df_merge = merge(sra,geo,by.x="Experiment",by.y="SRA.Accession")
#df_merge = merge(sra,geo,by.x="Sample_Name",by.y="Accession")
#df_merge = merge(sra,geo,by.x="Organism",by.y="Taxonomy")
#df_merge = merge(sra,geo,by.x="GEO_Accession..exp.",by.y="Accession")
dim(df_merge)  #[1] 768  41
dim(geo)    #[1] 768  12
dim(sra)    #[1] 768  30
# 41=30+12-1

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

函数概念及图形美化：
箱线图 - boxplot
五分位数 - fivenum
频数图 - histogram
密度图 - density

rm(list = ls())
options(stringsAsFactors = F)
sra <- read.table(file = "SraRunTable.txt",header = T,sep = ',')
dim(sra)
#箱线图
boxplot(sra$Bases)
#五分位数
fivenum(sra$Bases)
#频数图
hist(sra$Bases)
#密度图 用于估计随机变量概率密度函数的一种非参数方法。是一种连续型变量分布的有效方法
density(sra$Bases)
plot(density(sra$Bases))

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

rm(list = ls())
options(stringsAsFactors = F)
pd = read.csv('sample.csv')
pd$plate=unlist(lapply(pd$Title,function(x){
  strsplit(x,"[_]")[[1]][3]
})
)
# 用一个lapply函数提取出title分割后的第三列，然后在数据集中添加一列$plate
pd$plate
table(pd$plate)  #0048 0049    384  384

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

根据之前的题目对数据进行合并

sra = read.table(file = "SraRunTable.txt",header = T,sep = ',')
df_merge = merge(sra,geo,by.x="Experiment",by.y="SRA.Accession")
t.test(df_merge$Bases~df_merge$plate) # p=0.02256

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

# 箱线图
boxplot(df_merge$Bases~df_merge$plate)
# 按照plate建立两个新的数据集
df0048 = df_merge[df_merge$plate=='0048',]
df0049 = df_merge[df_merge$plate=='0049',]
hist(df0048$Bases)
hist(df0049$Bases)
#密度图
plot(density(df0048$Bases))
plot(density(df0049$Bases))

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

library(ggplot2)
# 箱线图
ggplot(df_merge,aes(x=plate, y= Bases))+
  geom_boxplot(fill= 'cornflowerblue',
               color="black")
# 频数图
ggplot(data=df_merge, aes(x=Bases ,fill=plate)) + 
  geom_histogram(bins=30,color= "black")+ 
  facet_wrap(~ plate)

image.png

# 密度图
ggplot(data=df_merge, aes(x=Bases ,fill=plate)) + 
  geom_density(alpha= .3)

image.png

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

library(ggpubr)
ggboxplot(data = df_merge,
          x = "plate", y = "Bases",
          fill = "plate")
gghistogram(data = df_merge,x = "Bases",
            fill = "plate")
ggdensity(data = df_merge,x = "Bases",
          fill = "plate")