5. 生信技能树——GEO转录组RNA_seq_GSE162550
和生信技能树GEO转录组“GSE150392“分析类似，唯一区别就是在数据处理和ID转换这一环节略微有区别
1.数据下载
最方便的是xena。可以网页下载，也可以用代码下载。
proj = "DHA"
2.生存信息与临床信息
这里仅仅是查看一下，到生存信息部分再整理。
library(GEOquery)
eSet = getGEO("GSE162550",destdir = ".",getGPL = F)
eSet = eSet[[1]]
exp = exprs(eSet)
pd = pData(eSet)
3.表达矩阵行名ID转换
dat = data.table::fread("GSE162550_gene_sample_count_with_symbol (3).xls.gz",data.table = F)k = dat$Symbol!="---";table(k)
dat = dat[k,]k2 = !duplicated(dat$Symbol);table(k2)
dat = dat[k2,]exp = dat[,-(1:3)]
rownames(exp) = dat$Symbol
exp = as.matrix(exp)
4.基因过滤
需要过滤一下那些在很多样本里表达量都为0或者表达量很低的基因。过滤标准不唯一。
过滤之前基因数量：
nrow(exp)
常用过滤标准1：
仅去除在所有样本里表达量都为零的基因
exp1 = exp[rowSums(exp)>0,]
nrow(exp1)
常用过滤标准2(推荐)：
仅保留在一半以上样本里表达的基因
exp = exp[apply(exp, 1, function(x) sum(x > 0) >= 0.5*ncol(exp)), ]
nrow(exp)
5.分组信息获取
根据样本ID的第14-15位，给样本分组（tumor和normal）
Group = rep(c("DMSO","DHA"),each = 3)
Group = factor(Group,levels = c("DMSO","DHA"))
table(Group)
6.保存数据
TCGA以外的数据没有clinical,surv，从下面代码里去掉。
save(exp,Group,proj,file = paste0(proj,".Rdata"))
7.三大R包差异分析
rm(list = ls())
load("DHA.Rdata")
table(Group)
#deseq2----
library(DESeq2)
colData <- data.frame(row.names =colnames(exp), condition=Group)
if(!file.exists(paste0(proj,"_dd.Rdata"))){dds <- DESeqDataSetFromMatrix(countData = exp,colData = colData,design = ~ condition)dds <- DESeq(dds)save(dds,file = paste0(proj,"_dd.Rdata"))
}
load(file = paste0(proj,"_dd.Rdata"))
class(dds)
res <- results(dds, contrast = c("condition",rev(levels(Group))))
#constrast
c("condition",rev(levels(Group)))
class(res)
DEG1 <- as.data.frame(res)
DEG1 <- DEG1[order(DEG1$pvalue),]
DEG1 = na.omit(DEG1)
head(DEG1)#添加change列标记基因上调下调
logFC_t = 2
pvalue_t = 0.05k1 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange < -logFC_t);table(k1)
k2 = (DEG1$pvalue < pvalue_t)&(DEG1$log2FoldChange > logFC_t);table(k2)
DEG1$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG1$change)
head(DEG1)#edgeR----
library(edgeR)dge <- DGEList(counts=exp,group=Group)
dge$samples$lib.size <- colSums(dge$counts)
dge <- calcNormFactors(dge) design <- model.matrix(~Group)dge <- estimateGLMCommonDisp(dge, design)
dge <- estimateGLMTrendedDisp(dge, design)
dge <- estimateGLMTagwiseDisp(dge, design)fit <- glmFit(dge, design)
fit <- glmLRT(fit) DEG2=topTags(fit, n=Inf)
class(DEG2)
DEG2=as.data.frame(DEG2)
head(DEG2)k1 = (DEG2$PValue < pvalue_t)&(DEG2$logFC < -logFC_t);table(k1)
k2 = (DEG2$PValue < pvalue_t)&(DEG2$logFC > logFC_t);table(k2)
DEG2$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))head(DEG2)
table(DEG2$change)
###limma----
library(limma)
dge <- edgeR::DGEList(counts=exp)
dge <- edgeR::calcNormFactors(dge)
design <- model.matrix(~Group)
v <- voom(dge,design, normalize="quantile")design <- model.matrix(~Group)
fit <- lmFit(v, design)
fit= eBayes(fit)DEG3 = topTable(fit, coef=2, n=Inf)
DEG3 = na.omit(DEG3)k1 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC < -logFC_t);table(k1)
k2 = (DEG3$P.Value < pvalue_t)&(DEG3$logFC > logFC_t);table(k2)
DEG3$change = ifelse(k1,"DOWN",ifelse(k2,"UP","NOT"))
table(DEG3$change)
head(DEG3)tj = data.frame(deseq2 = as.integer(table(DEG1$change)),edgeR = as.integer(table(DEG2$change)),limma_voom = as.integer(table(DEG3$change)),row.names = c("down","not","up"));tj
save(DEG1,DEG2,DEG3,Group,tj,file = paste0(proj,"_DEG.Rdata"))
8.可视化
library(ggplot2)
library(tinyarray)
exp[1:4,1:4]
# cpm 去除文库大小的影响
dat = log2(cpm(exp)+1)
pca.plot = draw_pca(dat,Group);pca.plot
save(pca.plot,file = paste0(proj,"_pcaplot.Rdata"))cg1 = rownames(DEG1)[DEG1$change !="NOT"]
cg2 = rownames(DEG2)[DEG2$change !="NOT"]
cg3 = rownames(DEG3)[DEG3$change !="NOT"]h1 = draw_heatmap(dat[cg1,],Group,n_cutoff = 2)
h2 = draw_heatmap(dat[cg2,],Group,n_cutoff = 2)
h3 = draw_heatmap(dat[cg3,],Group,n_cutoff = 2)v1 = draw_volcano(DEG1,pkg = 1,logFC_cutoff = logFC_t)
v2 = draw_volcano(DEG2,pkg = 2,logFC_cutoff = logFC_t)
v3 = draw_volcano(DEG3,pkg = 3,logFC_cutoff = logFC_t)library(patchwork)
(h1 + h2 + h3) / (v1 + v2 + v3) +plot_layout(guides = 'collect') &theme(legend.position = "none")ggsave(paste0(proj,"_heat_vo.png"),width = 15,height = 10)
5. 生信技能树——GEO转录组RNA_seq_GSE162550相关推荐
- 易生信九天的转录组分析培训班总结
易生信九天的转录组分析培训班第一期伴随着5个小时的考试在紧张中结束了.说是培训,倒不如研讨更确切些.在一个个问题的交流中学会转录组分析,效果远大于一人讲,自己练. 先分享两张现场的照片 前两天以集中讲 ...
- CancerSubtypes包的介绍（根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理）
CancerSubtypes包的介绍(根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理,感谢Jimmy老师!) 1. 引言 2. 数据处理 2.1 基本处理 2.1.1 通过检查数据分布来分 ...
- ProTICS包的介绍（根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理）
ProTICS包的介绍(根据生信技能树Jimmy老师分享的乳腺癌分子分型包资料整理,感谢Jimmy老师!) 1.设置环境 2.Part1的结果 3.Part2的结果 4.Part3的结果 5.相关函数 ...
- R：生信技能树学习笔记一
生信技能树小破站:R应该这样学1-4 1.查看已经安装的包的地址 .libPaths() 2.怎么查看函数用法 #在RStudio的右下角窗口的help可以看到 ?函数名 3.三个有用的函数 1.he ...
- R：生信技能树学习笔记二
生信技能树小破站:R应该这样学5-7 1.热图 rm(list=ls()) library(pheatmap) a1=rnorm(100) dim(a1)=c(5,20) #设置维度 pheatmap ...
- 生信技能树【代码大全搜录】
生信技能术代码大全: rm(list = ls()) options()$repos options()$BioC_mirror #options(BioC_mirror="https:// ...
- 生信技能树linux虚拟机,科学网—Windows10安装Linux子系统Ubuntu 20.04LTS，轻松使用生信软件，效率秒杀虚拟机 - 刘永鑫的博文...
很多优秀的生物信息学软件,如QIIME.QIIME 2.LEfSe等没有Windows版,而使用VirutalBox虚拟机不仅效率低,而且挂载外部硬盘和使用中也经常遇到各种问题,配置和使用详见 - 扩 ...
- TCGA学习笔记一(生信技能树概述版)
1.背景介绍 重要数据 外显子数据 表达数据 小RNA测序数据 拷贝数芯片 甲基化数据 蛋白质组学数据 临床信息 癌症背景知识 网页工具大全 GDC cbioportal:按照paper来分类的 UC ...
- 生信技能树课程记录笔记（七）20220531
一.数据框排序 法一:sort函数 默认升序. 例:sort(test$Sepal.Length) 法二:order函数 默认升序,返回数据下标组成的数组. 可以给向量排序,也可以给数据框排序 例:t ...
- 生信技能树课程记录笔记（三）20220526
matrix 矩阵--只允许一种数据 data.frame 数据框--每列只允许一种数据类型 list 列表--各种都可 依据生成函数判断/用class 或i函数判断一下数据类型 一.数据框 新建数据 ...
最新文章
- 怎么解决哈希冲突_从生日悖论谈哈希碰撞
- python画直方图成绩分析-Python数据分析：直方图及子图的绘制
- 9个tcpdump使用实例
- matlab 数字识别_在MATLAB中利用神经网络进行分类
- UPS不间断电源放电时间计算方法
- PHP使用fpdf生成pdf文件（含中文类）
- ASP.NET Core 3.0中支持AI的生物识别安全
- Python机器学习：决策树003使用信息熵寻找最优划分
- mysql timestamp 当前_技术分享 | MySQL 复制那点事 - Seconds_behind_Master 参数调查笔记
- 小程序php上传图片到服务器,关于微信小程序上传图片到服务器的代码
- QListView text动态显示
- 常用的向量矩阵求导公式
- 详细分析推荐系统和搜索引擎的差异陈运文
- 一条语句完成微信、支付宝支付，生成支付二维码
- 融云发送图片消息_融云开发者文档
- 三角形黑盒测试-Java Swing
- 从零开始的RISCV架构CPU设计（2）-CISC与RISC
- 牛客练习赛91A~D
- 崇州 鸡冠山 白塔湖 九龙沟 罨画池 陆游祠{组图及介绍}
- 华为路由器AR2200与部分交换机的密码破解