只能下载bam文件的10x单细胞转录组项目数据处理

通常情况下，我们希望一个文献在发表后可以公布它的10x单细胞转录组项目数据 ，并且是以3个表达量矩阵文件（barcodes.tsv.gz，matrix.mtx.gz，genes.tsv.gz或者features.tsv.gz）的格式给我们。因为下游处理的时候，一定要保证这3个文件同时存在，而且在同一个文件夹下面。这样可以直接读取文件夹的路径即可，走seurat流程进行单细胞降维聚类分群。这样的基础分析，也可以看基础10讲：

但是有一些文章并不会给我们表达量矩阵，而是原始数据，比如sra文件或者fastq文件，就需要自己走cellranger的定量流程啦。我们在单细胞天地多次分享过cellranger流程的笔记，大家可以自行前往学习，如下：

因为这个流程其实是需要10X单细胞转录组的fastq文件，而且呢，需要你的fastq文件的命名是有规则的！

但是，更麻烦的是有一些文章给的居然是bam文件，就是它其实自己走了cellranger流程但是并不给我们表达量矩阵文件，而是cellranger流程里面的star软件比对后的bam文件。比如文章：《Determinants of Response and Intrinsic Resistance to PD-1 Blockade in Microsatellite Instability–High Gastric Cance》，我们很容易通读全文，拿到其数据文件地址，项目编号是PRJEB40416，并且成功下载，如下所示：

$ ls -lh *bam|cut -d" " -f 5-
22G 5月  23 02:31 EP-72-B.possorted_genome_bam.bam
32G 5月  23 02:17 EP-72-F.possorted_genome_bam.bam
35G 5月  23 02:00 EP-75-F.possorted_genome_bam.bam
23G 5月  23 01:35 EP-76-B.possorted_genome_bam.bam
23G 5月  23 01:22 EP-76-F.possorted_genome_bam.bam
21G 5月  23 01:09 EP-77-B.possorted_genome_bam.bam
22G 5月  23 00:21 EP-77-F.possorted_genome_bam.bam
23G 5月  22 23:32 EP-78-B.possorted_genome_bam.bam
30G 5月  22 22:39 EP-78-F.possorted_genome_bam.bam

这个时候，就需要 走 cellranger的 bamtofastq流程，把下载的 bam转fastq文件，代码如下所示：

bin=/home/bakdata/x10/pipeline/cellranger-6.1.2/bin/cellranger
# 单个样品
$bin bamtofastq --nthreads 3  --traceback raw/EP-72-B.possorted_genome_bam.bam EP-72-B

每个样品都会输出如下所示的多个fq文件，它是有规律的哦 ：

$ ls -lh *gz|cut -d" " -f 5-
730M 5月  23 10:30 bamtofastq_S1_L001_I1_001.fastq.gz
738M 5月  23 10:35 bamtofastq_S1_L001_I1_002.fastq.gz
733M 5月  23 10:40 bamtofastq_S1_L001_I1_003.fastq.gz
737M 5月  23 10:45 bamtofastq_S1_L001_I1_004.fastq.gz
735M 5月  23 10:50 bamtofastq_S1_L001_I1_005.fastq.gz
257M 5月  23 10:52 bamtofastq_S1_L001_I1_006.fastq.gz
1.1G 5月  23 10:30 bamtofastq_S1_L001_R1_001.fastq.gz
1.1G 5月  23 10:35 bamtofastq_S1_L001_R1_002.fastq.gz
1.1G 5月  23 10:40 bamtofastq_S1_L001_R1_003.fastq.gz
1.1G 5月  23 10:45 bamtofastq_S1_L001_R1_004.fastq.gz
1.1G 5月  23 10:50 bamtofastq_S1_L001_R1_005.fastq.gz
397M 5月  23 10:52 bamtofastq_S1_L001_R1_006.fastq.gz
2.7G 5月  23 10:30 bamtofastq_S1_L001_R2_001.fastq.gz
2.7G 5月  23 10:35 bamtofastq_S1_L001_R2_002.fastq.gz
2.7G 5月  23 10:40 bamtofastq_S1_L001_R2_003.fastq.gz
2.6G 5月  23 10:45 bamtofastq_S1_L001_R2_004.fastq.gz
2.6G 5月  23 10:50 bamtofastq_S1_L001_R2_005.fastq.gz
1.1G 5月  23 10:52 bamtofastq_S1_L001_R2_006.fastq.gz

这个项目是5个病人的9个样品，每个样品都是几十个fq文件，而文件名字规律很明显，如果以下划线为分隔符，那么简单的介绍一下：

• 第2列是S1到S99这些选择，通常情况下，保证都是S1就足够了。
• 第3列是L001到L999这些选择，上面的例子里面都是L001
• 第5列是R1，R1，I1这样的3种情况，这个时候每个样品的I1是可以省略的，但是R1和R2必须存在并且完整。
• 第6列是001到999的选项，比如这个吗里面就是001到006的6种情况。

然后就正常走cellranger的定量流程即可，代码我已经是多次分享了。参考：

差不多几个小时就可以完成全部的样品的cellranger的定量流程，我简单检查了一下这5个病人的9个单细胞转录组样品走cellranger的定量流程，发现他们的质量并不好，比如：

质量并不好

结论就是这个样品它绝大部分的reads虽然是比对成功，而且也确实是在外显子区域，如下所示：

Reads Mapped to Genome 98.6%
Reads Mapped Confidently to Genome 97.1%
Reads Mapped Confidently to Intergenic Regions 8.1%
Reads Mapped Confidently to Intronic Regions 7.8%
Reads Mapped Confidently to Exonic Regions 81.1%
Reads Mapped Confidently to Transcriptome 79.3%
Reads Mapped Antisense to Gene 0.7%

但是它并不在最后确定的450个有效细胞里面， 绝大部分的比对成功的reads都落到了那些被cellranger的定量流程不认为是有效的细胞里面，相当于是白白浪费了哦。

这样的质量差有两个可能性，要么是high levels of ambient RNA ，就是说很多比对成功的reads都是环境污染，这个需要细致的处理，但是没办法解决细胞数量太少的问题。另外一个可能性是本次建库得到的单细胞普遍RNA数量都很少，所以被cellranger的定量流程不认为是有效的细胞，这个可以调整参数，加上--force-cells会显著改善。

拿这5个病人的9个单细胞转录组样品走cellranger的定量流程走降维聚类分群结果如下所示：

降维聚类分群

可以看到，仍然是能比较清晰的看到不同亚群的差异，文章里面是：

文章的降维聚类分群

很难去对比了，因为文章对其单细胞转录组数据描述太少了，没有阈值信息，没有软件参数，甚至没有细胞数量的直接描述。

我想，这就是为什么文章并不想给表达量矩阵给大家的原因吧，毕竟绝大部分小伙伴没办法做到下载原始数据并且处理它。

文末友情宣传

强烈建议你推荐给身边的博士后以及年轻生物学PI，多一点数据认知，让他们的科研上一个台阶：

• 数据挖掘（GEO,TCGA,单细胞）2022年6月场，快速了解一些生物信息学应用图表
• 生信入门课-2022年6月场，你的生物信息学第一课