这是生信技能树关于生信人的20个R语言习题的答案：
1 安装R包
数据包： ALL, CLL, pasilla, airway
软件包：limma，DESeq2，clusterProfiler
工具包：reshape2
绘图包：ggplot2
#################################

1.安装R包

#################################
source(“https://bioconductor.org/biocLite.R”)
options(BioC_mirror=“http://mirrors.ustc.edu.cn/bioc/”)
biocLite(c(“ALL”,“CLL”, “pasilla”, “airway”)) #数据包
biocLite(c(“limma”,“DESeq2”, “clusterProfiler”)) #软件包
install.packages(“reshape2”) #工具包
install.packages(“ggplot2”) #绘图包
#另外可以检测某个包是否存在，只有不存在时才会安装
if (! require (‘CLL’)){
options(BioC_mirror=“http://mirrors.ustc.edu.cn/bioc/”)
BiocInstaller::biocLite(‘CLL’,ask = F, suppressUpdates = T)
}
#意思就是设定安装镜像，然后隐藏安装信息，中间过程不询问

biocLite(“hgu95av2.db”)
library(hgu95av2.db)
ls(“package:hgu95av2.db”)#产生36个映射数据（探针id转为36种主流id）
capture.output(hgu95av2()) #将hgu95av2包含的具体信息输出为字符串【只用hgu95av2()可能输出，但不易整理】
#结果可以看到，每一个子集都是有keys-value构成的，就相当于list的结构。因此对于hgu95av2的操作都要使用as.list()
as.list(hgu95av2SYMBOL[1]) #变成列表进行查看第一组元素
#结果得到的就是
#$1000_at
#[1] “MAPK3”
#不懂hgu95av2SYMBOL是什么意思，就查找
?hgu95av2SYMBOL
#结果得到：hgu95av2SYMBOL is an R object that provides mappings between manufacturer identifiers and gene abbreviations【就是基因id与探针号的关系，这里显示的就是MAPK3这个基因对应的探针id是1000_at】

ids=toTable(hgu95av2SYMBOL)
length(unique(idssymbol))tail(sort(table(idssymbol)) tail(sort(table(idssymbol))tail(sort(table(idssymbol)))
table(sort(table(idssymbol)))plot(table(sort(table(idssymbol))) plot(table(sort(table(idssymbol)))plot(table(sort(table(idssymbol))))

table(rownames(exprSet) %in% idsprobeid)exprSet=exprSet[rownames(exprSet)probe_id) exprSet=exprSet[rownames(exprSet) %in% idsprobeid)exprSet=exprSet[rownames(exprSet)probe_id,]
dim(exprSet)

ids=ids[match(rownames(exprSet),idsprobeid),]head(ids)exprSet[1:5,1:5]tmp=by(exprSet,idsprobe_id),] head(ids) exprSet[1:5,1:5] tmp = by(exprSet,idsprobeid),]head(ids)exprSet[1:5,1:5]tmp=by(exprSet,idssymbol,function(x) rownames(x)[which.max(rowMeans(x))] )
probes = as.character(tmp)
exprSet=exprSet[rownames(exprSet) %in% probes ,]
dim(exprSet)

rownames(exprSet)=ids[match(rownames(exprSet),idsprobeid),2]exprSet[1:5,1:5]library(reshape2)exprSetL=melt(exprSet)colnames(exprSetL)=c(′probe′,′sample′,′value′)exprSetLprobe_id),2] exprSet[1:5,1:5] library(reshape2) exprSet_L=melt(exprSet) colnames(exprSet_L)=c('probe','sample','value') exprSet_Lprobeid),2]exprSet[1:5,1:5]library(reshape2)exprSetL=melt(exprSet)colnames(exprSetL)=c(′probe′,′sample′,′value′)exprSetLgroup=rep(group_list,each=nrow(exprSet))
head(exprSet_L)

ggplot2

library(ggplot2)
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print§
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_violin()
print§
p=ggplot(exprSet_L,aes(value,fill=group))+geom_histogram(bins = 200)+facet_wrap(~sample, nrow = 4)
print§
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()+facet_wrap(~sample, nrow = 4)
print§
p=ggplot(exprSet_L,aes(value,col=group))+geom_density()
print§
p=ggplot(exprSet_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
p=p+stat_summary(fun.y=“mean”,geom=“point”,shape=23,size=3,fill=“red”)
p=p+theme_set(theme_set(theme_bw(base_size=20)))
p=p+theme(text=element_text(face=‘bold’),axis.text.x=element_text(angle=30,hjust=1),axis.title=element_blank())
print§

mean,median,max,min,sd,var,mad

g_mean <- tail(sort(apply(exprSet,1,mean)),50)
g_median <- tail(sort(apply(exprSet,1,median)),50)
g_max <- tail(sort(apply(exprSet,1,max)),50)
g_min <- tail(sort(apply(exprSet,1,min)),50)
g_sd <- tail(sort(apply(exprSet,1,sd)),50)
g_var <- tail(sort(apply(exprSet,1,var)),50)
g_mad <- tail(sort(apply(exprSet,1,mad)),50)
g_mad
names(g_mad)

heatmap

library(pheatmap)
choose_gene=names(tail(sort(apply(exprSet,1,mad)),50))
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)

UpSetR

https://cran.r-project.org/web/packages/UpSetR/README.html
library(UpSetR)
g_all <- unique(c(names(g_mean),names(g_median),names(g_max),names(g_min),
names(g_sd),names(g_var),names(g_mad) ))
dat=data.frame(g_all=g_all,
g_mean=ifelse(g_all %in% names(g_mean) ,1,0),
g_median=ifelse(g_all %in% names(g_median) ,1,0),
g_max=ifelse(g_all %in% names(g_max) ,1,0),
g_min=ifelse(g_all %in% names(g_min) ,1,0),
g_sd=ifelse(g_all %in% names(g_sd) ,1,0),
g_var=ifelse(g_all %in% names(g_var) ,1,0),
g_mad=ifelse(g_all %in% names(g_mad) ,1,0)
)
upset(dat,nsets = 7)
pdata=pData(sCLLex)
group_list=as.character(pdata[,2])
group_list
dim(exprSet)
exprSet[1:5,1:5]

hclust

colnames(exprSet)=paste(group_list,1:22,sep=’’)

Define nodePar

nodePar <- list(lab.cex = 0.6, pch = c(NA, 19),
cex = 0.7, col = “blue”)
hc=hclust(dist(t(exprSet)))
par(mar=c(5,5,5,10))
plot(as.dendrogram(hc), nodePar = nodePar, horiz = TRUE)

PCA

library(ggfortify)
df=as.data.frame(t(exprSet))
df$group=group_list
autoplot(prcomp( df[,1:(ncol(df)-1)] ), data=df,colour = ‘group’)

t.test

dat = exprSet
group_list=as.factor(group_list)
group1 = which(group_list == levels(group_list)[1])
group2 = which(group_list == levels(group_list)[2])
dat1 = dat[, group1]
dat2 = dat[, group2]
dat = cbind(dat1, dat2)
pvals = apply(exprSet, 1, function(x){
t.test(as.numeric(x)~group_list)$p.value
})
p.adj = p.adjust(pvals, method = “BH”)
avg_1 = rowMeans(dat1)
avg_2 = rowMeans(dat2)
log2FC = avg_2-avg_1
DEG_t.test = cbind(avg_1, avg_2, log2FC, pvals, p.adj)
DEG_t.test=DEG_t.test[order(DEG_t.test[,4]),]
DEG_t.test=as.data.frame(DEG_t.test)
head(DEG_t.test)

DEG by limma

suppressMessages(library(limma))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))
rownames(design)=colnames(exprSet)
design
contrast.matrix<-makeContrasts(paste0(unique(group_list),collapse = “-”),levels = design)
contrast.matrix
##这个矩阵声明，我们要把progres.组跟stable进行差异分析比较
##step1
fit <- lmFit(exprSet,design)
##step2
fit2 <- contrasts.fit(fit, contrast.matrix) ##这一步很重要，大家可以自行看看效果
fit2 <- eBayes(fit2) ## default no trend !!!
##eBayes() with trend=TRUE
##step3
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput)
#write.csv(nrDEG2,“limma_notrend.results.csv”,quote = F)
head(nrDEG)

volcano plot

DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )
DEGchange=as.factor(ifelse(DEGchange = as.factor(ifelse(DEGchange=as.factor(ifelse(DEGP.Value < 0.05 & abs(DEGlogFC)>logFCcutoff,ifelse(DEGlogFC) > logFC_cutoff, ifelse(DEGlogFC)>logFCcutoff,ifelse(DEGlogFC > logFC_cutoff ,‘UP’,‘DOWN’),‘NOT’)
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
'\nThe number of up gene is ',nrow(DEG[DEGKaTeX parse error: Expected 'EOF', got '\nThe' at position 21: …e =='UP',]) , '\̲n̲T̲h̲e̲ ̲number of down …change ==‘DOWN’,])
)

g = ggplot(data=DEG, aes(x=logFC, y=-log10(P.Value), color=change)) +
geom_point(alpha=0.4, size=1.75) +
theme_set(theme_set(theme_bw(base_size=20)))+
xlab(“log2 fold change”) + ylab("-log10 p-value") +
ggtitle( this_tile ) + theme(plot.title = element_text(size=15,hjust = 0.5))+
scale_colour_manual(values = c(‘blue’,‘black’,‘red’)) ## corresponding to the levels(res$change)
print(g)

different P values

head(nrDEG)
head(DEG_t.test)
DEG_t.test=DEG_t.test[rownames(nrDEG),]
plot(DEG_t.test[,3],nrDEG[,1])
plot(DEG_t.test[,4],nrDEG[,4])
plot(-log10(DEG_t.test[,4]),-log10(nrDEG[,4]))

exprSet[‘GAPDH’,]
exprSet[‘ACTB’,]
exprSet[‘DLEU1’,]
library(ggplot2)
library(ggpubr)
my_comparisons <- list(
c(“stable”, “progres.”)
)
dat=data.frame(group=group_list,
sampleID= names(exprSet[‘DLEU1’,]),
values= as.numeric(exprSet[‘DLEU1’,]))
ggboxplot(
dat, x = “group”, y = “values”,
color = “group”,
add = “jitter”
)+
stat_compare_means(comparisons = my_comparisons, method = “t.test”)

heatmap

library(pheatmap)
choose_gene=head(rownames(nrDEG),25)
choose_matrix=exprSet[choose_gene,]
choose_matrix=t(scale(t(choose_matrix)))
pheatmap(choose_matrix)

生信技能树

生信人的20个R语言习题的答案相关推荐

生信学习——生信人的20个R语言习题（上）（附详细答案解读）
题目目录 1. 安装一些R包. 2. 了解ExpressionSet对象,比如CLL包里面就有data(sCLLex),找到它包含的元素,提取其表达矩阵(使用exprs函数),查看其大小. 3. 了解 ...
生信学习——生信人的20个R语言习题（下）（附详细答案解读）
题目目录 12. 理解统计学指标mean,median,max,min,sd,var,mad并计算出每个基因在所有样本的这些统计学指标,最后按照mad值排序,取top 50 mad值的基因,得到列表. ...
python perl 比较生信_科学网—生信人写程序1. Perl语言模板及配置 - 刘永鑫的博文...
科学网对Markdown排版支持较差,对格式不满意的用户请跳转至 CSDN 或微信阅读: 如果感觉文章对您有帮助,想继续阅读同类文章,请扫描下方二维码关注"生信宝典"公众号,每天接 ...
生信c语言,生信人的R使用
接下来介绍R语言: [生信技能树]生信人应该这样学R语言 R语言在你开始R之旅前,建议你看看下面这两个 1. 介绍R语言及Rstudio 了解R,Rstudio及R包;安装的包在packages中检 ...
生信人写程序2. Editplus添加Perl, Shell, R模板和语法高亮
https://www.editplus.com/ 前言 "工欲善其事必先利其器",生信工程师每天写代码.搭流程,而且要使用至少三门编程语言,没有个好集成开发环境(IDE,Inte ...
引用另一模板的宏_生信人值得拥有的编程模板Shell
前言 "工欲善其事必先利其器",生信工程师每天写代码.搭流程,而且要使用至少三门编程语言,没有个好集成开发环境(IDE,Integrated Development Environ ...
生信人值得拥有的编程模板-Shell
前言 "工欲善其事必先利其器",生信工程师每天写代码.搭流程,而且要使用至少三门编程语言,没有个好集成开发环境(IDE,Integrated Development Environ ...
生信人的自我修养：Linux 命令速查手册
标题:生信人的自我修养:Linux 命令速查手册目标:致力于为生信人打造一个完整的 Linux 命令速查手册作者:简佐义(jianzuoyi@qq.com) 版本:1.0 日期:2020-11-2 ...
生存分析系列教程（一）使用生信人工具盒进行生存分析
生信人工具盒是生信人团队的开发的一款软件,非常方便.下面我将演示一下如何通过这款软件进行生存分析.为了方便大家理解,形式依然是数据结构-操作-结果解读. 1. 表达矩阵与生存信息矩阵表达矩阵依然 ...

生信人的20个R语言习题的答案

1.安装R包

ggplot2

mean,median,max,min,sd,var,mad

heatmap

UpSetR

hclust

Define nodePar

PCA

t.test

DEG by limma

volcano plot

different P values

heatmap

生信人的20个R语言习题的答案相关推荐

最新文章

热门文章