吨吨吨的生信技能树r语言习题（3）

自从俄罗斯打乌克兰之后真的是让人乌心学习hhhh真的是见证历史了家人们。这个时候感觉不赶快翻墙看视频是对vpn的严重浪费hhhhh
看了一天油土鳖，quora，这个地球实在是太奇妙了，有各种称呼evil animal leaders的，有long live Xi的，有莫斯科游街protest的，各种意识形态各不停输出，简直是魔幻。 其实本来对乌克兰都没什么印象，唯一接触的是学校留学生的漂亮毛妹,但是现在有一种近距离看战争的感触。画面里炮弹连续的声音那么像老家大年三十狂放的鞭炮和烟花，感觉就跟看电影没什么两样。最新的视频里面，乌克兰文化部部长说他们靠的是一种spirit，然后泽连斯基狂打道德牌，很明显已经撑不住了。普皇的一袭万字宣言的震撼还留在脑子里，不仅值得川宝夸赞，要是传销洗脑也都是这个层次，我很难保证不加入。。。。

反正这个地球就是tmd玄幻，正正反反好好坏坏，对于我们这种普通人，可能一辈子也无法明白真相，不要接触政治，要不然真的会变得不幸🙃🙃🙃

可视化专题30题

Q1: 对RNAseq_expr的每一列绘制boxplot图

boxplot(RNAseq_expr)
library(ggplot2)
library(reshape2)
RNAseq_L=melt(RNAseq_expr)
colnames(RNAseq_L)=c('gene','sample','value')
p=ggplot(RNAseq_L,aes(x=sample,y=value))+geom_boxplot()
p

Q2: 对RNAseq_expr的每一列绘制density图

plot(density(RNAseq_expr))
p=ggplot(RNAseq_L,aes(value,col=sample))+geom_density()
p

Q3: 对RNAseq_expr的每一列绘制条形图

barplot(RNAseq_expr)

Q4: 对RNAseq_expr的每一列取log2后重新绘制boxplot图，density图和条形图
有太多0了，去掉全是0的部分,下面就广ggplot了

e=RNAseq_expr[apply(RNAseq_expr,1,function(x) sum(x>0)>1),]
e=log2(e+1)
boxplot(e)
e_L=melt(e)
colnames(e_L)=c('gene','sample','value')
p=ggplot(e_L,aes(x=sample,y=value))+geom_boxplot()
print(p)
p=ggplot(e_L,aes(value,col=sample))+geom_density()
print(p)
ggplot(e_L,aes(x = sample, y = value)) + geom_bar(stat = 'identity')

Q5: 对Q4的3个图里面添加 trt 和 untrt 组颜色区分开来

e_L$group=rep(RNAseq_gl,each=nrow(e))
p=ggplot(e_L,aes(x=sample,y=value,fill=group))+geom_boxplot()
print(p)
p=ggplot(e_L,aes(x=sample,y=value,fill=group))+geom_bar(stat = 'identity')
p
p=ggplot(e_L,aes(value,col=group))+geom_density()+facet_wrap(~group, nrow = 4)
p

Q6: 对RNAseq_expr的前两列画散点图并且计算线性回归方程
我电脑不知道为啥画不出来就只能这样了
y=-0.007 + 0.97x (r)^2 =0.959

 p <- ggplot(e2, aes(x = SRR1039508, y = SRR1039509)) +
+   geom_smooth(method = "lm", se=FALSE, 
+               color="black", formula = y ~ x) +
+   geom_point()
> p
> lm_eqn <- function(x){
+     m <- lm(SRR1039509 ~ SRR1039508, e2);
+     eq <- substitute(italic(y) == a + b %.% italic(x)*","~~italic(r)^2~"="~r2, 
+                      list(a = format(unname(coef(m)[1]), digits = 2),
+                           b = format(unname(coef(m)[2]), digits = 2),
+                           r2 = format(summary(m)$r.squared, digits = 3)))
+     as.character(as.expression(eq));
+ }
 p + geom_text(x = 25, y = 300, label = lm_eqn(e2), parse = TRUE)

先留着
R语言ggplot2散点图添加拟合曲线和回归方程的简单小例子 - 云+社区 - 腾讯云 (tencent.com)
R-ggpmisc|回归曲线添加回归方程，R2，方差表，香不香？ - 知乎 (zhihu.com)

Q7: 对RNAseq_expr的所有列两两之间计算相关系数，并且热图可视化。

M=cor(e)
heatmap(M)

Q8: 取RNAseq_expr第一行表达量绘制折线图

ggplot(e3,aes(x=sample, y=value, group =1)) + geom_line() + geom_point() +
+     theme(axis.text.x = element_text(angle=45, hjust=1, vjust=1,size=8 ))

Q9: 取RNAseq_expr表达量最高的10个基因的行绘制多行折线图

choose_gene=e[names(tail(sort(apply(e,1,sum)),10)),]
choose_gene=melt(choose_gene)
colnames(choose_gene)=c('gene','sample','value')
ggplot(choose_gene,aes(x=sample,y=value,colour=gene,group=gene)) + geom_line() + p=geom_point() +
+     theme(axis.text.x = element_text(angle=45, hjust=1, vjust=1,size=8 ),
+           legend.text = element_text(size = 8, face = "bold"),legend.key.size=unit(0.3,'cm'))

Q10: 一行行的运行 https://github.com/jmzeng1314/5years/blob/master/learn-R/tasks/2-chunjuan-600.R 代码

ggplot绘图

Q10: 一行行的运行：http://biotrainee.com/jmzeng/markdown/ggplot-in-R.html代码
这个就不展示了

生物信息学绘图

需要参考 https://github.com/jmzeng1314/GEO/blob/master/airway_RNAseq/DEG_rnsseq.R

** Q1: 一行行的运行：https://github.com/jmzeng1314/5years/blob/master/learn-R/tasks/top50ggplot.Rmd 代码**
这个是曾老师设置的50个美图的鉴赏

Q2: 对RNAseq_expr挑选MAD值最大的100个基因的表达矩阵绘制热图

choose_gene=e[names(tail(sort(apply(e,1,mad)),100)),]
library(pheatmap)
annotation_col =data.frame(RNAseq_gl)
row.names(annotation_col) <- colnames(choose_gene)
pheatmap(choose_gene,
         color=colorRampPalette(c("navy","white","firebrick3"))(100),
         border_color = NA,
         annotation_col = annotation_col)

Q3: 对RNAseq_expr进行主成分分析并且绘图

library("FactoMineR")#
library("factoextra") 
df=as.data.frame(t(RNAseq_expr))
dat.pca <- PCA(df, graph = FALSE)
fviz_pca_ind(dat.pca,
             geom.ind = "point", # show points only (nbut not "text")
             col.ind = RNAseq_gl, # color by groups
             # palette = c("#00AFBB", "#E7B800"),
             addEllipses = TRUE, # Concentration ellipses
             legend.title = "Groups"
)

Q4: 对RNAseq_expr进行差异分析并且绘制火山图

a=log2(RNAseq_expr+1)
suppressMessages(library(limma)) 
design <- model.matrix(~0+factor(RNAseq_gl))
colnames(design)=levels(factor(RNAseq_gl))
rownames(design)=colnames(a)
design
contrast.matrix<-makeContrasts(paste0(unique(RNAseq_gl),collapse = "-"),levels = design)
fit <- lmFit(a,design)
fit2 <- contrasts.fit(fit, contrast.matrix)
fit2 <- eBayes(fit2)  
tempOutput = topTable(fit2, coef=1, n=Inf)
nrDEG = na.omit(tempOutput) 
head(nrDEG)
DEG=nrDEG
logFC_cutoff <- with(DEG,mean(abs( logFC)) + 2*sd(abs( logFC)) )#设置cutoff
DEG$change = as.factor(ifelse(DEG$P.Value < 0.05 & abs(DEG$logFC) > logFC_cutoff,
                              ifelse(DEG$logFC > logFC_cutoff ,'UP','DOWN'),'NOT')
)
this_tile <- paste0('Cutoff for logFC is ',round(logFC_cutoff,3),
                    '\nThe number of up gene is ',nrow(DEG[DEG$change =='UP',]) ,
                    '\nThe number of down gene is ',nrow(DEG[DEG$change =='DOWN',])
)
this_tile
head(DEG)
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)

Q5: 对RNAseq_expr进行差异分析并且绘制(平均值VS变化倍数)图
刚才用的limma包，现在改用DEseq2包，二者差异：

但是这里用好像结果差还蛮大的。。。。。希望有霸霸解释一下

#<第一步构建对象>
colData <- data.frame(row.names=colnames(exprSet),group_list=group_list) #搞清谁是treat,谁是untreat
dds <- DESeqDataSetFromMatrix(countData = RNAseq_expr,
                              colData = colData,
                              design = ~ group_list)

#可以看到我们构造的dds对象有8个样本，共64102features
suppressMessages(dds2 <- DESeq(dds)) 
resultsNames(dds2)
res <-  results(dds2) 
resOrdered <- res[order(res$padj),]
resOrdered=as.data.frame(resOrdered)
head(resOrdered)
DEG <- as.data.frame(resOrdered)  
table(is.na(DEG))
DEG <-  na.omit(DEG) 
nrDEG <- DEG
DEseq_DEG <- nrDEG
write.csv(DEseq_DEG,file = 'DESeq_DEG.csv',quote = F)

Q6: 绘制其中一个差异基因在两个分组的表达量boxplot并且添加统计学显著性指标

agene=RNAseq_expr[rownames(DEG[100,]),] <我选了第100个
agene=as.data.frame(agene)
agene$gene=rownames(agene)
agene$group=group_list
library(ggpubr)
ggplot(data = agene,aes(y=agene,x=group)) + geom_boxplot() +stat_compare_means(method = "t.test")

Q7: 通过org.Hs.eg.db包拿到RNAseq_expr所有基因的染色体信息，绘制染色体的基因数量条形图

suppressMessages(library(org.Hs.eg.db))
CHR <- toTable(org.Hs.egCHR)
ensembl <- toTable(org.Hs.egENSEMBL)
ensembl_id <- data.frame(rownames(RNAseq_expr))
colnames(ensembl_id) <- c("ensembl_id")
eg <- merge(ensembl_id,ensembl,by="ensembl_id")
egc <- merge(eg,CHR,by="gene_id")
colnames(egc)

Q8: 在上面染色体的基因数量条形图并列叠加差异基因数量条形图

DEG_ensembl <- data.frame(rownames(nrDEG))
colnames(DEG_ensembl) <- c("ensembl_id")
b=merge(DEG_ensembl,egc,by='ensembl_id')
egc$DEG <- as.factor(ifelse(egc$ensembl_id %in% b$ensembl_id,'DEG','NOT'))
ggplot(egc,aes(x=chromosome,fill=DEG)) + 
+     geom_bar()+
+     theme_classic2()

Q9: 在oncolnc网页工具拿到CUL5基因在BRCA数据集的表达量及病人生存资料自行本地绘制生存分析图

ggbetweenstats(data = CUL5, x='Status', y='Expression')
suppressPackageStartupMessages(library(survminer))
suppressPackageStartupMessages(library(survival))
suppressPackageStartupMessages(library(ggstatsplot))
table(CUL5$Status)  #查看多少活得多少死的
CUL5$Status=ifelse(CUL5$Status=='Dead',1,0) #死了是1活着是0
sfit <- survfit(Surv(Days, Status)~Group, data=CUL5)#画生存曲线
sfit
summary(sfit)
ggsurvplot(sfit, conf.int=F, pval=TRUE)

Q10: 在xena网页工具拿到CUL5基因在BRCA数据集的表达量及病人的PAM50分类并且绘制分类的boxplot

CULP5<- read.table('denseDataOnlyDownload.tsv',sep = '\t',header = T,na.strings = '')
View(CULP5)
library(RColorBrewer)
ggplot(data = CULP5,aes(y=CUL5,x=PAM50_mRNA_nature2012,color=PAM50_mRNA_nature2012)) + geom_boxplot() +
+     theme_light() + stat_boxplot(geom = "errorbar",width=0.2)