单样本间的差异分析
分析 差异 样本
2023-06-13 09:16:31 时间
缘起
前两天，曾老师给了我一个8个样本8个组别的转录组数据，即每组只有一个样本的转录组数据。一看到这个数据，还是感到挺震惊的，毕竟作者这样太节省经费了。前面的推文中多次提到了1V1，想了想，此次就用之前的1V1转录组处理流程与原文简单比较下吧。虽然作者没有提供上游处理完的count矩阵，所幸，曾老师那有处理好的count矩阵。接下来，就让我们用count矩阵走自身1V1的流程，与作者的原文结果比较一下吧。
探究
今天，我们使用的转录组数据集源自2019年发表在Diabetes杂志上，文献名称为Sarm1 Gene Defificiency Attenuates Diabetic Peripheral Neuropathy in Mice。该数据集由8个样本组成，每个样本代表一个分组。
转录组数据集介绍
该数据集提交在ENA官网，其PRJ项目号是PRJNA540413。若有小伙伴想从上游开始处理，具体数据下载链接见 https://www.ebi.ac.uk/ena/browser/view/PRJNA540413。该数据集有8个组别，但是8个组别中分为两个部位，SC与SN，处理中含基因敲除鼠与WT鼠，以及STZ处理与药物处理。
正式分析
1.清空环境，加载R包
rm(list = ls())
library(edgeR)
library(DESeq2)
library(FactoMineR)
library(factoextra)
library(clusterProfiler)
library(org.Hs.eg.db)
library(org.Mm.eg.db)
library(stringr)
library(stringi)
library(tidyverse)
library(ggplot2)
library(patchwork)
library(pheatmap)
library(VennDiagram)
library(RColorBrewer)
library(patchwork)
library(ggplotify)
2.读取count数据
rawcount <- read.table("./all_trim.id.txt",header = T,sep="\t",row.names = 1)
colnames(rawcount)
## [1] "Chr" "Start" "End"
## [4] "Strand" "Length" "SRR8988890.sort.bam"
## [7] "SRR8988891.sort.bam" "SRR8988892.sort.bam" "SRR8988893.sort.bam"
## [10] "SRR8988894.sort.bam" "SRR8988895.sort.bam" "SRR8988896.sort.bam"
## [13] "SRR8988897.sort.bam"
rawcount[1:4,1:4]
## Chr
## ENSMUSG00000102693.1 chr1
## ENSMUSG00000064842.1 chr1
## ENSMUSG00000051951.5 chr1;chr1;chr1;chr1;chr1;chr1;chr1
## ENSMUSG00000102851.1 chr1
## Start
## ENSMUSG00000102693.1 3073253
## ENSMUSG00000064842.1 3102016
## ENSMUSG00000051951.5 3205901;3206523;3213439;3213609;3214482;3421702;3670552
## ENSMUSG00000102851.1 3252757
## End
## ENSMUSG00000102693.1 3074322
## ENSMUSG00000064842.1 3102125
## ENSMUSG00000051951.5 3207317;3207317;3215632;3216344;3216968;3421901;3671498
## ENSMUSG00000102851.1 3253236
## Strand
## ENSMUSG00000102693.1 +
## ENSMUSG00000064842.1 +
## ENSMUSG00000051951.5 -;-;-;-;-;-;-
## ENSMUSG00000102851.1 +
# 本身就是基因表达矩阵（无需降重与ID转换）；选择二分组的样本
rawcount=rawcount[,6:13]
keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
filter_count <- rawcount[keep,] #获得filter_count矩阵
express_cpm <- log2(cpm(filter_count)+ 1)
express_cpm[1:4,1:4] #获得cpm矩阵
## SRR8988890.sort.bam SRR8988891.sort.bam
## ENSMUSG00000051951.5 1.5845466 1.6597803
## ENSMUSG00000102331.1 0.5182788 0.6228564
## ENSMUSG00000025900.13 0.2167244 0.2165807
## ENSMUSG00000025902.13 3.7177481 3.7461175
## SRR8988892.sort.bam SRR8988893.sort.bam
## ENSMUSG00000051951.5 2.0250141 1.7654425
## ENSMUSG00000102331.1 0.4685779 0.8318200
## ENSMUSG00000025900.13 0.2196636 0.5655548
## ENSMUSG00000025902.13 3.9206494 4.4710626
#2.获取分组信息
#根据列的信息，提取分组信息
colnames(rawcount)
## [1] "SRR8988890.sort.bam" "SRR8988891.sort.bam" "SRR8988892.sort.bam"
## [4] "SRR8988893.sort.bam" "SRR8988894.sort.bam" "SRR8988895.sort.bam"
## [7] "SRR8988896.sort.bam" "SRR8988897.sort.bam"
group=rep(c("SC","SN"),each=4)
group #直接用部位作为group
## [1] "SC" "SC" "SC" "SC" "SN" "SN" "SN" "SN"
group_list=factor(group,levels = c("SC","SN"))
table(group_list)#检查一下组别数量
## group_list
## SC SN
## 4 4
3.绘制整体表达的箱线图与PCA图
## 01绘制整体表达的箱线图
exprSet <- express_cpm
class(express_cpm)
## [1] "matrix" "array"
## [1] "matrix" "array"
data <- data.frame(expression=c(exprSet),
sample=rep(colnames(exprSet),each=nrow(exprSet)))
p1 <- ggplot(data = data,aes(x=sample,y=expression,fill=sample))+ geom_boxplot() + theme(axis.text.x = element_text(angle = 90)) +
xlab(NULL) + ylab("log2(CPM+1)")+theme_bw()+ theme(axis.text.x = element_text(angle=90, hjust=1, vjust=.5))
p1
## 02绘制PCA图（此处p2直接按照部位画看看吧）
## 从PCA结果中可以看出，相较于干预方式，组织部位间的差异更显著。
dat <- express_cpm
dat <- as.data.frame(t(dat))
dat <- na.omit(dat)
dat$group_list <- group_list
dat_pca <- PCA(dat[,-ncol(dat)], graph = FALSE)#画图仅需数值型数据，去掉最后一列的分组信息
p2 <- fviz_pca_ind(dat_pca,
geom.ind = "point", # 只显示点，不显示文字
col.ind = dat$group_list, # 用不同颜色表示分组
palette = c("#00AFBB", "#E7B800"),
addEllipses = T, # 是否圈起来，少于4个样圈不起来
legend.title = "Groups") + theme_bw()
p1+p2
4.对八组中的两组进行差异分析
此处，就挑选样本号890与891结尾的两个样本SC-WT+Vehicle与SC-WT+STZ组样本进行差异分析吧。
## 取出SC-WT+Vehicle与SC-WT+STZ进行差异分析
filter_count <- filter_count[,c(1,2)] #调整成正常在前；疾病在后
head(filter_count)
## SRR8988890.sort.bam SRR8988891.sort.bam
## ENSMUSG00000051951.5 37 40
## ENSMUSG00000102331.1 8 10
## ENSMUSG00000025900.13 3 3
## ENSMUSG00000025902.13 225 230
## ENSMUSG00000102269.1 2 2
## ENSMUSG00000098104.1 7 5
a=str_split(rownames(filter_count),"\\.",simplify = T)[,1]
rownames(filter_count)=a
rawcount=filter_count
##转换ID
id2symbol <- bitr(rownames(rawcount),
fromType = "ENSEMBL",
toType = "SYMBOL",
OrgDb = org.Mm.eg.db)
rawcount <- cbind(GeneID=rownames(rawcount),rawcount) #多一列换成基因ID
rawcount <- merge(id2symbol,rawcount,
by.x="ENSEMBL",by.y="GeneID",all.y=T)
rawcount=rawcount[!duplicated(rawcount$SYMBOL),]
rawcount <- drop_na(rawcount) #去除含有NA值的行
rownames(rawcount)=rawcount$SYMBOL
rawcount=rawcount[,c(3,4)]
head(rawcount)
## SRR8988890.sort.bam SRR8988891.sort.bam
## Gnai3 732 601
## Cdc45 75 79
## H19 28 59
## Scml2 6 12
## Apoh 16 19
## Narf 996 801
exprSet <- rawcount
dim(exprSet)
## [1] 17964 2
# 查看分组信息和表达矩阵数据
# dim(exprSet)
# exprSet[1:4,1:4]
# group_list
# table(group_list)
# 加载包
library(DESeq2)
exprSet = exprSet[,1:2]
group=rep(c("vehicle","STZ"),each=1)
group
## [1] "vehicle" "STZ"
group_list=factor(group,levels = c("vehicle","STZ"))
table(group_list)#检查一下组别数量
## group_list
## vehicle STZ
## 1 1
group_list=group_list[1:2]
#加载包
library(edgeR)
#设置分组信息
exprSet <- DGEList(counts = exprSet, group = group_list)
bcv = 0.1#设置bcv为0.1
et <- exactTest(exprSet, dispersion=bcv^2)
DEG_edgeR=as.data.frame(topTags(et, n = nrow(exprSet$counts)))
head(DEG_edgeR)
## logFC logCPM PValue FDR
## Hipk2 3.390986 5.608619 2.310790e-43 4.151103e-39
## Myh1 3.174025 4.787663 2.208903e-35 1.984037e-31
## Hoxd10 -3.781282 2.709049 1.295668e-25 7.758462e-22
## Mpz -2.248826 8.213331 1.998350e-25 8.974590e-22
## Capn11 3.546979 2.741779 2.998023e-24 1.077130e-20
## Ogn -2.443415 3.702862 1.632260e-19 4.886986e-16
write.csv(DEG_edgeR, 'single-sample-deg-result.csv', quote = FALSE)
#筛选上下调，设定阈值
fc_cutoff <- 2
pvalue <- 0.01
DEG_edgeR$regulated <- "normal"
loc_up <- intersect(which(DEG_edgeR$logFC>log2(fc_cutoff)),
which(DEG_edgeR$PValue<pvalue))
loc_down <- intersect(which(DEG_edgeR$logFC < (-log2(fc_cutoff))),
which(DEG_edgeR$PValue<pvalue))
DEG_edgeR$regulated[loc_up] <- "up"
DEG_edgeR$regulated[loc_down] <- "down"
table(DEG_edgeR$regulated)
##
## down normal up
## 110 17637 217
library(EnhancedVolcano)
EnhancedVolcano(DEG_edgeR,
lab =rownames(DEG_edgeR),
x = 'logFC',
y = 'FDR')
5.看一下文章描述的几个显著变化基因
cg=c("Pvalb","Cox7a1","Cox6a2")
deg_fc=DEG_edgeR[cg,]
head(deg_fc)
## logFC logCPM PValue FDR regulated
## Pvalb 0.06243797 7.240778 0.7658528260 1.00000000 normal
## Cox7a1 0.02857483 2.302549 0.9565856240 1.00000000 normal
## Cox6a2 1.22077820 2.125091 0.0001561307 0.02018893 up
以上演示了1V1流程对8组中的2组单样本进行差异分析的结果。与原文相比的，我们的上调基因有217个，下调基因有110个，与原文具有一定的区别。但是由于作者没有提供差异基因，我们同样只能看一下作者展示的验证基因吧。
验证的差异基因中Pvalb、Cox7a1与Cox6a2中只有一个发生显著上调，与作者的原文具有一定的区别。这是为什么呢？为什么两者的分析结果存在不同呢？感兴趣的小伙伴们可以点评下。
除此之外，曾老师还提供了一个批量对8次差异分析结果进行差异分析的脚本。由于篇幅与时间问题，我们在下次再对其进行展示吧。
相关文章
- 一阶倒立摆分析_倒立摆受力分析
- 利用perl一键生成符合LEFse差异分析的Table表
- pandas.read_csv() 报错 OSError: Initializing from file failed，报错原因分析和解决方法
- siamfc代码解读_分析SiamFC
- 跟着存档教程动手学RNAseq分析（五）：DESeq2基因水平差异表达分析
- 都是FPKM进行差异分析，为啥差异感觉这么大呢？
- 差异分析分组构建到底谁在前面--关于limma包中model.matrix()的问题
- 【Linux 内核 内存管理】物理分配页 ⑤ ( get_page_from_freelist 快速路径调用函数源码分析 | 遍历备用区域列表 | 启用 cpuset 检查判定 | 判定脏页数量 )
- monocle2轨迹分析
- SQL Server中的分级差异分析（sqlserver中分差）
- 分析Oracle不同版本的特点与差异（oracle 版本区别）
- 范式的差异分析MySQL三范式与反三范式是数据库设计中常用的两种规范，本文将阐述两者间的差异及其使用中的优缺点
- datalist,Repeater和Gridview的区别分析
- java字符串拼接与性能分析详解
- python类继承用法实例分析
- javascript事件模型实例分析