单细胞测序分析上游cellranger软件入门介绍

公共数据库中的SRA 单细胞转录组数据究竟包含了哪些数据?CellRanger怎样利用10x平台下机数据进行下游一系列分析?这篇文章简单记录CellRanger 包括的主要分析步骤,纯理论。

SRA 原始数据转fastq

公共数据库的SRA 数据需要借助fastq-dump 转为fastq文件,然后进行质控、CellRanger定量等操作。相较于普通转录组数据,原始SRA数据会得到3个fastq文件,分别是Library 的Index(8bp)文件,包括长度为26bp 的Barcode(16bp)和UMI(10bp)的Read1文件,和测序reads文件。

conda install -c bioconda sra-tools ## 安装软件
wkd=/home/project/single-cell/MCC
cd $wkd/raw/P2586-4
cat SRR_Acc_List-2586-4.txt |while read i
do
time fastq-dump --gzip --split-files -A $i ${i}.sra && echo "** ${i}.sra to fastq done **"
done
### 单细胞数据参数为 --split-files 而不是 --split-3

image
image

i7 sample index (library barcode)
是加到Illumina测序接头上的,保证多个测序文库可以在同一个flow-cell上或者同一个lane上进行混合测序(multiplexed)。它的作用就是在CellRanger的mkfastq 功能中体现出来的,它自动识别样本index名称(例如:SA-GA-A1),将具有相同4种oligo的fq文件组合在一起,表示同一个样本。它保证了一个测序lane上可以容纳多个样本

Barcode
是10X特有的,用来区分GEMs,也就是对细胞做了一个标记。一般在拆分混样测序数据(demultiplexing)这个过程后进行操作,当然这也很符合原文的操作。

UMI
UMI就是Unique Molecular Identifier,由4-10个随机核苷酸组成,在mRNA反转录后,进入到文库中,每一个mRNA随机连上一个UMI,根据PCR结果可以计数不同的UMI,最终统计mRNA的数量。它的主要作用是,处理PCR 扩增偏差,因为起始文库很小时需要的PCR扩增次数就越多,因为越容易引入扩增误差。

fastq文件更名
为什么要更名?CellRanger 定量过程输入文件指定命名格式。
如何更名?下图格式:

image
# 比如,将原来的SRR7692286_1.fastq.gz改成SRR7692286_S1_L001_I1_001.fastq.gz
# 依次类推,将原来_2的改成R1,将_3改成R2
cat  SRR_Acc_List-9245-3.txt | while read i ;do (mv ${i}_1*.gz ${i}_S1_L001_I1_001.fastq.gz;mv ${i}_2*.gz ${i}_S1_L001_R1_001.fastq.gz;mv ${i}_3*.gz ${i}_S1_L001_R2_001.fastq.gz);done

fastq 文件质控
主要查看Q30比例。

# 以P2586-4为例
mkdir -p $wkd/qc
cd $wkd/qc
find $wkd/raw/P2586-4 -name '*R1*.gz'>P2586-4-id-1.txt
find $wkd/raw/P2586-4 -name '*R2*.gz'>P2586-4-id-2.txt
cat P2586-4-id-1.txt P2586-4-id-2.txt >P2586-4-id-all.txt

cat P2586-4-id-all.txt| xargs fastqc -t 20 -o ./

CellRanger 流程

它主要包括四个主要基因表达分析流程:

  • cellranger mkfastq : 它借鉴了Illumina的bcl2fastq ,可以将一个或多个lane中的混样测序样本按照index标签生成样本对应的fastq文件。即对下机数据base calling files转为fastq文件。
  • cellranger count :利用mkfastq生成的fq文件,进行比对(基于STAR)、过滤、UMI计数。利用细胞的barcode生成gene-barcode矩阵,然后进行样本分群、基因表达分析
  • cellranger aggr :接受cellranger count的输出数据,将同一组的不同测序样本的表达矩阵整合在一起,比如tumor组原来有4个样本,PBMC组有两个样本,现在可以使用aggr生成最后的tumor和PBMC两个矩阵,并且进行标准化去掉测序深度的影响
  • cellranger reanalyze :接受cellranger countcellranger aggr生成的gene-barcode矩阵,使用不同的参数进行降维、聚类。属于定制化分析。

它的结果主要是包含有细胞信息的BAM, MEX, CSV, HDF5 and HTML文件

CellRanger 软件安装及测试实战操作参考:CellRanger走起(四) CellRanger流程概览。包括安装详细指南及其他一些小技巧。

拆分数据 mkfastq、细胞定量 count、定量组合 aggr

mkfastq 拆分
目的:将每个flowcell 的Illumina sequencer's base call files (BCLs)转为fastq文件

特色: 它借鉴了Illumina出品的bcl2fastq,另外增加了:

  • 将10X 样本index名称与四种寡核苷酸对应起来,比如A1孔是样本SI-GA-A1,然后对应的寡核苷酸是GGTTTACT, CTAAACGG, TCGGCGTC, and AACCGTAA ,那么程序就会去index文件中将存在这四种寡核苷酸的fastq组合到A1这个样本
  • 提供质控结果,包括barcode 质量、总体测序质量如Q30、R1和R2的Q30碱基占比、测序reads数等
  • 可以使用10X简化版的样本信息表
### 两种使用方式
# 第一种
$ cellranger mkfastq --id=bcl \
                     --run=/path/to/bcl \
                     --samplesheet=samplesheet-1.2.0.csv
# 第二种
$ cellranger mkfastq --id=bcl \
                     --run=/path/to/bcl \
                     --csv=simple-1.2.0.csv
# 其中id指定输出目录的名称,run指的是下机的原始BCL文件目录
# 重要的就是测序lane、样本名称、index等信息

参考文章CellRanger走起(四) CellRanger流程概览 中解释了几种不同方式输出fastq文件的不同存放目录,可根据实际分析自行查找。

count 定量
这个过程是最重要的,它完成细胞与基因的定量,它将比对、质控、定量都包装了起来,内部流程很多,但使用很简单。每个版本要求的参数是不同的,尤其是V2与V3版本存在较大差异,这里先对V2进行了解。

# 这是示例,不是真实数据 #
cellranger count --id=sample345 \
                   --transcriptome=/opt/refdata-cellranger-GRCh38-1.2.0 \
                   --fastqs=/home/scRNA/runs/HAWT7ADXX/outs/fastq_path \
                   --sample=mysample \
                   --expect-cells=1000 \
                   --nosecondary
# id指定输出文件存放目录名
# transcriptome指定与CellRanger兼容的参考基因组
# fastqs指定mkfastq或者自定义的测序文件
# sample要和fastq文件的前缀中的sample保持一致,作为软件识别的标志
# expect-cells指定复现的细胞数量,这个要和实验设计结合起来
# nosecondary 只获得表达矩阵,不进行后续的降维、聚类和可视化分析(因为后期会自行用R包去做)

###输出文件
Outputs:
- Run summary HTML:                      /opt/sample345/outs/web_summary.html
- Run summary CSV:                       /opt/sample345/outs/metrics_summary.csv
- BAM:                                   /opt/sample345/outs/possorted_genome_bam.bam
- BAM index:                             /opt/sample345/outs/possorted_genome_bam.bam.bai
- Filtered gene-barcode matrices MEX:    /opt/sample345/outs/filtered_gene_bc_matrices
- Filtered gene-barcode matrices HDF5:   /opt/sample345/outs/filtered_gene_bc_matrices_h5.h5
- Unfiltered gene-barcode matrices MEX:  /opt/sample345/outs/raw_gene_bc_matrices
- Unfiltered gene-barcode matrices HDF5: /opt/sample345/outs/raw_gene_bc_matrices_h5.h5
- Secondary analysis output CSV:         /opt/sample345/outs/analysis
- Per-molecule read information:         /opt/sample345/outs/molecule_info.h5
- Loupe Cell Browser file:               /opt/sample345/outs/cloupe.cloupe

Pipestance completed successfully!

输出文件从上到下依次来看:

  • web_summary.html:官方说明 summary HTML file
  • metrics_summary.csv:CSV格式数据摘要
  • possorted_genome_bam.bam:比对文件
  • possorted_genome_bam.bam.bai:索引文件
  • filtered_gene_bc_matrices:是重要的一个目录,下面又包含了 barcodes.tsv.gz、features.tsv.gz、matrix.mtx.gz,是下游Seurat、Scater、Monocle等分析的输入文件
  • filtered_feature_bc_matrix.h5:过滤掉的barcode信息HDF5 format
  • raw_feature_bc_matrix:原始barcode信息
  • raw_feature_bc_matrix.h5:原始barcode信息HDF5 format
  • analysis:数据分析目录,下面又包含聚类clustering(有graph-based & k-means)、差异分析diffexp、主成分线性降维分析pca、非线性降维tsne
  • molecule_info.h5:下面进行aggregate使用的文件
  • cloupe.cloupe:官方可视化工具Loupe Cell Browser 输入文件

count 定量的重点和难点在于参考基因组索引的构建,可详细参考CellRanger走起(四) CellRanger流程概览,文章中还介绍了怎样自行构建参考基因组索引,及一比对过程的重要细节信息。

aggr
当处理多个生物学样本或者一个样本存在多个重复/文库时,最好的操作就是先分别对每个文库进行单独的count定量,然后将定量结果利用aggr组合起来.

  1. 得到count结果
  2. 构建Aggregation CSV
  3. 运行aggr
### step1
$ cellranger count --id=LV123 ...
... wait for pipeline to finish ...
$ cellranger count --id=LB456 ...
... wait for pipeline to finish ...
$ cellranger count --id=LP789 ...
... wait for pipeline to finish ...

## step2
# AGG123_libraries.csv
library_id,molecule_h5
LV123,/opt/runs/LV123/outs/molecule_info.h5
LB456,/opt/runs/LB456/outs/molecule_info.h5
LP789,/opt/runs/LP789/outs/molecule_info.h5
# 其中
# molecule_h5:文件molecule_info.h5 file的路径 

### step3
cellranger aggr --id=AGG123 \
                  --csv=AGG123_libraries.csv \
                  --normalize=mapped
# 结果输出到AGG123这个目录中

cellranger count的结果

count结果一般放在out目录下,主要有summary和analysis两大类

Summary
是一个HTML格式文件,包括实验捕获的细胞数目、检测到的基因数目、测序数据的产出与质量统计、参考基因组的比对情况。

image

几个指标可以关注一下:

  • 左上部分中,包括了reads数、barcodes数、UMI、index、Q30等统计值

  • 左下是reads比对的比例,包括基因间区、外显子、内含子,如果比对率太低(一般认为外显子的比对率要在60%以上)

  • 右上图是利用barcodes上的UMI标签分布来估计细胞数,绿/蓝色表示细胞,灰色表示背景,其中Y轴表示每个barcode对应的UMI数量,X轴是一定数量的UMI序列所对应barcode的数量,比如上图中有1000个barcodes含有10k个UMI,细胞过滤就是通过这个图来展示的。

  • 首先明确,barcode用来区分细胞,UMI用来区分转录本;其次,barcodes数量时要大于细胞数量的(以保证每个细胞都会有barcode来进行区分)如果图线陡降说明

另参考文章CellRanger走起(二) CellRanger out 结果 中还记录了一个手动构建物种GTF文件的案例尝试,及相应的检验方法,可以根据需要自行学习。

总之,无论是从公共数据库下载单细胞SRA数据,还是来自10x 平台的原始下机数据,或者先经过fastq-dump转为fastq文件,然后质控,然后CellRanger count 定量,或者首先自行拆分数据mkfastq,然后构建索引进行细胞定量,及可选的定量组合等。CellRanger的总体流程见本文记录,但是其中的各个细节仍需要自行实践总结。

作者:新欣enjoy
链接:https://www.jianshu.com/p/5157ab9f6977

©著作权归作者所有,转载或内容合作请联系作者
  • 序言:七十年代末,一起剥皮案震惊了整个滨河市,随后出现的几起案子,更是在滨河造成了极大的恐慌,老刑警刘岩,带你破解...
    沈念sama阅读 151,511评论 1 330
  • 序言:滨河连续发生了三起死亡事件,死亡现场离奇诡异,居然都是意外死亡,警方通过查阅死者的电脑和手机,发现死者居然都...
    沈念sama阅读 64,495评论 1 273
  • 文/潘晓璐 我一进店门,熙熙楼的掌柜王于贵愁眉苦脸地迎上来,“玉大人,你说我怎么就摊上这事。” “怎么了?”我有些...
    开封第一讲书人阅读 101,595评论 0 225
  • 文/不坏的土叔 我叫张陵,是天一观的道长。 经常有香客问我,道长,这世上最难降的妖魔是什么? 我笑而不...
    开封第一讲书人阅读 42,558评论 0 190
  • 正文 为了忘掉前任,我火速办了婚礼,结果婚礼上,老公的妹妹穿的比我还像新娘。我一直安慰自己,他们只是感情好,可当我...
    茶点故事阅读 50,715评论 3 270
  • 文/花漫 我一把揭开白布。 她就那样静静地躺着,像睡着了一般。 火红的嫁衣衬着肌肤如雪。 梳的纹丝不乱的头发上,一...
    开封第一讲书人阅读 39,672评论 1 192
  • 那天,我揣着相机与录音,去河边找鬼。 笑死,一个胖子当着我的面吹牛,可吹牛的内容都是我干的。 我是一名探鬼主播,决...
    沈念sama阅读 31,112评论 2 291
  • 文/苍兰香墨 我猛地睁开眼,长吁一口气:“原来是场噩梦啊……” “哼!你这毒妇竟也来了?” 一声冷哼从身侧响起,我...
    开封第一讲书人阅读 29,837评论 0 181
  • 序言:老挝万荣一对情侣失踪,失踪者是张志新(化名)和其女友刘颖,没想到半个月后,有当地人在树林里发现了一具尸体,经...
    沈念sama阅读 33,417评论 0 228
  • 正文 独居荒郊野岭守林人离奇死亡,尸身上长有42处带血的脓包…… 初始之章·张勋 以下内容为张勋视角 年9月15日...
    茶点故事阅读 29,928评论 2 232
  • 正文 我和宋清朗相恋三年,在试婚纱的时候发现自己被绿了。 大学时的朋友给我发了我未婚夫和他白月光在一起吃饭的照片。...
    茶点故事阅读 31,316评论 1 242
  • 序言:一个原本活蹦乱跳的男人离奇死亡,死状恐怖,灵堂内的尸体忽然破棺而出,到底是诈尸还是另有隐情,我是刑警宁泽,带...
    沈念sama阅读 27,773评论 2 234
  • 正文 年R本政府宣布,位于F岛的核电站,受9级特大地震影响,放射性物质发生泄漏。R本人自食恶果不足惜,却给世界环境...
    茶点故事阅读 32,253评论 3 220
  • 文/蒙蒙 一、第九天 我趴在偏房一处隐蔽的房顶上张望。 院中可真热闹,春花似锦、人声如沸。这庄子的主人今日做“春日...
    开封第一讲书人阅读 25,827评论 0 8
  • 文/苍兰香墨 我抬头看了看天上的太阳。三九已至,却和暖如春,着一层夹袄步出监牢的瞬间,已是汗流浃背。 一阵脚步声响...
    开封第一讲书人阅读 26,440评论 0 180
  • 我被黑心中介骗来泰国打工, 没想到刚下飞机就差点儿被人妖公主榨干…… 1. 我叫王不留,地道东北人。 一个月前我还...
    沈念sama阅读 34,523评论 2 249
  • 正文 我出身青楼,却偏偏与公主长得像,于是被迫代替她去往敌国和亲。 传闻我的和亲对象是个残疾皇子,可洞房花烛夜当晚...
    茶点故事阅读 34,583评论 2 249

推荐阅读更多精彩内容