10X文库送测序后,从测序公司拿到的测序数据是fastq格式的,要经过linux上跑cellranger程序,得到表达矩阵,才能做后面的功能分析。这里就是讲一下如何跑cellranger。
biomamba没有做视频,但是提供了文字版学习资料:单细胞分析的最上游——处理Fastq文件:cellranger。链接如下:
https://mp.weixin.qq.com/s?__biz=MzAwMzIzOTk5OQ==&mid=2247484923&idx=1&sn=b5876af14fbee68d1e9db4b0f70cd1c8&chksm=9b3f7cabac48f5bdb3aac7d20201d89121a83720b41a10017564fedbb8c4a13df0da92daffd7&scene=21#wechat_redirect
一、下载安装cellranger
https://support.10xgenomics.com/single-cell-gene-expression/software/downloads/latest
在表格中填写具体信息。之后会转到下载界面。可以直接点击:Download - Linux 64-bit – 600 MB下载,再上传到服务器。也可以在服务器上直接用curl或者wget命令下载。点击右下角红色数字,也可以下载老版本的cellranger。我选择下载稍微老一些的版本文件名是cellranger-6.1.2.tar。要把这个文件放到software文件夹下,而不是直接就在家目录下,因为最好要将软件归类放置。
在服务器上进入到software文件夹下,输入:tar -xzvf cellranger-6.1.2.tar.gz 来解压缩。会生成一个新的文件夹名为cellranger-6.1.2。
二、在北鲲云服务器里设置环境变量,添加系统路径
首先通过WinSCP软件登录北鲲云服务器,用WinSCP比直接从网页界面操作要更方便,因为.bashrc等隐藏文件,在软件里可以点击直接打开,如同文本一样操作保存就行。可是网页界面就不能直接打开编辑隐藏文件。
在.bashrc文件里面加一句:export PATH=/home/cloudam/software/cellranger-6.1.2:$PATH。之后到家目录下运行:source .bashrc,来激活环境。
在任何目录下,输入:cellranger,都会显示如下界面,表明安装cellranger成功。
三、准备reference
一般常用的就是人和小鼠的,在前面下载cellranger的界面下面就有下载链接。
下载文件:refdata-gex-GRCh38-2020-A.tar.gz和refdata-gex-mm10-2020-A.tar.gz,放到reference文件夹下。
在服务器进入到reference文件夹下解压就行:tar -xzvf refdata-gex-GRCh38-2020-A.tar.gz ,tar -xzvf refdata-gex-mm10-2020-A.tar.gz。会生成相应种属的文件夹,其中包括参考基因组序列、gtf文件以及star的索引文件等。
四、准备fastq格式的测序文件
公司测序回来的文件都有一定的格式,Illumina测序仪下机FASTQ命名类似Ery_S1_L004_R2_001.fastq.gz(V2有12个文件)
比如下图
1,最前面Ery是样本名,和填写测序单上样本名一样。
2,S后跟的数字与样本在sampleSheet中的顺序一致,从1开始。同一个样本名可能对应多个S,因为会有一个样本有多个index的情况,在10X单细胞V2版本的实验里,一个样本就对应4个index,也就会依次从S1-S4。
3,L后面表明在那个line上。
4,同一个样本同一个index有三个fastq文件,I1为index,R1时barcode和UMI。R2才是测序read,文件最大。
五、cellranger对fastq文件进行分析
cellranger有好几个命令,最核心的是cellranger count。
cellranger count这个命令是把fastq文件中的序列比对到参考转录组上并产生一个以.cloupe为结尾的文件以便在loupe cell browser上分析,同时会产生多个与目前主流分析软件兼容的文件以便进一步分析。
cellranger count
--id=run_count_1kpbmcs \
--fastqs=/mnt/home/user.name/yard/run_cellranger_count/pbmc_1k_v3_fastqs \
--sample=pbmc_1k_v3 \
--transcriptome=/mnt/home/user.name/yard/run_cellranger_count/refdata-cellranger-GRCh38-3.0.0
--nosecondary
说明
--id是自己起的,将来会生成这个id名的文件夹,分析结果统统在里面。
--fastqs是fastq数据的具体保存路径。
--sample是S1前面那个样本名。
--transcriptome是reference的路径
--nosecondary是
六、建立slurm脚本,运行
建立slurm脚本:
#!/bin/bash
#SBATCH --output=cellranger.out
#SBATCH --error=cellranger.err
#SBATCH --mail-type=end
#SBATCH --mail-user=zmeraner@126.com
project=~/singlecell #项目文件夹
cellranger count --id=Li1_cellranger --fastqs=$project/fastq --sample=Li1 --transcriptome=/home/cloudam/reference/refdata-gex-mm10-2020-A
因为10X官网显示,cellranger的运行条件为:
- 8-core Intel or AMD processor (16 cores recommended)
- 64GB RAM (128GB recommended)
- 1TB free disk space
- 64-bit CentOS/RedHat 7.0 or Ubuntu 14.04; See the [10x Genomics OS Support]
Note: Cell Ranger v6.1 was the last version that supported CentOS/RedHat 6 or Ubuntu 12.04
北鲲云服务器好像硬盘只有200G,有点儿少。我的数据有36G。不知道可以运行不?
用sinfo
命令查看可以选择的队列。CPU分区命名规则为c-核心数-每核心内存大小,如c-8-4:表示单节点规格为8核,每核心有4G内存,即节点规格为8核32G。
输入:
sbatch -p c-16-4 cellranger.slurm
提交作业输入:
squeue
查看作业运行情况JOBID:作业号。ST:状态 (R:运行中;CF:配置中;PD:排队中)。
提交之后用squeue查看,先显示为CF,几分钟后显示为PD,又过了两分钟左右显示为R。
运行了大约21.5个小时之后,结果出来了。
Li1_cellranger文件夹下有个out文件夹。下载这个out文件夹到本地电脑,里面的文件很多,有些文件可以进行质控,有些是可以用cloupe打开看分群特征的,有些是可以用做其他分析的输入文件。
参考教程:https://zhuanlan.zhihu.com/p/390516422?ivk_sa=1024320u
1,web_summary.html这个文件打开就可以看到这个样本的质控信息。
2,cloupe文件可用cloupe软件打开。
3,filtered_feature_bc_matrix文件夹中是过滤过后的表达矩阵,可以对接searat等后期分析
4,raw_feature_bc_matrix文件夹是未过滤的,没啥用。
5,possorted_genome_bam.bam是比对的bam文件,包含每个reads的,所以这个文件挺大。
七、cellranger count结果解读之summary
左上角Sequencing模块
Number of Reads总reads pair的量
Valid Barcodes是有效的barcode——表示包了beads合格的液滴比率。Valid UMIs是有效的UMI——表示合格的RNA序列比率。实验没问题的话,这两个参数都能达到大于95%以上。
Q30的几个参数,一般情况都应该大于85%,说明测序质量不错。
右上角cells模块
Estimated Number of Cells是预估的细胞数。根据曲线图,上机每个样本都有几万个barcode,但是真正包入细胞的,才能有足够多的UMI为蓝色部分,其他都是灰色部分没有啥UMI的液滴。
Fraction Reads in Cells表示在确定为cell的barcode中的reads占到总reads的比率,不低于80%才好。低于70%的话,认为实验有问题或者数据质量不好。
Mean Reads per Cell表示每个细胞平均reads,一般只要在20-30K reads/cell应该就够了。
Median Genes per Cell为每个细胞检测到的基因数量的中位数,大于1000更好,有利于后面分群。如果小于500认为可能不太可靠。
Total Genes Detected鉴定到的基因总数,这个参数没啥太大意义,和物种细胞类型相关,有些细胞表达的基因种类就是比较少。有些组织复杂度高,细胞种类丰富,那基因种类也就多。常见的在1.2-2万之间。
Sequencing saturation是测序饱和度,只要在80%以上就完全可以了,再高就浪费测序量了,60%-80%就可以了。
左下角mapping模块
为比对到各个不同位置上的比率。包括全基因组上,基因间区,外显子,内含子,转录本区等。
Reads Mapped Confidently to Genome这个比对率一般都能到85%以上。
Reads Mapped Confidently to Exonic Regions应该在60%以上。
右下角sample模块
样本信息,包括名称,试剂版本,比对使用的reference,cellranger版本等。
八、cellranger count结果解读之analysisi
1,t-SNE Projection 分群情况
用 t-SNE算法分群的两个图。每个点儿代表一个细胞。左图为每个细胞中含的UMI数量。右图为分群图。
2,Top Features by Cluster (Log2 fold-change, p-value)
按照上图分群后,不同群之间的差异表达基因列表。
可以在这个列表中类似excel操作排序,查看各个群的基因表达情况。比如cluster1列中L2FC值越高的基因,表明cluster1里这些基因比所有其他群的表达都要多,那么应该可以从L2FC高的基因里面找找有没有这一群的marker。。
3,Sequencing Saturation饱和度评估图
和summary里面的Sequencing saturation参数相对应。
4,Median Genes per Cell
和summary里面Median Genes per Cell参数也是相对应的。
上两个图都是对reads抽样,观察不同抽样条件下检测到的转录本数量占检测到的所有转录本的比例,并绘制曲线。发现抽样越多,饱和度越高,每个细胞基因数的中位数也越高。只要饱和度大于80%都是很不错了。因为基本就可以代表整个样本了。
下图我这个曾经做的实验,测的就过多了,饱和度都到98%了。当时预计细胞6000多个细胞,但实际上只捕获到600来个细胞,所以就测多了,呜呜。