image.png

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.png

Drop-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.3/Drop-seq_tools-2.5.3.zip
unzip Drop-seq_tools-2.5.3.zip

安装 bgzip

bgzip是一个常用的将文本压缩成 BGZF 格式的命令行工具，通常用于对大规模的基因组数据进行压缩，比如 VCF 文件。

以下是安装 bgzip 的步骤：

1. 安装 zlib 库：bgzip 依赖于 zlib 库，需要先安装该库。在 Ubuntu 系统下可以使用以下命令安装：

sudo apt-get update
sudo apt-get install zlib1g-dev

在其他 Linux 发行版中，可以使用相应的包管理器安装 zlib 库。

2. 下载安装 htslib 库：bgzip 是 htslib 库的一部分，需要先安装 htslib 库。可以从 htslib 的官方网站（https://github.com/samtools/htslib/releases）下载最新版本的源代码，解压后进入该目录，使用以下命令进行编译和安装：

make
sudo make install

如果编译过程中出现错误，可能需要安装一些其他的库和工具，如 OpenSSL、pkg-config 等。可以根据提示安装相应的库和工具。

3. 安装 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个文件

在分析流程之前，需要准备好参考基因组和转录本注释文件，这些文件可以从多个资源中获取，例如：

1. Ensembl（https://www.ensembl.org/index.html）和UCSC（https://genome.ucsc.edu/）等数据库提供了大量的参考基因组和注释文件，可以根据需要下载。

2. NCBI提供的 RefSeq 数据库（https://www.ncbi.nlm.nih.gov/refseq/）也提供了参考基因组和注释文件，可以在 NCBI 的 FTP 网站（ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/）上找到并下载。

3. 如果需要使用人类基因组的参考文件，可以考虑使用 GATK 团队提供的资源包（https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle），其中包含了多个版本的参考基因组和注释文件。

4. 10X Genomics 也提供了已经处理好的基因组和注释文件。在其官网可以获得。

一般来说，参考基因组文件是一个 fasta 格式的文件，包含了染色体序列和对应的基因组坐标信息。转录本注释文件可以是 gtf 或 gff3 格式的文件，包含了基因和转录本的注释信息，以及它们在基因组上的位置信息。

因为这里涉及到一些meta信息，我这里从Ensembl下载，下载sm_primary_assemble的fasta，或者 primary_assemble.fa。从头构建一套用于dropseq 分析的文件。

image-20230312011633163.png

mkdir 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 包含了以下几个步骤：

1. TagBamWithReadSequenceExtended：将每个read的UMI和barcode信息添加到bam文件的read名中。
2. FilterBAM：用于过滤低质量的reads和reads mapping到基因组外的区域。
3. TrimStartingSequence：用于去除 read 的起始序列。
4. PolyATrimmer：用于去除read的3'端的polyA序列。
5. TagReadWithInterval:TagReadWithInterval函数将每个比对结果中的Read映射到参考基因组上的转录本进行分类，并将分类结果添加到比对结果的SAM格式文件中。
6. TagReadWithGeneFunction:根据所属基因的注释信息（如基因名、基因类型、外显子位置等），将注释信息添加到比对结果的SAM格式文件中的XF标签中。
7. DetectBeadSubstitutionErrors:检测Drop-seq实验中可能存在的珠子替换错误（Bead Substitution Error），以保证单细胞RNA测序的准确性。
8. 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的输出文件以及基因表达矩阵文件。

具体操作步骤如下：

1. 安装Seurat包。在R控制台中输入以下代码：

install.packages("Seurat")

2. 读入基因表达矩阵文件。如果基因表达矩阵文件是以.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"表示使用制表符分隔。

3. 将基因表达矩阵文件转换为Seurat对象。使用CreateSeuratObject函数将基因表达矩阵文件转换为Seurat对象。例如，可以使用以下代码将其转换：

seuratObj <- CreateSeuratObject(counts)

到此，将基因表达矩阵文件转换为Seurat对象后，可以对其进行下游分析，例如数据质控、细胞聚类、细胞亚群分析、基因差异表达分析等。

再后续的步骤就不在这里写了。

报错收集

1. 在使用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.