找不到差异就删样品吗

生信技能树

发布于 2024-12-24 18:38:46

7800

代码可运行

文章被收录于专栏：生信技能树生信技能树

运行总次数：0

代码可运行

转录组授课结束后，就是单细胞实战了。老规矩，仍然是让学员们自己找到感兴趣的的文献和数据集就开始跑我们的的授课代码，其中一个小伙伴提出来了个奇怪的数据集：GSE155179

对应的文章是：《Single-cell transcriptome analysis of human oocyte ageing》. J Cell Mol Med 2021 Jul; 之所以说这个数据集奇怪，因为它确实是单细胞，但是表达量矩阵看起来就跟传统的bulk转录组测序后的表达量矩阵一模一样，看起来非常像一个2分组的传统的bulk转录组测序：

GSM4696891 OOA
GSM4696892 OOD
GSM4696893 OOE
GSM4696894 OOF
GSM4696895 OOG
GSM4696896 OOH

GSM4696897 YOB
GSM4696898 YOC
GSM4696899 YOD
GSM4696900 YOF
GSM4696901 YOG
GSM4696902 YOH

两个不同流派的单细胞转录组技术

Smart-seq2和10x这两个单细胞技术是现在初学者进入单细胞领域最需要掌握的，它们代表着单细胞的两个全然不同的发展策略。

绝大部分的技术原理介绍会从单细胞悬浮液制备到测序细节面面俱到，其实并不那么的初学者友好。给大家推荐了一个高度精炼的综述，这个综述于2020年9月发表在《Experimental & Molecular Medicine》杂志，标题是：《Single-cell sequencing techniques from individual to multiomics analyses》，链接是：https://www.nature.com/articles/s12276-020-00499-2

首先呢，smart-seq2技术依赖于C1这个仪器，每次都是96个细胞一起测序，每个细胞的测序量这个综述可能是写错了，应该是1M-10M为佳，不太可能是100-1000个M，最重要的是它可以覆盖到整个RNA分子的全长测序，每个细胞都是独立的测序，独立的fastq文件哦。所以定量过程基本上等同于传统的bulk转录组测序数据分析，因为细胞数量非常少，细胞数量就等同于样品数量。
然后呢，对于10X技术单细胞转录组呢，每次可以测好几千的细胞，每个细胞只需要5-10K的reads，而且仅仅是测RNA分子的一段即可，全部的细胞都混合在一起是一个fastq文件，虽然说有barcode可以区分，可以拆分成为不同细胞的表达量矩阵。因为每个样品的细胞数量动辄就好几千，甚至几十万上百万的细胞。

流程不能混用

虽然都是单细胞技术，但是上面的数据集：GSE155179给出来的表达量矩阵文件，很明显没办法走我们授课的Seurat的降维聚类分群。应该是两分组的传统的bulk转录组测序后的表达量矩阵的差异分析！

fs=list.files('GSE155179_RAW/', full.names = T)
fs
# [1] "GSE155179_RAW//GSM4696891_OOA_counts.txt.gz"
# [2] "GSE155179_RAW//GSM4696892_OOD_counts.txt.gz"
# [3] "GSE155179_RAW//GSM4696893_OOE_counts.txt.gz"
# [4] "GSE155179_RAW//GSM4696894_OOF_counts.txt.gz"
# [5] "GSE155179_RAW//GSM4696895_OOG_counts.txt.gz"
# [6] "GSE155179_RAW//GSM4696896_OOH_counts.txt.gz"
# [7] "GSE155179_RAW//GSM4696897_YOB_counts.txt.gz"
# [8] "GSE155179_RAW//GSM4696898_YOC_counts.txt.gz"
# [9] "GSE155179_RAW//GSM4696899_YOD_counts.txt.gz"
# [10] "GSE155179_RAW//GSM4696900_YOF_counts.txt.gz"
# [11] "GSE155179_RAW//GSM4696901_YOG_counts.txt.gz"
# [12] "GSE155179_RAW//GSM4696902_YOH_counts.txt.gz"
library(data.table)
tmp= fread(fs[1],data.table = F) 
gid=fread(fs[1],data.table = F)[,1]
head(gid)
length(gid)
length(unique(gid))

rawcount = do.call(cbind,
                   lapply(fs, function(x){
                     fread(x,data.table = F)[,7]
                   }))
rawcount[1:4,1:4] 
dim(rawcount)
kp=!duplicated(gid);table(kp)
rawcount = rawcount[kp,]
gid=gid[kp]
keep <- rowSums(rawcount>0) >= floor(0.75*ncol(rawcount))
table(keep)
rawcount=rawcount[keep,]
dim(rawcount)
gid=gid[keep] 
head(gid) 
rawcount[1:4,1:4]
rownames(rawcount) = gid
if(T){
  gid=gsub('\\.[0-9]*','',gid)
  head(gid) 
  
  library(AnnoProbe)
  ids=annoGene( gid ,
                ID_type = 'ENSEMBL',species = 'human'
  )
  head(ids)
  ids=ids[!duplicated(ids$SYMBOL),]
  
  rawcount=rawcount[match(ids$ENSEMBL , gid),]
  rownames(rawcount) = ids$SYMBOL
}
rawcount[1:4,1:4]
#     [,1] [,2] [,3] [,4]
# AL627309.1  946  251  304  534
# AP006222.1    3    5    2    1
# OR4F29        1    3    2    3
# MTND1P23    220  297  117  354
rawcount['GAPDH',]

接下来就可以基于这个转录组的count矩阵使用金标准算法（DESeq2，edgeR，limma-voom）计算差异基因。首先需要看看质量控制3张图，如下所示：

质量控制3张图

基于这个质量控制3张图，可以推测如果是简简单单差异分析，肯定是符合统计学显著的基因数量会很少很少：

合统计学显著的基因数量会很少很少

但是，很明显文章里面的并不是这样的.文献里面的差异分析首先就不再是我们看得到6vs6啦，而是6vs3的组合，可以看到差异非常明显：

差异非常明显

也就是说，作者进行了样品层面的挑挑拣拣，而且文章给出来了一些基因的名字：

From these DEGs, 17 hub genes were identified including 12 upregulated genes (UBE2C, UBC, CDC34, UBR1, KIF11, ASF1B, PRC1, ESPL1, GTSE1, EXO1, UBA1, KIF4A) and 5 downregulated genes (UBA52, UBE2V2, SKP1, CCNB1, MAD2L1).

在我们的差异分析结果里面，这些基因仍然是可以看到很明显的差异的，有冲突的基因很少很少，就UBA52这个基因比较奇怪：

>   DEG_limma_voom[downregulated_genes,c(1,2,5)]
           logFC   AveExpr adj.P.Val
UBE2C -4.5402492 0.6400363 0.8436287
UBC   -2.7104666 4.3991196 0.8902490
CDC34 -0.9931482 5.6388869 0.8902490
UBR1  -2.0307141 3.9672935 0.8902490
KIF11 -2.2117590 3.2498928 0.8902490
ASF1B -1.4288795 3.4573283 0.9252201
PRC1  -0.9834946 5.7479743 0.8902490
ESPL1 -2.8060765 2.3913733 0.8902490
GTSE1 -1.4277600 6.5825099 0.8782001
EXO1  -4.3751189 1.3859940 0.8463446
UBA1  -0.5402465 3.8810789 0.9773120
KIF4A -0.7650145 8.9052580 0.8902490
>   DEG_limma_voom[upregulated_genes,c(1,2,5)]
              logFC   AveExpr adj.P.Val
UBA52  -0.008606384  6.595008 0.9991876
UBE2V2  1.208057751  5.042271 0.8902490
SKP1    0.536661437  7.856673 0.9361513
CCNB1   0.806690898 11.187017 0.8681044
MAD2L1  0.826284405  9.489264 0.8825524

也就是说，并不需要像作者那样的对样品进行删减，也可以拿到跟做的类似的结果。

那么，问题来了，作者为什么要冒险删除样品呢，做出这样的举动有什么好处呢？让我们看看文章作者是如何描述这个过程的：

作者为什么要冒险删除样品呢

我们也可以在pca上面看看这个样品的情况；

  load(file = 'symbol_matrix.Rdata') 
  # https://mp.weixin.qq.com/s/bGclqc3qbpvm6rvwQtDWTg
  colnames(symbol_matrix)
  symbol_matrix[1:4,1:4] ## 基因名字的样品，矩阵
  dat = log2(edgeR::cpm(symbol_matrix)+1)
  cg=names(tail(sort(apply(dat,1,sd)),1000))#apply按行（'1'是按行取，'2'是按列取）取每一行的方差，从小到大排序，取最大的1000个
  exp=t(dat[cg,])
  library("FactoMineR")#画主成分分析图需要加载这两个包
  library("factoextra")   
  dat.pca <- PCA(exp , graph = FALSE) 
 fviz_pca_ind(dat.pca,
                     geom.ind = c( "point", "text" ), # show points only (nbut not "text")
                     col.ind = group_list, # color by groups
                     palette = "Dark2",repel = T,
                     addEllipses = TRUE, # Concentration ellipses
                     legend.title = "Groups")+ 
    theme(plot.title = element_text(size=12,hjust = 0.5))

从pca图看，有点难啊，如何合理的删除3个样品？

合理的删除3个样品？

如果合理的删除了3个样品然后就重新走一次差异分析流程，但是我做不到，我看了看文章里面的pca图，确实是很容易删除3个样品，虽然不一定合理：

文章里面的pca图

为什么不选择10x呢

因为这个数据集没有走我们默认的10x单细胞转录组，所以我们的单细胞授课代码失效了，但幸运的是前面的表达量芯片和转录组测序授课知识点仍然是有效。然后学员们就七嘴八舌讨论了为什么这个研究不选择10x呢，每个样品动辄就可以好几千甚至一万多细胞。

我问了人工智能大模型，原来是因为每个女性的卵母细胞（oocyte）数量是有限的。女性的卵母细胞数量在出生时就已经确定，并且随着年龄的增长而逐渐减少。女性在出生时拥有大约100万到200万个原始生殖细胞（primordial germ cells），这些细胞在胎儿发育过程中会逐渐减少。到了青春期，一个女性可能只剩下大约30万到50万个卵母细胞。

然而，并非所有的卵母细胞都能成熟并被排卵。在每个月经周期中，通常只有一个卵母细胞会成熟并被排出，尽管在某些情况下可能会有多个。随着时间的推移，女性的生育能力会因为卵母细胞数量的减少而降低，特别是在35岁之后，卵母细胞的质量和数量都会显著下降，这会导致生育能力下降和流产风险增加。也就是说大部分卵子通过细胞凋亡机制而自行退化，一般最后是400多到500个左右卵泡发育成熟并排卵。卵子数量是有限的，女性一生中排卵的数量通常不到五百个，一般是在四百多。仅仅占卵泡总数的百分之零点一左右。

这样的话，细胞数量本来就很有限，就无所谓高通量的10x单细胞技术了。

学徒作业

有一个（Nat Aging 2024）的文章《Spatiotemporal transcriptomic changes of human ovarian aging and the regulatory role of FOXP1》，数据集是GSE202601. ，它引用了这个GSE155179. 样品分组也是类似的：

GSM8077816 Ovary,young-1,scRNAseq
GSM8077817 Ovary,young-3,scRNAseq
GSM8077818 Ovary,young-4,scRNAseq
GSM8077819 Ovary,middle-2,scRNAseq
GSM8077820 Ovary,middle-3,scRNAseq
GSM8077821 Ovary,middle-4,scRNAseq
GSM8077822 Ovary,old-1,scRNAseq
GSM8077823 Ovary,old-2,scRNAseq
GSM8077824 Ovary,old-3,scRNAseq

大家试试看，对GSE202601. 进行降维聚类分群，找到里面的oocyte，看看它的差异分析结果，跟GSE155179. 的是否一致？

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

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

芯片