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分子。
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 的步骤：
- 安装 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的fasta，或者 primary_assemble.fa。从头构建一套用于dropseq 分析的文件。
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 包含了以下几个步骤：
-
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.
