样品分组都弄错了，数据分析结论还可靠吗

生信技能树

发布于 2024-05-08 10:10:41

840

发布于 2024-05-08 10:10:41

酒香也怕巷子深：好像大家不怎么知道我们生信技能树团队有一个生物信息学入门课，详见；生物信息学马拉松授课（买一得五） ，你值得拥有！

五年前大家做单细胞转录组的时候，三五万左右一个样品的收费让绝大部分课题组都望而却步，一般来说那个时候做少量的几个样品就足够发cns或者子刊级别的文章啦。

但是，很多人珍惜自己的数据，或者说没什么能力把一两个样品的单细胞转录组数据的分析结果编造出来一个生物学故事，卖不出去啊，就只能说是持续贬值了。现在做一个单细胞转录组就五千块钱左右了，时光荏苒岁月如梭啊！现如今做单细胞转录组项目，基本上就跟十年前的普通转录组进化趋势一模一样的，需要有好的实验设计，起码依据自己的生物学背景搞两个分组吧，比如用药前后，转移与否的，生存好坏的不同样品，每个分组起码三五个样品吧。但是这样的话，就有一个很不起眼的风险，就是分组可能会对全局的单细胞转录组表达量没什么特别大的影响，比如取两个分组病人，它们未来是有淋巴结转移这样的结局事件的区别，但是你现在取样他们都是病人的原位肿瘤，这样的话，未来的转移与否并不直接显现在现在的单细胞转录组表达量矩阵的差异上面。

更麻烦的是，有些时候大家有可能自己就把样品分组弄混了，比如2024年1月的文章：《Unveiling inflammatory and prehypertrophic cell populations as key contributors to knee cartilage degeneration in osteoarthritis using multi-omics data integration》，它对应的数据集是：https://www.ncbi.nlm.nih.gov/geo/query/acc.cgi?acc=GSE255460 有如下所示的19个样品：

GSM8072834 C1_Condition_Control_Rep1
GSM8072835 C2_Condition_Control_Rep2
GSM8072836 C3_Condition_Control_Rep3
GSM8072837 OA1-1_Condition_NWB_Rep1
GSM8072838 OA1-2_Condition_WB_Rep1
GSM8072839 OA2-1_Condition_NWB_Rep2
GSM8072840 OA2-2_Condition_WB_Rep2
GSM8072841 OA3-1_Condition_NWB_Rep3
GSM8072842 OA3-2_Condition_WB_Rep3
GSM8072843 OA4-1_Condition_NWB_Rep4
GSM8072844 OA4-2_Condition_WB_Rep4
GSM8072845 OA5-1_Condition_NWB_Rep5
GSM8072846 OA5-2_Condition_WB_Rep5
GSM8072847 OA6-1_Condition_NWB_Rep6
GSM8072848 OA6-2_Condition_WB_Rep6
GSM8072849 OA7-1_Condition_NWB_Rep7
GSM8072850 OA7-2_Condition_WB_Rep7
GSM8072851 OA8-1_Condition_NWB_Rep8
GSM8072852 OA8-2_Condition_WB_Rep8

分组信息也是很明显：19 cartilage samples from 8 OA donors, and 3 non-OA control donors，而且作者给出来了做出来的降维聚类分群结果文件：

GSE255460_RAW.tar 1.0 Gb (http)(custom) TAR (of MTX, TSV)
GSE255460_metadata.csv.gz 1.9 Mb (ftp)(http) CSV
GSE255460_sc_counts.txt.tar.gz 266.3 Mb (ftp)(http) TAR

我们需要把上面的 GSE255460_RAW.tar 1.0 Gb 文件解压后，整理出下面的结构，然后就可以批量读取到Seurat里面做降维聚类分群，很容易构建成为Seurat对象。仍然是走常规的单细胞转录组降维聚类分群代码，可以看链接: https://pan.baidu.com/s/1bIBG9RciAzDhkTKKA7hEfQ?pwd=y4eh ，基本上大家只需要读入表达量矩阵文件到r里面就可以使用Seurat包做全部的流程啦：

整理出下面的结构

很容易对比一下自己的降维聚类分群结果，以及文章给出来了的GSE255460_metadata.csv.gz 1.9 Mb里面的降维聚类分群，如下所示的代码：

  library(stringr)
  phe2=data.table::fread('GSE255460_metadata.csv.gz',
                         data.table = F)
  head(phe2[,c(1,9:12)])
  phe2$barcodes = str_split(phe2$V1,'[.]',simplify = T)[,2]
  length(unique(  phe2$barcodes))
  rownames(phe2)=paste(phe2$ID, phe2$barcodes,sep = '_')
  table(phe2$ID)
   
  phe=sce.all.int@meta.data
  head(phe[,c(1,2)])
  phe$orig.ident =gsub('^GSM[0-9]*_','',phe$orig.ident)
  table(phe$orig.ident) 
  phe$barcodes = str_split(rownames(phe),'_',simplify = T)[,3]
  length(unique(  phe$barcodes))
  rownames(phe)=paste(phe$orig.ident, phe$barcodes,sep = '_')

文章给出来了比较详细的单细胞亚群的生物学名字，我不想浪费时间，就不给了：

>   head(phe2[,c( 9:12)])
                         celltype sample trait group
OA1-1_AAACCCACATCAACCA-1     HomC     S1    OA   NWB
OA1-1_AAACCCACATGTTCAG-1     ProC     S1    OA   NWB
OA1-1_AAACCCAGTCAGGTGA-1     RegC     S1    OA   NWB
OA1-1_AAACCCAGTCTGCAAT-1     HomC     S1    OA   NWB
OA1-1_AAACCCATCAGACTGT-1     ProC     S1    OA   NWB
OA1-1_AAACCCATCTGCGGAC-1     ProC     S1    OA   NWB

>   head(phe[,c(1,2)])
                      orig.ident nCount_RNA
C1_AAACCCAGTATGCGTT-1         C1       2716
C1_AAACCCAGTTTCGTGA-1         C1      12216
C1_AAACCCATCTTCTGGC-1         C1      18764
C1_AAACGAAGTCCTCATC-1         C1      27109
C1_AAACGAAGTTTGACAC-1         C1      17225
C1_AAACGAATCTCGTCAC-1         C1      11918

这个时候我们两个人的细胞表型矩阵里面的行名是可以对应的，因为都是样品名字+细胞的barcode，但是可以看到的是同一个样品在我们两个人的细胞表型矩阵里面的细胞数量居然是差异很大：

>   table(phe2$ID)

   C1    C2    C3 OA1-1 OA1-2 OA2-1 OA2-2 OA3-1 OA3-2 OA4-1 OA4-2 OA5-1 OA5-2 OA6-1 
 1345  8183  6979 11731  8994  5263 10457  5004  5245  9807 12006  6458  5063  6003 
OA6-2 OA7-1 OA7-2 OA8-1 OA8-2 
 8519  5935  4690  6577  7637 
>   table(phe$orig.ident) 

   C1    C2    C3 OA1-1 OA1-2 OA2-1 OA2-2 OA3-1 OA3-2 OA4-1 OA4-2 OA5-1 OA5-2 OA6-1 
 1571 10638  7017 13873 11018  5570 10612  6096  8740  7277  8725  5572  5774 10348 
OA6-2 OA7-1 OA7-2 OA8-1 OA8-2 
12591  6731  5130  6173  5187

这样的问题在于，如果是取交集，就会很可怜：

  ids=intersect(rownames(phe),rownames(phe2))
  table(phe[ids,'orig.ident']) 
  g1=phe[ids,]$RNA_snn_res.1
  g2=phe2[ids,]$celltype
  tmp =  table(g1 ,g2 );tmp
  gplots::balloonplot(
    table(g1 ,g2 )
  )

对应关系

可以看到TMD对应关系并不是很好，这个时候有两个可能性，首先是本来就是单细胞降维聚类分群结果不稳定，亚群直接是没有清晰度界限，分群完全取决于代码运行的参数。另外一个可能性是这交集是错位的。如果是完全错位，理论上它不应该是有一点点比较好的对应关系，既然是有对应成功的，说明是单细胞降维聚类分群结果不稳定。

而且我们两个人的细胞表型矩阵里面的细胞数量都是14万，但是只有5万是可以交集的，看了看，基本上就前面的7个样品是交集成功了：

>   table(phe[ids,'orig.ident']) 

   C1    C2    C3 OA1-1 OA1-2 OA2-1 OA2-2 OA3-1 OA3-2 OA4-1 OA4-2 OA5-1 OA5-2 OA6-1 
 1345  8183  6966 11730  8994  5260 10455     5    14    16    31    11    15    18 
OA6-2 OA7-1 OA7-2 OA8-1 OA8-2 
   25    16     5     8    11

让我想起来了之前遇到了另外一个文章里面的错误，就是不知道10x单细胞转录组样品和fastq文件的对应关系，可以做同样的探索：

  list1_barcodes = split(phe$barcodes,phe$orig.ident)
  list2_barcodes = split(phe2$barcodes,phe2$ID)
  names(list1_barcodes) 
  tmp = do.call(rbind,
          lapply(names(list1_barcodes), function(i){
            unlist(lapply(names(list1_barcodes), function(j){
              #i=names(list1_barcodes)[1];j=i
              sample1_barcodes <- list1_barcodes[[i]]
              sample1_length <- as.numeric(length(sample1_barcodes))
              sample2_barcodes <- list2_barcodes[[j]]
              sample2_length <- as.numeric(length(sample2_barcodes))
              both <- intersect(sample1_barcodes,sample2_barcodes)
              both_barcodes_length <- as.numeric(length(both))
              rate <- as.numeric(2*both_barcodes_length/(sample1_length+sample2_length))
              rate
            }))
          }))
  rownames(tmp)=names(list1_barcodes)
  colnames(tmp)=names(list1_barcodes)
  pheatmap::pheatmap(tmp,
                     show_rownames = T,show_colnames = T,
                     cluster_rows = F,cluster_cols = F)

很明显可以看到确实是前面的7个样品是交集成功了：

前面的7个样品是交集成功了

最后的对应关系如下所示：

Var1 Var2 value

C1 C1 0.922496570644719
C2 C2 0.869560597205249
C3 C3 0.995427264932838
OA1-1 OA1-1 0.916263083893142
OA1-2 OA1-2 0.898860683589846
OA2-1 OA2-1 0.971106803286255
OA2-2 OA2-2 0.992453367506763

OA5-1 OA3-1 0.945915279878971
OA5-2 OA3-2 0.951992013794355
OA6-1 OA4-1 0.972860332423716
OA6-2 OA4-2 0.975891368866122
OA7-1 OA5-1 0.976723026764728
OA7-2 OA5-2 0.992249583047189
OA3-1 OA6-1 0.989503264732623
OA3-2 OA6-2 0.986152152500145
OA8-1 OA7-1 0.979517674264949
OA8-2 OA7-2 0.949478586615369
OA4-1 OA8-1 0.949473076367836
OA4-2 OA8-2 0.9333822271116

这样就很有意思了，作者提供的分组信息也是很明显：19 cartilage samples from 8 OA donors, and 3 non-OA control donors，文章里面的描述如下所示：

the OA samples were split into weight-bearing (WB) cartilage (n=62 611 chondrocytes) from the medial condyle of the femur
non-weight-bearing (NWB) cartilage (n=56 778 chondro- cytes) from the intact lateral condyle of the femur to serve as internal controls

文章里面的取样示意图如下所示：