R_basics_Biotrainee

1.R脚本

rm(list = ls())
#gc()
options(stringsAsFactors = F)
library(data.table)
library(tidyr)
library(dplyr)
library(pwr)
library(ggplot2)
library(ggpubr)
#get the working directory
getwd()
#新建6个向量，基于不同的“原子类型”
v1 <- letters[1:6];class(v1)
v2 <- c(1.2,3.4,5,6.7);class(v2)
v3 <- 1L:6L;class(v3)
v4 <- c(T,F,T,T,F,F);class(v4)
v5 <- 2i + 1;class(v5)
v6 <- charToRaw("Biotrainee");class(v6) 

#新建一些数据结构,数据框data1，矩阵data2，三维数组data3，列表data4
data1 <- data.frame(c1 = paste0("gene",1:4),
                    c2 = 1:4,
                    c3 = 2:5,
                    c4 = 3:6,
                    row.names = paste0("r",1:4));data1
class(colnames(data1))
data2 <- matrix(seq(24,2,-2),nrow = 4,dimnames = list(paste0("r",seq(4)),paste0("c",seq(3))));data2
data3 <- array(1:12,c(2,3,2));data3
data4 <- list(data1,data2,data3);data4

#get the subset of data frame & matrix
(sub1_data1 <- data1[1:2,])
(sub2_data1 <- data1[,2:3])
(sub1_data2 <- data2[1:2,])
(sub2_data2 <- data2[,2:3])

#加载内置数据集`rivers`并加以描述
d <- rivers;head(d,20)
class(d)
length(d)

#统计描述,极值、中位数、四分位数等信息
(summary(d))

#GSE111229
name_folder <- "RunInfoTable"
name_fold <- "SraRunTable.txt"
runinfo <- data.frame(fread(paste(name_folder,name_fold,sep = "/")))
dim(runinfo);colnames(runinfo)

#The class of each column in `runinfo`
sapply(runinfo,class)

#Read `sample`
sample <- data.frame(fread("sample.csv"))
dim(sample);colnames(sample)

#The class of each column in `sample`
sapply(sample,class)

#merge
merge_GSE111229 <- merge(runinfo,sample,by.x = "Sample_Name",by.y = "Accession")


#########################################
##############Visualization##############
#########################################
mbs <- runinfo$MBases
#Boxplot of MBases
boxplot(mbs)
#fivenum of MBases
fivenum(mbs)
#Hist of MBases
hist(mbs,col = "orange",ylab = "freq of MBases",
     main = "MBases of RunInfoTable--GSE111229")
#Density of MBases
density(mbs)
plot(density(mbs))

######################################
###############分组统计################
######################################
#merge_GSE111229内拆分
sep_merge <- separate(merge_GSE111229,Title,c("sep1","sep2","plate","sep3"),sep = "_",remove = F)
select_merge <- select(sep_merge,c("MBases","plate"))#简化数据框，只取两列
colnames(select_merge)
#查看分组
table(sep_merge$plate)##根据plate，分为两组，数量相同，各占一半

#显著性水平为0.05，以plate分组，检验MBases是否有显著差异
t.test(select_merge$MBases~select_merge$plate)
#结论：p=0.02<0.05,拒绝原假设，即MBases两组数据均值有显著差异

###############分组作图################
boxplot(select_merge$MBases~select_merge$plate)
#ggplot2
ggplot(select_merge,aes(x=plate,y=MBases))+geom_boxplot()
#ggpubr
ggboxplot(select_merge, x = "plate", y = "MBases", color = "plate", palette = "jco", add = "jitter")

#随机取384条
samp_merge <- sample_frac(select_merge,0.5)
#与前面plate合并成新数据框
comb_sample <- rbind(samp_merge,select_merge);nrow(comb_sample)

2.运行结果及解析

运行结果-part1

> rm(list = ls())
> #gc()
> options(stringsAsFactors = F)
> library(data.table)
> library(tidyr)
> library(dplyr)
> library(pwr)
> library(ggplot2)
> library(ggpubr)
> #get the working directory
> getwd()
[1] "E:/mingliya/Biotrainee/R_basics"
> #新建6个向量，基于不同的“原子类型”
> v1 <- letters[1:6];class(v1)
[1] "character"
> v2 <- c(1.2,3.4,5,6.7);class(v2)
[1] "numeric"
> v3 <- 1L:6L;class(v3)
[1] "integer"
> v4 <- c(T,F,T,T,F,F);class(v4)
[1] "logical"
> v5 <- 2i + 1;class(v5)
[1] "complex"
> v6 <- charToRaw("Biotrainee");class(v6) 
[1] "raw"
> #新建一些数据结构,数据框data1，矩阵data2，三维数组data3，列表data4
> data1 <- data.frame(c1 = paste0("gene",1:4),
+                     c2 = 1:4,
+                     c3 = 2:5,
+                     c4 = 3:6,
+                     row.names = paste0("r",1:4));data1
      c1 c2 c3 c4
r1 gene1  1  2  3
r2 gene2  2  3  4
r3 gene3  3  4  5
r4 gene4  4  5  6
> class(colnames(data1))
[1] "character"
> data2 <- matrix(seq(24,2,-2),nrow = 4,dimnames = list(paste0("r",seq(4)),paste0("c",seq(3))));data2
   c1 c2 c3
r1 24 16  8
r2 22 14  6
r3 20 12  4
r4 18 10  2
> data3 <- array(1:12,c(2,3,2));data3
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

> data4 <- list(data1,data2,data3);data4
[[1]]
      c1 c2 c3 c4
r1 gene1  1  2  3
r2 gene2  2  3  4
r3 gene3  3  4  5
r4 gene4  4  5  6

[[2]]
   c1 c2 c3
r1 24 16  8
r2 22 14  6
r3 20 12  4
r4 18 10  2

[[3]]
, , 1

     [,1] [,2] [,3]
[1,]    1    3    5
[2,]    2    4    6

, , 2

     [,1] [,2] [,3]
[1,]    7    9   11
[2,]    8   10   12

> #get the subset of data frame & matrix
> (sub1_data1 <- data1[1:2,])
      c1 c2 c3 c4
r1 gene1  1  2  3
r2 gene2  2  3  4
> (sub2_data1 <- data1[,2:3])
   c2 c3
r1  1  2
r2  2  3
r3  3  4
r4  4  5
> (sub1_data2 <- data2[1:2,])
   c1 c2 c3
r1 24 16  8
r2 22 14  6
> (sub2_data2 <- data2[,2:3])
   c2 c3
r1 16  8
r2 14  6
r3 12  4
r4 10  2

运行结果--part2

#加载内置数据集`rivers`并加以描述
> d <- rivers;head(d,20)
 [1]  735  320  325  392  524  450 1459  135  465  600  330  336  280  315  870  906  202  329  290 1000
> class(d)
[1] "numeric"
> is.vector(d)
[1] TRUE
> length(d)
[1] 141
> #统计描述,极值、中位数、四分位数等信息
> (summary(d))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  135.0   310.0   425.0   591.2   680.0  3710.0 
> #GSE111229
> name_folder <- "RunInfoTable"
> name_fold <- "SraRunTable.txt"
> runinfo <- data.frame(fread(paste(name_folder,name_fold,sep = "/")))
> dim(runinfo)
[1] 768  31       
> #The class of each column in `runinfo`
> sapply(runinfo,class)
         BioSample         Experiment             MBases             MBytes                Run 
       "character"        "character"          "integer"          "integer"        "character" 
        SRA_Sample        Sample_Name         Assay_Type       AssemblyName         AvgSpotLen 
       "character"        "character"        "character"        "character"          "integer" 
        BioProject        Center_Name            Consent DATASTORE_filetype DATASTORE_provider 
       "character"        "character"        "character"        "character"        "character" 
        InsertSize         Instrument      LibraryLayout   LibrarySelection      LibrarySource 
         "integer"        "character"        "character"        "character"        "character" 
          LoadDate           Organism           Platform        ReleaseDate          SRA_Study 
       "character"        "character"        "character"        "character"        "character" 
               age          cell_type       marker_genes        source_name             strain 
       "character"        "character"        "character"        "character"        "character" 
            tissue 
       "character" 

> #Read `sample`
> sample <- data.frame(fread("sample.csv"))
> dim(sample)
[1] 768  12
> #The class of each column in `sample`
> sapply(sample,class)
          Accession               Title         Sample.Type            Taxonomy            Channels 
        "character"         "character"         "character"         "character"           "integer" 
           Platform              Series Supplementary.Types Supplementary.Links       SRA.Accession 
        "character"         "character"         "character"         "character"         "character" 
            Contact        Release.Date 
        "character"         "character" 

> #######merge#######
> merge_GSE111229 <- merge(runinfo,sample,by.x = "Sample_Name",by.y = "Accession")

运行结果--part3

> ###############Visualization############
> mbs <- runinfo$MBases
> #Boxplot of MBases
> boxplot(mbs)
> #fivenum of MBases
> fivenum(mbs)
[1]  0  8 12 16 74
> #Hist of MBases
> hist(mbs,col = "orange",ylab = "freq of MBases",
+      main = "MBases of RunInfoTable--GSE111229")
> #Density of MBases
> density(mbs)

Call:
    density.default(x = mbs)

Data: mbs (768 obs.);   Bandwidth 'bw' = 1.423

       x                y            
 Min.   :-4.269   Min.   :0.0000000  
 1st Qu.:16.366   1st Qu.:0.0000353  
 Median :37.000   Median :0.0003001  
 Mean   :37.000   Mean   :0.0121039  
 3rd Qu.:57.634   3rd Qu.:0.0142453  
 Max.   :78.269   Max.   :0.0665647  
> plot(density(mbs))

> ###############分组统计################
> #merge_GSE111229内拆分
> sep_merge <- separate(merge_GSE111229,Title,c("sep1","sep2","plate","sep3"),sep = "_",remove = F)
> select_merge <- select(sep_merge,c("MBases","plate"))#简化数据框，只取两列
> colnames(select_merge)
[1] "MBases" "plate" 
> #查看分组
> table(sep_merge$plate)##根据plate，分为两组，数量相同，各占一半

0048 0049 
 384  384 

> #显著性水平为0.05，以plate分组，检验MBases是否有显著差异
> t.test(select_merge$MBases~select_merge$plate)

    Welch Two Sample t-test

data:  select_merge$MBases by select_merge$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 

> #结论：p=0.02<0.05,拒绝原假设，即MBases两组数据均值有显著差异
> ###############分组作图################
> boxplot(select_merge$MBases~select_merge$plate)
> #ggplot2
> ggplot(select_merge,aes(x=plate,y=MBases))+geom_boxplot()
> #ggpubr
> ggboxplot(select_merge, x = "plate", y = "MBases", color = "plate", palette = "jco", add = "jitter")
> #随机取384条
> samp_merge <- sample_frac(select_merge,0.5)
> #与前面plate合并成新数据框
> comb_sample <- rbind(samp_merge,select_merge);nrow(comb_sample)
[1] 1152

p.s.如有误望大家不吝赐教，提出错误和建议。作者：Leah 如需转载，请注明出处。