DropSeq 简介
Drop-seq技术最初由麻省理工学院的McCarroll实验室开发,于2015年首次发表在《Cell》杂志上。自此以后,这种技术已被广泛应用于许多研究领域,包括神经科学、发育生物学和肿瘤学等。
Drop-seq是一种单细胞RNA测序技术,通过在微滴中捕获单个细胞并进行RNA扩增,从而获得单个细胞的转录组数据。该技术通过微滴分离单个细胞并将细胞裂解,随后在微滴中添加反转录酶和一种称为“barcode beads”的特殊珠子,这些珠子上有一个独特的序列标识符。这些珠子能够被反转录酶转录RNA,并在RNA末端添加一段序列,从而将每个RNA分子与其来源细胞进行标记。然后,通过PCR扩增这些标记的RNA分子,将其构建成文库并进行高通量测序,以获得单个细胞的转录组数据。
DropSeq技术与其他单细胞RNA测序技术相比,具有高通量、高灵敏度和高特异性等优点,能够在细胞数量有限的情况下实现单细胞RNA测序,并广泛应用于单细胞转录组学研究。
DropSeq 文库结构:
read 1中包含12bp的barcode,8bp的UMI,read2中是捕获的RNA分子。
figure-2D.pngDrop-seq数据前期分析
软件安装
下面列出的软件如果你已经安装了,则跳过此步骤。
安装DropSeq Tools
DropSeq 数据通常使用DropSeq Tools做基础分析,可以进行 reads 到 barcodes 和 UMIs 的mapping、去除低质量 reads、去除 PCR 重复、生成 gene expression 矩阵等操作。
DropSeq Tools github下载最新版地址:broadinstitute/Drop-seq
wget https://github.com/broadinstitute/Drop-seq/releases/download/v2.5.2/Drop-seq_tools-2.5.1.zip
unzip Drop-seq_tools-2.5.1.zip
安装 bgzip
bgzip
是一个常用的将文本压缩成 BGZF 格式的命令行工具,通常用于对大规模的基因组数据进行压缩,比如 VCF 文件。
以下是安装 bgzip
的步骤:
- 安装 zlib 库:
bgzip
依赖于 zlib 库,需要先安装该库。在 Ubuntu 系统下可以使用以下命令安装:
sudo apt-get update
sudo apt-get install zlib1g-dev
在其他 Linux 发行版中,可以使用相应的包管理器安装 zlib 库。
- 下载安装 htslib 库:
bgzip
是 htslib 库的一部分,需要先安装 htslib 库。可以从 htslib 的官方网站(https://github.com/samtools/htslib/releases)下载最新版本的源代码,解压后进入该目录,使用以下命令进行编译和安装:
make
sudo make install
如果编译过程中出现错误,可能需要安装一些其他的库和工具,如 OpenSSL、pkg-config 等。可以根据提示安装相应的库和工具。
- 安装 bgzip:安装完 htslib 库后,即可使用以下命令安装
bgzip
:
git clone https://github.com/samtools/tabix.git
cd tabix && make
sudo make install
安装完成后,可以使用 bgzip
命令进行文件压缩和解压缩。
安装STAR
DropSeq使用STAR软件作为比对软件,所以需要下载安装STAR
在latest release页面找到最新版的STAR,download即可
cd /path/to/Software/star/
wget https://github.com/alexdobin/STAR/releases/download/2.7.10b/STAR_2.7.10b.zip
unzip STAR_2.7.10b.zip
准备2个文件
在分析流程之前,需要准备好参考基因组和转录本注释文件,这些文件可以从多个资源中获取,例如:
-
Ensembl(https://www.ensembl.org/index.html)和UCSC(https://genome.ucsc.edu/)等数据库提供了大量的参考基因组和注释文件,可以根据需要下载。
-
NCBI提供的 RefSeq 数据库(https://www.ncbi.nlm.nih.gov/refseq/)也提供了参考基因组和注释文件,可以在 NCBI 的 FTP 网站(ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/)上找到并下载。
-
如果需要使用人类基因组的参考文件,可以考虑使用 GATK 团队提供的资源包(https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle),其中包含了多个版本的参考基因组和注释文件。
-
10X Genomics 也提供了已经处理好的基因组和注释文件。在其官网可以获得。
一般来说,参考基因组文件是一个 fasta 格式的文件,包含了染色体序列和对应的基因组坐标信息。转录本注释文件可以是 gtf 或 gff3 格式的文件,包含了基因和转录本的注释信息,以及它们在基因组上的位置信息。
因为这里涉及到一些meta信息,我这里从Ensembl下载,下载sm_primary_assemble的faata,或者 primary_assemble.fa。从头构建一套用于dropseq 分析的文件。
image-20230312011633163.pngmkdir reference # 随便创建个文件夹放一下
wget https://ftp.ensembl.org/pub/release-109/fasta/homo_sapiens/dna/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz
wget https://ftp.ensembl.org/pub/release-109/gtf/homo_sapiens/Homo_sapiens.GRCh38.109.gtf.gz
下载完成后需要解压,
gzip -d Homo_sapiens.GRCh38.109.gtf.gz
gzip -d Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa.gz
在下载参考基因组和注释文件后,可以将它们作为参数传递给 create_Drop-seq_reference_metadata.sh
脚本,并根据脚本提示的要求提供相应的文件路径即可。
生成一些 meta 文件
在 Drop-seq 流程中,需要使用参考基因组和转录本注释文件(gtf)对测序数据进行比对和定量。这个脚本的作用就是帮助用户生成这些必要的metadata文件。
create_Drop-seq_reference_metadata.sh \
-n Homo_sapiens_genome_annotation \
-r /reference/Homo_sapiens.GRCh38.dna_sm.primary_assembly.fa \
-s GRCh38 \
-g /reference/Homo_sapiens.GRCh38.109.gtf \
-d /software/biosoftware/Drop-seq_tools-2.5.1 \
-o . \
-t . \
-a /usr/local/bin/STAR \
-b /usr/local/bin/bgzip \
-i /usr/local/bin/samtools
# 注意:STAR,bgzip,samtools 如果在环境变量中,就可以省略这三个参数。
create_Drop-seq_reference_metadata.sh
脚本通常会为参考基因组创建 .dict
文件,该文件包含参考基因组的字典信息,是 Drop-seq 流程所需的重要参考文件之一。会重新生成一下fasta,gtf,一些intervals, 还会构建STAR index。
在 create_Drop-seq_reference_metadata.sh
脚本中,picard CreateSequenceDictionary
命令用于创建 .dict
文件。该命令会读取参考基因组的 .fasta
文件,并基于参考基因组序列的名称、长度、描述等信息生成 .dict
文件。需要注意的是,这里使用了 Picard 工具中的 CreateSequenceDictionary
命令来创建 .dict
文件,而不是使用 STAR 工具创建。这是因为 Picard 工具在处理字典信息时更加严格和精确,可以确保 Drop-seq 流程的正确性。
这一步生成的这些文件要用在后面alignment中作为参数,而不用刚刚下载的原始reference。
生成表达矩阵
第一步:fqtobam
fastq 文件必须要转换为Picard-queryname-sorted BAM 文件,也就是用picard,按照queryname 排序后的bam文件。
DropSeq tools 文件夹中带了picard.jar,那就不用自己安装了。路径在:Drop-seq_tools-2.5.1/3rdParty/picard/picard.jar
使用:
java -jar Drop-seq_tools-2.5.1/3rdParty/picard/picard.jar FastqToSam \
F1=file1.R1.fq.gz \
F2=file1.R2.fq.gz \
O=unaligned_read_pairs.bam \
SM=My_Sample
# 上述写法虽然可以运行,但是picard更新了语法,
# 我们与时俱进,用下面的新语法更时尚:
java -jar picard.jar FastqToSam \
--FASTQ file1.R1.fq.gz \
--FASTQ2 file1.F2.fq.gz \
--OUTPUT unaligned_read_pairs.bam \
--SAMPLE_NAME My_Sample
第二步:alignment
Drop-seq_alignment.sh
是一个all-in-one的脚本,它将 Drop-seq 实验中的多个处理步骤整合到一个脚本中,简化了数据处理的流程。
具体来说,Drop-seq_alignment.sh
包含了以下几个步骤:
-
TagBamWithReadSequenceExtended
:将每个read的UMI和barcode信息添加到bam文件的read名中。 -
FilterBAM
:用于过滤低质量的reads和reads mapping到基因组外的区域。 -
TrimStartingSequence
:用于去除 read 的起始序列。 -
PolyATrimmer
:用于去除read的3'端的polyA序列。 -
TagReadWithInterval
:TagReadWithInterval函数将每个比对结果中的Read映射到参考基因组上的转录本进行分类,并将分类结果添加到比对结果的SAM格式文件中。 -
TagReadWithGeneFunction
:根据所属基因的注释信息(如基因名、基因类型、外显子位置等),将注释信息添加到比对结果的SAM格式文件中的XF标签中。 -
DetectBeadSubstitutionErrors
:检测Drop-seq实验中可能存在的珠子替换错误(Bead Substitution Error),以保证单细胞RNA测序的准确性。 -
DetectBeadSynthesisErrors
:检测每个UMI的错配率,并将潜在的bead合成错误的UMI标记为bad_umi。
用法:
/Drop-seq_tools-2.5.1/Drop-seq_alignment.sh \
-g /path/to/STAR \
-r /path/to/Homo_sapiens_genome_annotation.fasta.gz \
-d /software/biosoftware/Drop-seq_tools-2.5.1 \
-o . \
-t . \
-s /usr/local/bin/STAR \
unaligned_read_pairs.bam
注意:
-g
参数用的是create_Drop-seq_reference_metadata.sh
生成的STAR 的index的路径
-r
参数用的是create_Drop-seq_reference_metadata.sh
生成的新的fasta.gz
Drop-seq_alignment.sh
脚本的最终输出是经过处理后的 BAM 文件,其中每条 read 都带有 UMIs 和 barcodes 信息。但是这个 BAM 文件并不是最终的基因表达矩阵文件。
第三步: 表达矩阵的生成
要想得到表达矩阵,需要用 DigitalExpression
工具来从 BAM 文件中提取每个基因的表达量信息,生成一个基因表达矩阵文件。具体来说,可以运行以下命令来生成基因表达矩阵文件:
/Drop-seq_tools-2.5.1/DigitalExpression \
--INPUT input.bam \
--OUTPUT output.matrix.txt \
--SUMMARY output.summary.txt \
--MIN_NUM_GENES_PER_CELL 200 \
--TMP_DIR .
其中,
`[input BAM file]` 是 `Drop-seq_alignment.sh` 生成的 BAM 文件,
`[output directory]` 是基因表达矩阵文件的输出路径,
`[output summary file]` 是 DigitalExpression 运行过程的日志文件,
`[number of cores to use]` 是使用的 CPU 核数,
`[minimum read mapping quality]` 是过滤 BAM 文件中的低质量 read 的参数,
`[cell barcode tag name]` 和 `[UMI tag name]` 是 BAM 文件中保存 cell barcode 和 UMI 信息的 tag 名称,
`[gene barcode tag name]` 是 Drop-seq 实验中加在 cDNA bead 上的 barcode tag 名称,
`[maximum edit distance to allow]` 是允许的最大编辑距离,
`[single or dual]` 指定是单端测序还是双端测序,
`[logging level]` 是日志输出的详细程度。
在运行 DigitalExpression 后,会生成一个基因表达矩阵文件和一个日志文件,基因表达矩阵文件中每行对应一个基因,每列对应一个细胞,记录了每个细胞中每个基因的表达量。
Drop-seq数据下游分析:
使用Seurat进行后续分析
Seurat支持多种数据格式的导入,包括10x Genomics的输出文件、Drop-seq的输出文件以及基因表达矩阵文件。
具体操作步骤如下:
- 安装Seurat包。在R控制台中输入以下代码:
install.packages("Seurat")
- 读入基因表达矩阵文件。如果基因表达矩阵文件是以.tab或.csv格式存储的,可以使用
read.table
或read.csv
函数读入。例如,如果基因表达矩阵文件名为counts.tab
,可以使用以下代码将其读入:
counts <- read.table("counts.tab", header = TRUE, row.names = 1, sep = "\t")
其中,header = TRUE
表示第一行是列名,row.names = 1
表示第一列是行名,sep = "\t"
表示使用制表符分隔。
- 将基因表达矩阵文件转换为Seurat对象。使用
CreateSeuratObject
函数将基因表达矩阵文件转换为Seurat对象。例如,可以使用以下代码将其转换:
seuratObj <- CreateSeuratObject(counts)
到此,将基因表达矩阵文件转换为Seurat对象后,可以对其进行下游分析,例如数据质控、细胞聚类、细胞亚群分析、基因差异表达分析等。
再后续的步骤就不在这里写了。
报错收集
- 在使用create_Drop-seq_reference_metadata.sh时报错Unrecognized VM option 'UseParallelOldGC'
Runtime.totalMemory()=2181038080
Unrecognized VM option 'UseParallelOldGC'
Did you mean '(+/-)UseParallelGC'? Error: Could not create the Java Virtual Machine.
Error: A fatal exception has occurred. Program will exit.
ERROR: non-zero exit status 1 executing /software/biosoftware/Drop-seq_tools-2.5.1/FilterGtf
GTF=/reference/Homo_sapiens.GRCh38.109.gtf SEQUENCE_DICTIONARY=./Homo_sapiens_genome_annotation.dict
OUTPUT=./Homo_sapiens_genome_annotation.gtf
解释:UseParallelOldGC
是一个Java虚拟机参数,用于启用并行的老年代垃圾收集器。然而,该参数并不是所有的Java虚拟机版本都支持。因此,建议尝试将该参数从命令中删除,使用默认的垃圾收集器,或者使用其他可用的垃圾收集器。将命令中的-XX:+UseParallelOldGC
选项更改为-XX:+UseParallelGC
,并再次运行命令即可。
参考文献:
Salomon R, Kaczorowski D, Valdes-Mora F, Nordon RE, Neild A, Farbehi N, Bartonicek N, Gallego-Ortega D. Droplet-based single cell RNAseq tools: a practical guide. Lab Chip. 2019 May 14;19(10):1706-1727. doi: 10.1039/c8lc01239c. PMID: 30997473.
网友评论