guangzhou_biotrainee_taskI_junior初级

1、首先是软件和R package安装：

rm(list = ls()) 
options()$repos 
options()$BioC_mirror
options(BioC_mirror="https://mirrors.ustc.edu.cn/bioc/")
options("repos" = c(CRAN="https://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
options()$repos 
options()$BioC_mirror

以上都是镜像设置相关操作，其中的ustc.edu是设置成中科大镜像，下面的tuna.tsinghua.edu.cns是设置成清华镜像，之所以要设置镜像，是因为直接从CRAN下载R包可能速度较慢。同样的，当一个R包安装不顺利的时候，也可以考虑换个镜像去安装。
关于镜像设置，其实鼠标操作更简单：
Tools>Global Options>Packages>Change，选择中国国内镜像，一共有6个可选，一般清华的是最稳定的。
这里顺便说一下R 包安装的三种方式，

(1)从CRAN安装库：

install.packages("BiocManager")

也可以从Rstudio界面用鼠标点击下载，效果是一样的，如图：


R package install.png

(2)从Bioconductor安装：

(!requireNamespace("BiocManager", quietly = TRUE)) 
install.packages("BiocManager") 
BiocManager::install("biomaRt", version = "3.8")#安装biomaRt

(3)github上R包安装：

devtools::install_github("ggrothendieck/sqldf")#如果没有devtools需要先安装devtools,如下：
     install.packages("devtools")
     library(devtools)

大部分包我都已经安装过，只有2个需要安装：

library("FactoMineR")
library("factoextra")
library(GSEABase)
library(GSVA)
BiocManager::install('GSVA')
library(clusterProfiler)
library(genefu)#这个包不太好装，看网，win7比win10貌似好装点
if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")
BiocManager::install("genefu")

library(ggplot2)
library(ggpubr)
library(hgu133plus2.db)
library(limma)
library(org.Hs.eg.db)
library(pheatmap)

2、下面开始做题：

1、打开 Rstudio 告诉我它的工作目录。#告诉我在你打开的rstudio里面 getwd() 代码运行后返回的是什么？
getwd()
2、新建6个向量，基于不同的原子类型。（重点是字符串，数值，逻辑值）

a1<-c("a","b","c")
a2<-c(1:3)#数值型向量（integer）  
a2_2<-c(2.1,4.2,5.5,6.7,9.3,10.6)#数值型向量（小数）    
a3<-c("hello","R","nihao")#字符串向量
a4<-  a1<a2#通过逻辑运算赋值逻辑型向量
a4_2 <- c(T,F,F,F,T,F)#直接赋值逻辑值向量   
a5<-factor(c("stageI","stageII","stageIII"))
a6<-as.numeric(a2)+1
a7 <- c(1,2,3,"yes")#这里是看到别人这么写的，ps，向量不是只能有一种类型吗？
a8 <- rep('e',times=3)
a1
a2
a2_2
a3
a4
a4_2
a5
a6
a7
class(a8)
a8

3、新建一些数据结构，比如矩阵，数组，数据框，列表等重点是数据框，矩阵）

b1<-matrix(LETTERS[1:24],nrow=4,ncol = 6)
b1
b1_1<- matrix(data=c(1:20),nrow=4,ncol=5) 
b1_1
b2<-data.frame(a1,a2)
b2
b3<-matrix(1:10,nrow=5)
b3
b3_0<-c(rep(a1,2),rep(a2,2))
b3_0
b3<-as.data.frame(matrix(b3_0,ncol = 6))
b3
b4<-matrix(rep(c(a1,a2),3),nrow=3,byrow = F)
b4
b5 <- list(Name <- c("Liming","Xiaofang","Xiaohong","Lihua","Shanshan"),Gender <- c("M","F","M","M"),score <- matrix(1:12,3,4))
b5 #创建列表

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

b1
b1[1,3]
b1[,c(4,6)]#通过位置索引

colnames(b1)<-c(1:6)
b1
b1[,colnames(b1)>3&colnames(b1)!=5]#通过逻辑值索引

dim(b1)
tmp1 <- c(T,F,T,F)
tmp2 <- c(F,F,F,T,F,T)
b1[tmp1,tmp2]#通过逻辑值索引

colnames(b1)<-paste0("col_",c(1:6))
b1
b1[,c("col_4","col_6")]#通过名字索引
b1

5、使用data函数来加载R内置数据集 rivers 描述它。并且可以查看更多的R语言内置的数据集：
R语言系列：探索R自带数据包(这里列举的是可能会用到和遇到的数据，包括一些以前见过的数据包）"https://mp.weixin.qq.com/s/dZPbCXccTzuj0KkOL7R31g

data(rivers)
head(rivers)
rivers
data(euro)#欧元汇率，长度为11，每个元素都有命名
euro
data(mtcars)#32辆汽车在11个指标上的数据
data(esoph)   ##法国的一个食管癌病例对照研究    
esoph    
data(iris)     # #3种鸢尾花形态数据
data(longley)    #强共线性的宏观经济数据   
data(LifeCycleSavings)    #50个国家的存款率
data(PlantGrowth)     #三种处理方式对植物产量的影响
data(quakes)     #1000次地震观测数据（震级>4）
data(sleep)     #两药物的催眠效果
data(ToothGrowth)     #VC剂量和摄入方式对豚鼠牙齿的影响
data(ChickWeight)    #饮食对鸡生长的影响
data(CO2)    #耐寒植物CO2摄取的差异
data(Indometh)    #某药物的药物动力学数据
data(Theoph)    #茶碱药动学数据

#时间序列数据
data(ldeaths)    #1974-1979年英国每月支气管炎、肺气肿和哮喘的死亡率
data(fdeaths)    #前述死亡率的女性部分
data(mdeaths)    #前述死亡率的男性部分
data(USAccDeaths)    #1973-1978美国每月意外死亡人数
data(EuStockMarkets)    #多变量时间序列。欧洲股市四个主要指标的每个工作日记录，共1860条记录

6、下载 https://www.ncbi.nlm.nih.gov/sra?term=SRP133642 里面的 RunInfo Table 文件读入到R里面，了解这个数据框，多少列，每一列都是什么属性的元素。（参考B站生信小技巧获取runinfo table） 这是一个单细胞转录组项目的数据，共768个细胞，如果你找不到RunInfo Table 文件，可以点击下载，然后读入你的R里面也可以。
下载方法，NCBI---SRA---搜索SRP133642，点击“Send results to Run selector"，进来后再点击”RunInfo Table“即可下载。

sra<-read.table('SraRunTable.txt',header = T,sep = '\t')
dim(sra)
str(sra)
colnames(sra)

7、下载 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，也可以点击下载。
下载sample.csv:下载并保存GEO数据(第一次下载失败，第二次成功)

library(GEOquery)
GSE_name = 'GSE111229'
options( 'download.file.method.GEOquery' = 'libcurl' ) #windows系统
gset <- getGEO( 'GSE111229', destdir = '.',getGPL = F )
save( gset, file = 'gset.Rdata' ) #掌握保存和加载Rdata的命令

由于gset是列表，故将其转为可操作的数据结构Gset

load("gset.Rdata")
Gset <- gset[[1]]

用GEOquery里的pdata函数获取样本信息
看一下pdata的结构，很明显是数据框

pdata<-pData(Gset)
dim(pdata)#766行，47列，发现比手动下载的列数多
#dim查看行列 colnames查看列名

手动下载，载入数据https://www.ncbi.nlm.nih.gov/geo/browse/
两种方法都可以，这里我用的是手动下载的数据

sample<-read.csv('sample.csv',header = T,sep = ',')
head(sample)
class(sample)
View(sample)
dim(sample)
colnames(sample)

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

sra_merge<-merge(sra,sample,by.x='Sample_Name',by.y='Accession')

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

boxplot(sra$MBases)
fivenum(sra$MBases)
hist(sra$MBases)
plot(density(sra$MBases))

9、把前面读取的样本信息表格的样本名字根据下划线分割看第3列元素的统计情况。第三列代表该样本所在的plate, 根据plate把关联到的 RunInfo Table 信息的MBases列分组检验是否有统计学显著的差异

sample
head(sra_merge)
sra1<-sra_merge[,c("MBases","Title")]
head(sra1)
save(sra1,file = 'input.Rdata')
load(file = 'input.Rdata')
class(sra1)
str(sra1)
sra1$Title=as.character(sra1$Title)
sra1

strsplit(sra1$Title,"_")[[1]][3]
sra1

plate <- unlist(lapply(sra1$Title,function(x)
  {  x
  strsplit(x,'_')[[1]][3]
}))
table(plate)
class(plate)
head(plate)

str(plate)

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

sra1$plate<-plate
head(sra1)
t.test(sra1[,1] ~ plate)
[1] t = 2.3019, df = 728.18, p-value = 0.02162。有显著差异。

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

其实基础作图都可以不用花费太多时间去学习，练练玩玩即可，大部分功能都可以用ggplot来解决。
sra1$plate<-as.factor(plate)
sra1$plate<-as.numeric(plate)#这两个变量类型的转化需根据后面分组作图的要求确定要不要运行。
str(sra1$plate)
hist(sra1[,1],main='boxplot',freq = F, breaks = "sturges")
boxplot(sra1[,1]~plate)

plate0 <- sra1[,c("MBases","plate")]
plate1 <- plate0[1:384,1]
hist(plate1)          #此步为0048的频数图
plate2<- plate0[385:768,1]
hist(plate2)   #此步为0049的频数图

density(plate1)
plot(density(plate1))    #此步为0048的密度图
plot(density(plate2))   #*此步为0049的密度图#这里的plate1和plate2赋值原先没有想到，感觉稍微麻烦了点，暂时没想到更好的办法进行分组作图，当然boxplot不需要赋值即可进行分组。*

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

library(ggplot2)
colnames(sra1)
ggplot(sra1,aes(x=plate,y=MBases))+geom_boxplot()+theme_classic()#我喜欢用这个主题，背景比较干净

频数图：柱状图和直方图是很类似的，只是柱状图更适合画分类变量，
直方图更适合把连续型的数据按照一个个等长的分区（bin）来切分，然后计数，画柱状图。

sra1$plate<-as.factor(plate)

ggplot(sra1,aes(x=MBases,fill=plate))+geom_histogram()#与下面一句等价
ggplot(sra1)+geom_histogram(aes(x=MBases,fill=plate))#直方图最容易，提供一个变量x，画出数据分布，并可以根据其他变量给他填充颜色。
ggplot(sra1)+geom_histogram(aes(x=MBases,fill=plate),position = "dodge")#并列排列
ggplot(sra1)+geom_histogram(aes(x=MBases,fill=plate),position = "fill")#使用position="fill"，按照相对比例来画

画成柱状图是这个样子：

ggplot(sra1)+geom_bar(aes(x=MBases,y=1), stat="identity")#通过stat参数，可以让geom_bar按指定高度画图

密度图

ggplot(sra1,aes(x=MBases,fill=plate))+geom_density()
ggplot(sra1)+geom_density(aes(x=MBases, colour=plate,fill=plate))#color指定线条颜色，fill指定图片颜色填充

其他人的答案（boxplot一样）：(此时需要先把plate转换为numeric变量)

sra1$plate<-as.numeric(plate)
ggplot(sra1,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()

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

ggpubr号称专做出版级的图，确实做出来比较好看。
library(ggpubr)
#boxplot:
p <- ggboxplot(sra1, x="plate", y="MBases", color = "plate", 
               palette = c("#00AFBB", "#E7B800", "#FC4E07"), 
               add = "jitter", shape="plate")#增加了jitter点，点shape由plate映射
p
#换个颜色(别人答案)
p <- ggboxplot(sra1, x = "plate", y = "MBases",
               color = "plate", palette = "jco",
               add = "jitter")
# Add p-value
p + stat_compare_means(method = 't.test')
#再换个颜色
ggboxplot(sra1, x="plate", y="MBases", color = "plate", 
          palette = "aaas",add = "jitter") + stat_compare_means(method = "t.test")

#箱线图+分组形状+统计：
my_comparisons <- list(c("0.5", "1"))
p+stat_compare_means(comparisons = my_comparisons)+ #不同组间的比较 
  stat_compare_means(label.y = 60)

#内有箱线图的小提琴图+星标记：
ggviolin(sra1, x="plate", y="MBases", fill = "plate", 
         palette = c("#00AFBB", "#E7B800", "#FC4E07"), 
         add = "boxplot", add.params = list(fill="white"))+ 
  stat_compare_means(comparisons = my_comparisons, label = "p.signif")+#label这里表示选择显著性标记（星号） 
  stat_compare_means(label.y = 50)

#直方图：
sra1$plate<-as.factor(plate)
gghistogram(sra1, x="MBases", fill = "plate",palette = c("#f4424e", "#41a6f4"))
gghistogram(sra1, x="MBases", add = "mean", rug = TRUE, color = "plate", fill = "plate",
            palette = c("#00AFBB", "#E7B800"))#带有均值线和边际地毯线的直方图

#密度图：
ggdensity(sra1, x="MBases", fill = "plate",  color = "plate", add = "mean",palette = c("#f4424e", "#41a6f4"))
ggdensity(sra1, x="MBases", add = "mean", rug = TRUE, color = "plate", fill = "plate",
          palette = c("#00AFBB", "#E7B800"))#带有均值线和地毯线的密度图

其他人：

hist_0048
ggplot(plate0,aes(x=MBases))+geom_histogram(bins = 40,binwidth=0.01,color="blue")

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

这个题目没太明白意思，主要是练习sample随机取样函数的用法。
col_4<-sample(nrow(sra1), 384)#用到随机抽样函数sample
col_4<-sample(sra1$MBases,384)
sra1$col_4<-col_4
head(sra1)
sra1<-sra1[,c(3,1,4,2)]
sra1

别人的答案一：https://blog.csdn.net/weixin_34240520/article/details/91199550

sra2 <- sample(nrow(sra1), 384)
class(sra2)
sra3 <- sra1[sra2,][, c(3, 2, 1)]
dim(sra3)

别人的答案二：https://www.jianshu.com/p/448ec60bdd4a

data <-  sra1[sample(nrow(sra1),384),]
data <- data[,c(3,2,4,1)]
str(data)
dim(data)

merge一般都是都是有相同的列或者行才会用的。

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

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

data1 <- data.frame(plate2,plate1)
data2 <- data.frame(plate3,data1)
最后编辑于 ：
©著作权归作者所有,转载或内容合作请联系作者
  • 人面猴
    序言：七十年代末，一起剥皮案震惊了整个滨河市，随后出现的几起案子，更是在滨河造成了极大的恐慌，老刑警刘岩，带你破解...
    沈念sama阅读 219,490评论 6赞 508
  • 死咒
    序言：滨河连续发生了三起死亡事件，死亡现场离奇诡异，居然都是意外死亡，警方通过查阅死者的电脑和手机，发现死者居然都...
    沈念sama阅读 93,581评论 3赞 395
  • 救了他两次的神仙让他今天三更去死
    文/潘晓璐 我一进店门，熙熙楼的掌柜王于贵愁眉苦脸地迎上来，“玉大人，你说我怎么就摊上这事。” “怎么了？”我有些...
    开封第一讲书人阅读 165,830评论 0赞 356
  • 道士缉凶录：失踪的卖姜人
    文/不坏的土叔 我叫张陵，是天一观的道长。 经常有香客问我，道长，这世上最难降的妖魔是什么？ 我笑而不...
    开封第一讲书人阅读 58,957评论 1赞 295
  • 港岛之恋（遗憾婚礼）
    正文 为了忘掉前任，我火速办了婚礼，结果婚礼上，老公的妹妹穿的比我还像新娘。我一直安慰自己，他们只是感情好，可当我...
    茶点故事阅读 67,974评论 6赞 393
  • 恶毒庶女顶嫁案：这布局不是一般人想出来的
    文/花漫 我一把揭开白布。 她就那样静静地躺着，像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上，一...
    开封第一讲书人阅读 51,754评论 1赞 307
  • 城市分裂传说
    那天，我揣着相机与录音，去河边找鬼。 笑死，一个胖子当着我的面吹牛，可吹牛的内容都是我干的。 我是一名探鬼主播，决...
    沈念sama阅读 40,464评论 3赞 420
  • 双鸳鸯连环套：你想象不到人心有多黑
    文/苍兰香墨 我猛地睁开眼，长吁一口气：“原来是场噩梦啊……” “哼！你这毒妇竟也来了？” 一声冷哼从身侧响起，我...
    开封第一讲书人阅读 39,357评论 0赞 276
  • 万荣杀人案实录
    序言：老挝万荣一对情侣失踪，失踪者是张志新（化名）和其女友刘颖，没想到半个月后，有当地人在树林里发现了一具尸体，经...
    沈念sama阅读 45,847评论 1赞 317
  • 护林员之死
    正文 独居荒郊野岭守林人离奇死亡，尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 37,995评论 3赞 338
  • 白月光启示录
    正文 我和宋清朗相恋三年，在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 40,137评论 1赞 351
  • 活死人
    序言：一个原本活蹦乱跳的男人离奇死亡，死状恐怖，灵堂内的尸体忽然破棺而出，到底是诈尸还是另有隐情，我是刑警宁泽，带...
    沈念sama阅读 35,819评论 5赞 346
  • 日本核电站爆炸内幕
    正文 年R本政府宣布，位于F岛的核电站，受9级特大地震影响，放射性物质发生泄漏。R本人自食恶果不足惜，却给世界环境...
    茶点故事阅读 41,482评论 3赞 331
  • 男人毒药：我在死后第九天来索命
    文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹，春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 32,023评论 0赞 22
  • 一桩弑父案，背后竟有这般阴谋
    文/苍兰香墨 我抬头看了看天上的太阳。三九已至，却和暖如春，着一层夹袄步出监牢的瞬间，已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 33,149评论 1赞 272
  • 情欲美人皮
    我被黑心中介骗来泰国打工， 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留，地道东北人。 一个月前我还...
    沈念sama阅读 48,409评论 3赞 373
  • 代替公主和亲
    正文 我出身青楼，却偏偏与公主长得像，于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子，可洞房花烛夜当晚...
    茶点故事阅读 45,086评论 2赞 355