非常有幸在入门生信的路上遇到了生信技能树,在Jimmy老师的指导下终于找到了方向,今后就是埋头苦干的阶段啦。
做完高级题后才回来整理初级题(http://www.bio-info-trainee.com/3793.html)的代码,加了一些新东西,对原来没有完全理解的参考代码也有了更深的理解。
文末有彩蛋!
rm(list = ls())
options(stringsAsFactors = F)
a=c(1,2,3) #数值
b=c('a','b','c') #字符串
c=c(T,F,T) #逻辑值
str(a)
str(b)
str(c) #确认
getwd() #返回当前工作空间
mat=matrix(a,nrow = 3)
df=data.frame(name = b, age = a)
pts = list(x = cars[,1], y = cars[,2]) #来自list()的help文档
df$x=1
df$y=1:3 #之前的数据框太小了,随便加两列
df[c(1,3),]
df[,c(2,4)] #类矩阵操作
df[c(2,4)]
df$age
df$y #其他的方法
data(rivers)
class(rivers)
length(rivers)
str(rivers)
head(rivers) #用这些函数试图了解rivers,发现它是一个长度为141的数值型向量
#打开链接后,Send results to Run selector - RunInfo Table即可
runinfo = read.table("SraRunTable.txt", header = T, sep = '\t')
#打开链接后,More - Accession list truncated, click here to browse through all related public accessions即可
saminfo = read.csv('sample.csv')
dim(saminfo)
str(saminfo)
colnames(saminfo) #这是一个12列,768行的数据框。行名是样本名,列名是样本信息类型
d=merge(runinfo,saminfo,by.x = 'Sample_Name',by.y = 'Accession') #根据runinfo的Sample_Name和saminfo的Accession合并
e=d[,c("MBases","Title")] #由于后续只需要用到MBases和Title这两个信息,可以将暂时不用的数据弃去
save(e,file = 'input.Rdata')
rm(list = ls())
options(stringsAsFactors = F)
load(file = 'input.Rdata')
boxplot(e$MBases)
plot(fivenum(e$MBases))
hist(e$MBases)
plot(density(e$MBases))
# e[,2]
# plate=unlist(lapply(e[,2],function(x){
# # x=e[1,2]
# x
# strsplit(x,'_')[[1]][3]
# }))
#以上为答案中的写法
lis=strsplit(e$Title, "_")
lis[[1]][3]
lapply(strsplit(e$Title, "_"), function(x){x[3]}) #写下面的语句之前先对每一部分进行摸索
plate = unlist(lapply(strsplit(e$Title, "_"), function(x){x[3]})) #利用apply对每一行进行操作
table(plate) #查看分组情况
e$plate=plate #将分组信息加入到数据中
t.test(e$MBases~plate)
boxplot(e[,1]~plate)
par(mfrow=c(2,1),mar=c(2,2,2,2)) #调整画板参数
hist(e[which(e$plate=='0048'),1],main = '0048')
hist(e[which(e$plate=='0049'),1],main = '0049')
plot(density(e[which(e$plate=='0048'),1]),main = '0048')
plot(density(e[which(e$plate=='0049'),1]),main = '0049')
Rplot1.png
Rplot2.png
Rplot3.png
library(ggplot2)
colnames(e)
ggplot(e,aes(x=plate,y=MBases))+geom_boxplot()
suppressPackageStartupMessages(library(ggpubr))
p <- ggboxplot(e, x = "plate", y = "MBases",
color = "plate", palette = "jco",
add = "jitter")
p + stat_compare_means(method = 't.test') #在图上添加p-value信息
sp = e[sample(1:768, size = 384),]
data.frame(group = sp$plate, MBases = sp$MBases, Title = sp$Title)
Rplot4.png
Rplot5.png
最后,向大家隆重推荐生信技能树的一系列干货!
- 生信技能树超级VIP入场券:https://mp.weixin.qq.com/s/K8WiX5C7BYC3WonM2VcdHQ
- B站公益74小时生信工程师教学视频合辑:https://mp.weixin.qq.com/s/IyFK7l_WBAiUgqQi8O7Hxw
网友评论