单样本间的差异分析

生信菜鸟团

发布于 2023-01-05 21:10:52

2K00

代码可运行

文章被收录于专栏：生信菜鸟团生信菜鸟团

运行总次数：0

代码可运行

缘起

前两天，曾老师给了我一个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次差异分析结果进行差异分析的脚本。由于篇幅与时间问题，我们在下次再对其进行展示吧。

本文参与腾讯云自媒体同步曝光计划，分享自微信公众号。

原始发表：2022-12-28，如有侵权请联系 cloudcommunity@tencent.com 删除

count

本文分享自生信菜鸟团微信公众号，前往查看

如有侵权，请联系 cloudcommunity@tencent.com 删除。

本文参与腾讯云自媒体同步曝光计划，欢迎热爱写作的你一起参与！

count

登录后参与评论

0 条评论

热度