2023.4生信马拉松day7-R语言综合应用

options("repos" = c(CRAN="http://mirrors.tuna.tsinghua.edu.cn/CRAN/"))
if(!require(tidyr))install.packages("tidyr",update = F,ask = F)
if(!require(dplyr))install.packages("dplyr",update = F,ask = F)
if(!require(stringr))install.packages('stringr',update = F,ask = F)
if(!require(tibble))install.packages('tibble',update = F,ask = F)
library(tidyr)
library(dplyr)
library(stringr)
library(tibble)

rm(list = ls())
if(!require(stringr))install.packages('stringr')
library(stringr)

x <- "The birch canoe slid on the smooth planks."
x
### 1.检测字符串长度
str_length(x)
length(x)

### 2.字符串拆分
str_split(x," ")
x2 = str_split(x," ")[[1]];x2  #按空格拆分x，注意没有赋值就没有发生过！

y = c("jimmy 150","nicker 140","tony 152")
str_split(y," ")
str_split(y," ",simplify = T)   #简化拆分结果，简化成返回矩阵而不是列表

str_sub(x,5,9)  #提取x的第5到9个字符

str_detect(x2,"h")  #判断x的每个字符串含不含有某个字母或者多个字母的组合；
str_starts(x2,"T")    #判断是否以某个元素开头；
str_ends(x2,"e")   #判断是否以某个元素结尾；

x2
str_replace(x2,"o","A")
str_replace_all(x2,"o","A")

x
str_remove(x," ")
str_remove_all(x," ")

test <- iris[c(1:2,51:52,101:102),]
rownames(test) =NULL # 去掉行名，NULL是“什么都没有”
test

# arrange，数据框按照某一列排序
library(dplyr)
arrange(test, Sepal.Length) #从小到大
arrange(test, desc(Sepal.Length)) #从大到小

# distinct，数据框按照某一列去重复
distinct(test,Species,.keep_all = T)

mutate(test, new = Sepal.Length * Sepal.Width)
#问题：新增列之后，test这个数据框是5列还是6列（有没有发生改变）？——没有，因为没有赋值就没有发生过！！！

# 连续的步骤
# 1.多次赋值，产生多个中间的变量
x1 = select(iris,-5)  #取iris除了第5列的所有列
x2 = as.matrix(x1)
x3 = head(x2,50) #取x2的前50行
pheatmap::pheatmap(x3) #用x3画热图

# 2. 嵌套，代码不易读
pheatmap::pheatmap(head(as.matrix(select(iris,-5)),50))

# 3.管道符号传递，简洁明了
iris %>%
  select(-5) %>%
  as.matrix() %>%
  head(50) %>% 
  pheatmap::pheatmap()

i = -1
if (i<0) print('up')
if (i>0) print('up')

if(!require(tidyr)) install.packages('tidyr')  #根据一个包是否library成功来决定要不要安装这个包

# 1.读取group.csv，从第二列中提取圈出来的信息
g = read.csv("group.csv")
g2 = g[,2]
g2_s = str_split(g2," ",simplify = T);g2_s
g2_s[,4]
#老师的答案：
library(stringr)
a = read.csv("group.csv")
g = str_split(a$title," ",simplify = T)
g[,4]

# 2.如何把上一题结果中的Control和Vemurafenib改成全部小写？搜索一下
tolower(g2_s[,4])
str_to_lower(g[,4])  #两个函数均可

# 3.加载test1.Rdata,按照symbol列给ids数据框去重复,注意，要让ids数据真正发生修改。
load("test1.Rdata")
ids_d = distinct(ids,symbol,.keep_all = T)
dim(ids)
dim(ids_d)

#示例：
save(pd,exp,gpl,file = "step1output.Rdata")

rm(list = ls())
load("step1output.Rdata")

i =1
if (i>0){
  print('+')
} else {
  print("-")
}

i = 1
ifelse(i>0,"+","-")

x = rnorm(3)
x
ifelse(x>0,"+","-")

samples = c("tumor1","tumor2","tumor3","normal1","normal2","normal3")
k1 = str_detect(samples,"tumor");k1
ifelse(k1,"tumor","normal")

k2 = str_detect(samples,"normal");k2
ifelse(k2,"normal","tumor")

i = 0
if (i>0){
  print('+')
} else if (i==0) {
  print('0')
} else if (i< 0){
  print('-')
}

ifelse(i>0,"+",ifelse(i<0,"-","0"))  #嵌套式的ifelse：不符合大于零的条件，就再进行一步判断；

# 1.加载deg.Rdata,根据a、b两列的值，按照以下条件生成向量x：
#a< -1 且b<0.05,则x对应的值为down；
#a>1 且b<0.05,则x对应的值为up；
#其他情况，x对应的值为no；
#统计up、down、no各出现了多少次

#我的答案：
rm(list = ls())
load("deg.Rdata")
dim(deg)
colnames(deg)
k1 = deg$a
k2 = deg$b
k3=ifelse(k1 < -1,ifelse(k2 < 0.05,"down","no"),ifelse(k1 > 1,ifelse(k2 < 0.05,"up","no"),"no"))
table(k3)

#老师的答案：
load("deg.Rdata")
k1 = deg$a< -1 & deg$b<0.05;table(k1)
k2 = deg$a>1 & deg$b<0.05;table(k2)s
x = ifelse(k1,"down",ifelse(k2,"up","no"))
table(x)

for( i in 1:4){
  print(paste0("the current number is ",i))
}

#批量画图
par(mfrow = c(2,2))  #把画板分成2*2的四块使其能够同时放四张图；
for(i in 1:4){
  print(paste0("the current column is ",colnames(iris)[i]))
  plot(iris[,i],col = iris[,5])
}

#批量装包
pks = c("tidyr","dplyr","stringr")
for(g in pks){
  print(paste0("package ",g," is installing"))
  if(!require(g,character.only = T))  #表示a是变量名，代表包
    install.packages(g,ask = F,update = F) #避免被一个需要更新的包卡住后面的代码

}

# 生成一个表达矩阵
set.seed(10086)  #为了让模拟分析的结果可重现，给rnorm设计一个随机数种子，保证它每次生成的随机数都是那一组；
exp = matrix(rnorm(18),ncol = 6) #新建一个矩阵，把18个随机数分成6列；
exp = round(exp,2)  #取小数点后两位，不写的话默认取整
rownames(exp) = paste0("gene",1:3)
colnames(exp) = paste0("test",1:6)
exp[,1:3] = exp[,1:3]+1 #给exp的1-3列加一
exp

#处理数据

library(tidyr)
library(tibble)
library(dplyr)
dat = t(exp) %>%   #转置
  as.data.frame() %>%   #变成数据框
  rownames_to_column() %>%   #把行名变成一列，
  mutate(group = rep(c("control","treat"),each = 3))
#宽变长函数的最新版本：
pdat = dat%>% 
  pivot_longer(cols = starts_with("gene"), #把gene开头的列转换掉
               names_to = "gene",  #新的列名叫gene
               values_to = "count")  #由原来的值转换得到的那一列的列名叫count

#数据转换好就可以画图了
library(ggplot2)
p = ggplot(pdat,aes(gene,count))+
  geom_boxplot(aes(fill = group))+
  theme_bw()
p
p + facet_wrap(~gene,scales = "free") #还可以根据gene分面

rm(list = ls())
test<- iris[1:6,1:4]
apply(test, 2, mean)
apply(test, 1, sum)

#等价于
colMeans(test)
rowSums(test)

par(mfrow = c(2,2))
apply(iris[,1:4], 2, plot)
#如果有写不下的参数，可以继续写在apply的括号里
apply(iris[,1:4], 2, plot,col = iris[,5])

#或者也可以自定义函数
jimmy <- function(g){
  plot(g,col = iris[,5])
}
par(mfrow = c(2,2))
apply(iris[,1:4], 2, jimmy)

lapply(1:4,rnorm)
#批量画图
lapply(1:4, function(i){
  plot(iris[,i],col = iris[,5])
})
#批量装包
pks = c("tidyr","dplyr","stringr")
qa = function(g){
  print(paste0("package ",g," is installing"))
  if(!require(g,character.only = T))
    install.packages(g,ask = F,update = F)
}
lapply(pks, qa)

test1 <- data.frame(name = c('jimmy','nicker','Damon','Sophie'), 
                    blood_type = c("A","B","O","AB"))
test1
test2 <- data.frame(name = c('Damon','jimmy','nicker','tony'),
                    group = c("group1","group1","group2","group2"),
                    vision = c(4.2,4.3,4.9,4.5))
test2
library(dplyr)
inner_join(test1,test2,by="name")
right_join(test1,test2,by="name")
full_join(test1,test2,by="name")
semi_join(test1,test2,by="name")
anti_join(test1,test2,by="name")

#我的答案：
load("test2.Rdata")
test = as.data.frame(test)
test$gene_var <- apply(test,1,var)  #计算每行的方差并添加到test
library(tidyr)
library(tibble)
test = rownames_to_column(test,var="gene_name") 
test = arrange(test,gene_var)
tail(test$gene_name,1000)

#老师的答案
load("test2.Rdata")
x=names(tail(sort(apply(test,1,var)),1000)) #看到就一行的时候我好气呀……

#1. 加载test1.Rdata，将两个数据框按照probe_id列连接在一起，按共同列取交集
#2. 找出logFC最小的10个基因和logFC最大的10个基因(symbol列就是基因名)

#我的答案：
rm(list = ls())
load("test1.Rdata")
library(dplyr)
i = inner_join(dat,ids,by="probe_id")

min10 = head(sort(i$logFC),10)
max10 = tail(sort(i$logFC),10)
min10_gene = i$symbol[i$logFC %in% min10]
max10_gene = i$symbol[i$logFC %in% max10]
min10_gene
max10_gene

#老师的答案：
rm(list = ls())
load('test1.Rdata')
library(dplyr)
x = merge(dat,ids,by = "probe_id")
x2 = inner_join(dat,ids,by = "probe_id")
dim(x)
dim(x2)

x = arrange(x,logFC)
head(x$symbol,10)
tail(x$symbol,10)
head(x$logFC,10)