单细胞多组学数据分析不会分析

生信技能树jimmy

发布于 2023-08-31 03:17:00

34200

代码可运行

文章被收录于专栏：单细胞天地单细胞天地

运行总次数：0

代码可运行

不开玩笑，只需要你设计实验收集样品花钱在公司测序了即可，数据分享到公开数据库后就可以发表在《 BMC Genomic Data》杂志啦，比如2023年8月7号见刊的《The single-cell transcriptome and chromatin accessibility datasets of peripheral blood mononuclear cells in Chinese holstein cattle》，实验设计蛮简单的：

单细胞多组学技术：single-cell RNA sequencing (scRNA-seq) and single-cell sequencing assay for transposase-accessible chromatin (scATAC-seq).
四个分组： four whole-blood treatments (no, 2 h, 4 h, and 8 h LPS)
得到的细胞数量：7,107 (no), 9,174 (2 h), 6,741 (4 h), and 3,119 (8 h) cells

因为脂多糖(LPS)组成革兰氏阴性细菌的外膜，其暴露可导致牛的局部或全身炎症水平升高，所以本实验设计就是使用脂多糖(LPS)看LPS 作为慢性炎症的关键介质调节免疫应答。

全文没有一个图，但是数据是实打实的公开了，GSE225962，如下所示的样品：

GSM7061075 no LPS, scRNA-seq
GSM7061076 2 h LPS, scRNA-seq
GSM7061077 4 h LPS, scRNA-seq
GSM7061078 8 h LPS, scRNA-seq
GSM7061079 no LPS, scATAC-seq
GSM7061080 2 h LPS, scATAC-seq
GSM7061081 4 h LPS, scATAC-seq
GSM7061082 8 h LPS, scATAC-seq

虽然说给的文件有点奇怪：

GSM7061075_C.scRNAexpression.txt.gz 15.5 Mb
GSM7061076_T1.scRNAexpression.txt.gz 17.8 Mb
GSM7061077_T2.scRNAexpression.txt.gz 15.0 Mb
GSM7061078_T3.scRNAexpression.txt.gz 8.0 Mb
GSM7061079_C.scATACgenes-barcode.matrix.txt.gz 1.2 Mb
GSM7061080_T1.scATACgenes-barcode.matrix.txt.gz 10.5 Mb
GSM7061081_T2.scATACgenes-barcode.matrix.txt.gz 8.7 Mb
GSM7061082_T3.scATACgenes-barcode.matrix.txt.gz 5.4 Mb

单细胞转录组矩阵很容易读取并且降维聚类分群，然后看看随着脂多糖(LPS)处理时间段变化的基因，通路以及细胞亚群，但是单细胞ATAC数据作者给出来的文件应该是不够的，可能是需要去 PRJNA938112 里面下载原始数据后进行处理啦。感兴趣的可以试试看：

scATAC-seq4: scATAC-seq上游分析
scATAC-seq3:常用工具—SnapATAC简介
scATAC-seq1：由转录组到表观组
scATAC-seq2: scATAC-seq技术原理

单细胞ATAC-seq

同样的，单细胞ATAC-seq也是上下游独立开，走在Linux系统的cellranger-atac软件进行上游分析，然后走R语言里面的下游统计可视化即可。

这里我们拿刚刚发表的范文举例：西湖大学和浙江大学的科研团队合作的，发表于2022年6月的，发表在Cell Discovery 的文章：《Single-cell multiomics analysis reveals regulatory programs in clear cell renal cell carcinoma》，链接是：https://www.nature.com/articles/s41421-022-00415-0

是关于 clear cell renal cell carcinoma (ccRCC) 的肿瘤微环境的单细胞多组学，包括：

single-cell RNA sequencing (scRNA-seq)
single-cell assay for transposase-accessible chromatin sequencing (scATAC-seq)

数据在PRJNA768891，需要自己下载其测序数据，如下所示：

$ ls -lh *gz|cut -d" " -f5-
28G 8月   2 15:14 SRR16213608_S1_L001_R1_001.fastq.gz
28G 8月   2 15:18 SRR16213608_S1_L001_R2_001.fastq.gz
29G 8月   2 15:23 SRR16213609_S1_L001_R1_001.fastq.gz
29G 8月   2 15:27 SRR16213609_S1_L001_R2_001.fastq.gz
42G 8月   2 15:33 SRR16213610_S1_L001_R1_001.fastq.gz
42G 8月   2 15:40 SRR16213610_S1_L001_R2_001.fastq.gz
41G 8月   2 15:46 SRR16213611_S1_L001_R1_001.fastq.gz
34G 8月   2 15:51 SRR16213611_S1_L001_R2_001.fastq.gz
39G 8月   2 15:57 SRR16213612_S1_L001_R1_001.fastq.gz
34G 8月   2 16:02 SRR16213612_S1_L001_R2_001.fastq.gz
44G 8月   2 16:09 SRR16213613_S1_L001_R1_001.fastq.gz
37G 8月   2 16:15 SRR16213613_S1_L001_R2_001.fastq.gz
37G 8月   2 16:21 SRR16213614_S1_L001_R1_001.fastq.gz
33G 8月   2 16:26 SRR16213614_S1_L001_R2_001.fastq.gz

简单的看 PRJNA768891描述信息，可以知道前面的3个是 (scATAC-seq) ，后面的4个是(scRNA-seq) ，需要分开独立走cellranger流程，所以需要自行去cellranger官网注册并且下载软件和数据库文件，我们一直在做单细胞，所以下面的数据库文件是不同时间段下载的不同版本：

   34M 3月   1 17:06 aspera-connect-3.7.4.147727-linux-64.tar.gz
  972M 3月   1 17:03 cellranger-4.0.0.tar.gz
  955M 3月   1 17:03 cellranger-5.0.1.tar.gz
  1.1G 3月   1 17:03 cellranger-6.0.0.tar.gz
  528M 3月  24 10:50 cellranger-atac-2.1.0.tar.gz

 14G 5月   3  2021 refdata-cellranger-arc-GRCh38-2020-A-2.0.0.tar.gz
   11G 3月   1 17:04 refdata-gex-GRCh38-2020-A.tar.gz
  9.9G 3月   1 17:05 refdata-gex-GRCh38-and-mm10-2020-A.tar.gz
  9.7G 3月   1 17:06 refdata-gex-mm10-2020-A.tar.gz
  
  
  3.3M 3月   1 17:03 refdata-cellranger-vdj-GRCh38-alts-ensembl-4.0.0.tar.gz
3.3M 3月  24 15:20 refdata-cellranger-vdj-GRCh38-alts-ensembl-5.0.0.tar.gz
3.5M 3月  24 15:20 refdata-cellranger-vdj-GRCm38-alts-ensembl-5.0.0.tar.gz

大家作为初学者，可以一次性下载最新版即可。这个步骤需要的计算资源比较大，所以我们一次性跑完结果给大家哈。

交流群

因为这个Cell Discovery 的文章：《Single-cell multiomics analysis reveals regulatory programs in clear cell renal cell carcinoma》，非常贴心的整理了其全套单细胞多组学下游分析R语言代码给大家，在：GitHub (https://github.com/Dragonlongzhilin/RenalTumor).

我下载并且解压看了看，还是有很多可取之处，所以组建交流群号召大家一起解读一下这些代码，而且我们提供这个文章附带的PRJNA768891数据集的上游分析结果给大家哈。

就是你不做单细胞多组学，打开这套代码也可以观摩一下里面的自定义函数技巧：

    41 Combined.P.FC.R
     105 Dot.plot.R
      12 Filter.gene.R
      43 IDConvert.R
     125 Integrate.scRNA.scATAC.R
      91 Plot_colorPaletters.R
     294 analysis.diff.survival.TCGA.R
      40 doubletDetect.R
      31 plot.violin.R
      29 ratio.plot.R
      18 variableFeatureSelection.R

因为单细胞转录组我们已经有了大量教程（4个系列两百多个笔记），大家很容易得到文章里面的降维聚类分群和生物学命名，如下所示：

降维聚类分群和生物学命名

可以看到主要是5个淋巴系免疫亚群：

CD4+ (CD4, IL7R, CD3D, CD3E)
CD8+ T cells (CD8A, CD8B, CD3D, CD3E)
Treg (FOXP3, IL2RA)
natural killer (NK)/natural killer T (NKT) cells (KLRD1, GNLY)
B cells (MS4A1/CD20, CD79A)

以及5个髓系细胞

macrophages (CSF1R, CD68, CD163),
monocytes (S100A12, FCGR3A/CD16),
mast cells (TPSAB1, KIT)
Dendritic cell
Neutrophil

其实就是我们一直确定的肿瘤单细胞降维聚类分群规则啦，第一层次分群也是肿瘤上皮细胞，基质细胞（内皮，成纤维等），免疫细胞（淋巴系的T, NK 和B，髓系的巨噬和单核）。

所以本次交流群我们会重点解读单细胞ATAC-seq的代码哦，如下所示：

     208 1.mergeData.R
      79 2.cluster.R
     118 3.1.integrate.scATAC.scRNA3000.R
     189 3.2.AssignCellType.R
     541 4.1.callPeak&DAR.R
     213 4.2.peakAnnotation.genomicLocation.R
      70 4.3.peak.annotated.geneHancer.R
     138 5.1.motif.enrichment.R
     297 5.2.motif.analysis.R 
      39 5.3.2.drug.plot.R
     124 6.1.cis-coassessibility.R
     204 6.2.ccans.annotated.genomicLocation.R
     122 7.1.All.TF.target.R
     321 7.2.Tumor..TFs.regulatoryNetwork.R
     480 7.Tumor.TFs.regulatoryNetwork.R
     552 8.1.Immune.CD8T.R
     536 8.1.Immune.CD8T.TF.target.R
     507 8.2.Immune.Macrophage.R
     525 8.2.Immune.Macrophage.TF.target.R
      86 Endothelial.R
     814 ExtraDataAnalysis.R
      76 NK.analysis.R