win10下的RNA测序(二)
写在前面
这部分讲的是根据上游得到的counts数做差异分析以及通路富集
读入文件
options(stringsAsFactors = F)
rm(list = ls())
pwd = "E:\\process\\"
setwd(pwd)
library(readr)
library(dplyr)
data = read.csv("exp.txt",header = T,sep = "\t",fill=T,stringsAsFactors=F,skip = 1)
data = data[,c(1,7:15)]# 提取探针以及bam文件对应的counts数
先将探针转换为基因ID
这边转化ID,刚好之前写了个教程https://editor.csdn.net/md/?articleId=115487976,本着用代码处理的原则,则尝试了新方法
library("clusterProfiler")
library(dplyr)
gene = data$Geneid# %>% as.data.frame()
df = data.frame(gene)
#BiocManager::install("AnnotationDbi")
library("AnnotationDbi")
library("org.Mm.eg.db")
df$symbol <- mapIds(org.Mm.eg.db,keys=gene,column="SYMBOL",keytype="ENSEMBL",multiVals="first")
df = na.omit(df)
提取有对应基因的表达谱
如果一个基因对应多个探针的话,取最大值
## 2.1 subset the probe correspondant to gene_symbols
data = subset(data,Geneid%in%df$gene)
rownames(data) = data$Geneid
df = df[match(df$gene,rownames(data)),]
data = data[,-1]
## 2.2 a gene correspondant different probe,you can choose mean or max
data_finally = apply(data,2,function(x){tapply(x,df$symbol,max)
})
差异分析
# 加载包
library(edgeR)
library(limma)# 设定组别
group_list = c(rep("PBS",3),rep("yangxing",3),rep("yinxing",3))
design <- model.matrix(~0+factor(group_list))
colnames(design)=levels(factor(group_list))## 3.1 filter some low express counts and transform to log2
dge <- DGEList(counts=data_finally)
keep <- filterByExpr(dge, design)
dge <- dge[keep,,keep.lib.sizes=FALSE]
dge <- calcNormFactors(dge)
logCPM <- cpm(dge, log=TRUE, prior.count=3)
rownames(design)=colnames(dge)## 3.2 construct contrast matrix and search degs
contrast_class=c("yangxing-PBS","yinxing-PBS","yangxing-yinxing")## 3.3 define function to search degs
search_limma_degs=function(data_subset,group_list,contrast_class){design <- model.matrix(~0+factor(group_list))colnames(design)=levels(factor(group_list))head(design)exprSet=as.matrix(data_subset)mode(exprSet)="numeric"rownames(design)=colnames(exprSet)head(design)# 2.2 contrast.matrixcontrast.matrix<-makeContrasts(contrasts = contrast_class,levels = design)# 3.search the degs # 3.1 degs ##step1,lmfit起过滤的手续fit <- lmFit(exprSet,design)##step2fit2 <- contrasts.fit(fit, contrast.matrix)fit2 <- eBayes(fit2) ## default no trend !!!##eBayes() with trend=TRUEdata_degs=list()for (i in 1:length(contrast_class)){print(i)tempOutput1 = topTable(fit2, coef=i, number =Inf)data_degs[[i]] = na.omit(tempOutput1) }return(data_degs)}
# search degs
data_degs=search_limma_degs(data_subset = logCPM,group_list,contrast_class)
deg_yangxing_PBS = data_degs[[1]]
deg_yinxing_PBS = data_degs[[2]]
deg_yangxing_yinxing = data_degs[[3]]
通路富集与火山图以及相应图表的保存
valcano_pathway = function(nrDEG,i,contrast_class){
df=nrDEG
attach(nrDEG)
df$y= -log10(adj.P.Val)
df$state=ifelse(df$adj.P.Val>0.05,'NO',ifelse( df$logFC >2,'UP',ifelse( df$logFC < -2,'DOWN','NO') )
)
a = table(df$state)
df$state1<-factor(df$state,levels = names(a),labels = paste(names(a),as.numeric(a)))
df$name=rownames(df)
head(df)
library(ggplot2)
library(ggpubr)# valcano
ggplot(df,aes(x=logFC,y=y))+geom_point(aes(color=state1))+labs(x="log2FC",y="-log10FDR")+#增加阈值线:分别对应FDR=0.05,|log2FC|=1geom_hline(yintercept=-log10(0.05),linetype=4)+geom_vline(xintercept=c(-2,2),linetype=4)+theme_classic() + ggtitle(contrast_class[i]) +theme(plot.title = element_text(hjust = 0.5)) +labs(color="Pvalue<0.05\n|log2FC|>2")
ggsave(filename = paste(contrast_class[i],"valcano.png"))
write.csv(df,file = paste(contrast_class[i],".csv"))
# pathway
if (T) {df$ENTREZID <- bitr(df$name, fromType = "SYMBOL",toType = c( "ENTREZID"),OrgDb = org.Mm.eg.db)$ENTREZIDgene_up= df[df$state == 'UP','ENTREZID'] gene_down=df[df$state == 'DOWN','ENTREZID'] ##### up go <- enrichGO(gene_up, OrgDb = "org.Mm.eg.db", ont="all") library(ggplot2)library(stringr)barplot(go, split="ONTOLOGY",font.size =10)+ facet_grid(ONTOLOGY~., scale="free") + scale_x_discrete(labels=function(x) str_wrap(x, width=50))+ggsave(filename = paste(contrast_class[i],".gene_up_GO_all_barplot.png"))enrichKK <- enrichKEGG(gene = gene_up,organism = 'mmu',pvalueCutoff = 0.1,qvalueCutoff =0.1)enrichKK=DOSE::setReadable(enrichKK, OrgDb='org.Mm.eg.db',keyType='ENTREZID')barplot(enrichKK,showCategory=20)ggsave(filename = paste(contrast_class[i],".gene_up_kegg_all_barplot.png"))dotplot(enrichKK)ggsave(filename = paste(contrast_class[i],".gene_up_kegg_all_dotplot.png"))write.csv(go@result,file = paste(contrast_class[i],"gene_up_GO.csv"))write.csv(enrichKK@result,file = paste(contrast_class[i],"gene_up_KEGG.csv"))###### down go <- enrichGO(gene_down, OrgDb = "org.Mm.eg.db", ont="all") library(ggplot2)library(stringr)barplot(go, split="ONTOLOGY",font.size =10)+ facet_grid(ONTOLOGY~., scale="free") + scale_x_discrete(labels=function(x) str_wrap(x, width=50))+ggsave(filename = paste(contrast_class[i],".gene_down_GO_all_barplot.png"))enrichKK <- enrichKEGG(gene = gene_down,organism = 'mmu',pvalueCutoff = 0.1,qvalueCutoff =0.1)enrichKK=DOSE::setReadable(enrichKK, OrgDb='org.Mm.eg.db',keyType='ENTREZID')barplot(enrichKK,showCategory=20)ggsave(filename = paste(contrast_class[i],".gene_down_kegg_all_barplot.png"))dotplot(enrichKK)ggsave(filename = paste(contrast_class[i],".gene_down_kegg_all_dotplot.png"))write.csv(go@result,file = paste(contrast_class[i],"gene_down_GO.csv"))write.csv(enrichKK@result,file = paste(contrast_class[i],"gene_down_KEGG.csv"))}
}
valcano_pathway(nrDEG = deg_yangxing_PBS,1,contrast_class = contrast_class)
valcano_pathway(nrDEG = deg_yinxing_PBS,2,contrast_class = contrast_class)
valcano_pathway(nrDEG = deg_yangxing_yinxing,3,contrast_class = contrast_class)
结果展示
通路富集的另外食用方法
通路富集也可以通过kobas使用
http://kobas.cbi.pku.edu.cn/retrieve/visualization/?taskid=c74679b1b6e04162bf1dccbf43672018&app=gene_list使用
win10下的RNA测序(二)相关推荐
- win10下的RNA测序(一)
参考资料 大多数的下载软件都可以根据这篇博客下载:https://www.jianshu.com/p/6a4855023330 feature_count下载与用法 :https://blog.csd ...
- (备注)win10下wifi和蓝牙二合一模块网络冲突问题
做个备忘录 1.问题来由 这次的问题来自自组笔记本时买的wifi和蓝牙二合一模块. 为了让工作了6.7年的电脑重获新生,萌生自己升级的想法,主板,cpu,内存,固态都更新了一遍,均没有出现问题,表现良 ...
- HC32L110(一) HC32L110 芯片介绍和Win10下DAP-Link, ST-Link, J-Link方式的烧录
目录 HC32L110(一) HC32L110芯片介绍和Win10下的烧录 HC32L110(二) HC32L110在Ubuntu下的烧录 HC32L110(三) HC32L110的GCC工具链和VS ...
- win10下编译OpenCV的微信二维码库给Dotnet使用
文章目录 前言 一.编译OpenCV和opencv_contrib 二.制作Dotnet可以调用的dll 第一步:创建C++空项目并添加一个类 第二步:配置OpenCV环境 第三步:将C++代码编译成 ...
- win10下从源码运行Cura——环境配置(二次开发准备工作)
win10下从源码运行Cura(二次开发准备工作) win10下从源码运行Cura 配置过程(踩坑指南) win10下从源码运行Cura 本人小白一只,近来需要了解Cura,故自己摸爬滚打学习.网上对 ...
- 从引力波探测到RNA测序,AI如何加速科学发现
来源:AI科技评论 本文约6600字,建议阅读5分钟 本文还展示了多科学领域共同面临的挑战和应对策略,希望通过集成和加速的机器学习解决方案为科学发现提供更多示例和灵感. 本篇综述报告主要讨论了机器学习 ...
- 重磅综述:三万字长文读懂单细胞RNA测序分析的最佳实践教程 (原理、代码和评述)
原文链接: https://www.embopress.org/doi/10.15252/msb.20188746 主编评语 这篇文章最好的地方不只在于推荐了工具,提供了一套分析流程,更在于详细介绍了 ...
- 【R语言】Splatter,一个用于简单模拟单细胞RNA测序数据的R包
Splatter是一个用于模拟单细胞RNA测序数据的R包,本文概述并介绍Splatter的功能 一.参数功能 名称 功能 说明 可以通过splatEstimate函数估计 备注 nGenes -> ...
- 单细胞RNA测序研究的实验设计指南(部分阅读)
本篇内容只了解了数据处理与数据分析两个方面,其余方向与计算的关系不是很大,故没有学习.该篇论文的出版时间为2018,我们以流程了解为主,方法新颖性比较小. 目录 数据处理 Normalization ...
最新文章
- python request_python爬虫03 | 那个叫做 Urllib 的库让我们的 python 假装是浏览器
- 基于蔡氏混沌电路进行非线性共振探究
- word文档内容如何防止被复制
- 2014年全国计算机等级一级考试复习资料,2014年全国计算机等级一级考试复习资料..doc...
- BBC News 2012-02-07
- STM32F2系列系统时钟默认配置
- [Flexbox] Using order to rearrange flexbox children
- 批处理延迟sleep应用
- ubuntu下安装ettercap
- 由大脑工作原理,探讨向菩萨求聪明的灵验的科学原理
- 蓝屏代码大全 蓝屏全攻略
- 在Windows Xp上实现Ubuntu主题风格!
- 非计算机专业怎么准备蓝桥杯,大三接触算法,用寒假时间准备蓝桥杯,如何提高成绩?...
- Android View绘制流程
- 如何给计算机硬盘解除密码忘了,电脑硬盘加密忘记密码是怎么处理?
- 地信实验一利用矢量化软件AutoCAD对栅格文件矢量化
- networkx常用函数总结(持续更新)
- 3U VPX接口卡学习资料第288篇:基于FMC接口的Kintex-7 XC7K325T PCIeX4 3U VPX接口卡 数据采集IO卡 软件无线电处理平
- 59.Mongoose
- How to deal with blurred picture