DropSeq 簡(jiǎn)介
Drop-seq技術(shù)最初由麻省理工學(xué)院的McCarroll實(shí)驗(yàn)室開(kāi)發(fā)邻梆，于2015年首次發(fā)表在《Cell》雜志上爆侣。自此以后熄云，這種技術(shù)已被廣泛應(yīng)用于許多研究領(lǐng)域膜赃，包括神經(jīng)科學(xué)豌注、發(fā)育生物學(xué)和腫瘤學(xué)等伤塌。
Drop-seq是一種單細(xì)胞RNA測(cè)序技術(shù)，通過(guò)在微滴中捕獲單個(gè)細(xì)胞并進(jìn)行RNA擴(kuò)增轧铁，從而獲得單個(gè)細(xì)胞的轉(zhuǎn)錄組數(shù)據(jù)每聪。該技術(shù)通過(guò)微滴分離單個(gè)細(xì)胞并將細(xì)胞裂解，隨后在微滴中添加反轉(zhuǎn)錄酶和一種稱為“barcode beads”的特殊珠子属桦，這些珠子上有一個(gè)獨(dú)特的序列標(biāo)識(shí)符熊痴。這些珠子能夠被反轉(zhuǎn)錄酶轉(zhuǎn)錄RNA，并在RNA末端添加一段序列聂宾，從而將每個(gè)RNA分子與其來(lái)源細(xì)胞進(jìn)行標(biāo)記果善。然后系谐，通過(guò)PCR擴(kuò)增這些標(biāo)記的RNA分子巾陕，將其構(gòu)建成文庫(kù)并進(jìn)行高通量測(cè)序讨跟，以獲得單個(gè)細(xì)胞的轉(zhuǎn)錄組數(shù)據(jù)。
DropSeq技術(shù)與其他單細(xì)胞RNA測(cè)序技術(shù)相比鄙煤，具有高通量晾匠、高靈敏度和高特異性等優(yōu)點(diǎn)，能夠在細(xì)胞數(shù)量有限的情況下實(shí)現(xiàn)單細(xì)胞RNA測(cè)序梯刚，并廣泛應(yīng)用于單細(xì)胞轉(zhuǎn)錄組學(xué)研究凉馆。
DropSeq 文庫(kù)結(jié)構(gòu)：
read 1中包含12bp的barcode，8bp的UMI亡资，read2中是捕獲的RNA分子澜共。
Drop-seq數(shù)據(jù)前期分析
軟件安裝
下面列出的軟件如果你已經(jīng)安裝了，則跳過(guò)此步驟锥腻。
安裝DropSeq Tools
DropSeq 數(shù)據(jù)通常使用DropSeq Tools做基礎(chǔ)分析嗦董，可以進(jìn)行 reads 到 barcodes 和 UMIs 的mapping瘦黑、去除低質(zhì)量 reads、去除 PCR 重復(fù)幸斥、生成 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是一個(gè)常用的將文本壓縮成 BGZF 格式的命令行工具甲葬，通常用于對(duì)大規(guī)模的基因組數(shù)據(jù)進(jìn)行壓縮，比如 VCF 文件演顾。
以下是安裝 bgzip 的步驟：
- 安裝 zlib 庫(kù)：
bgzip依賴于 zlib 庫(kù)隅居，需要先安裝該庫(kù)。在 Ubuntu 系統(tǒng)下可以使用以下命令安裝：
sudo apt-get update
sudo apt-get install zlib1g-dev
在其他 Linux 發(fā)行版中胎源，可以使用相應(yīng)的包管理器安裝 zlib 庫(kù)。
- 下載安裝 htslib 庫(kù)：
bgzip是 htslib 庫(kù)的一部分涕蚤，需要先安裝 htslib 庫(kù)宪卿。可以從 htslib 的官方網(wǎng)站（https://github.com/samtools/htslib/releases）下載最新版本的源代碼万栅，解壓后進(jìn)入該目錄佑钾，使用以下命令進(jìn)行編譯和安裝：
make
sudo make install
如果編譯過(guò)程中出現(xiàn)錯(cuò)誤，可能需要安裝一些其他的庫(kù)和工具烦粒，如 OpenSSL休溶、pkg-config 等代赁。可以根據(jù)提示安裝相應(yīng)的庫(kù)和工具兽掰。
- 安裝 bgzip：安裝完 htslib 庫(kù)后芭碍，即可使用以下命令安裝
bgzip：
git clone https://github.com/samtools/tabix.git
cd tabix && make
sudo make install
安裝完成后，可以使用 bgzip 命令進(jìn)行文件壓縮和解壓縮孽尽。
安裝STAR
DropSeq使用STAR軟件作為比對(duì)軟件窖壕，所以需要下載安裝STAR
在latest release頁(yè)面找到最新版的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
準(zhǔn)備2個(gè)文件
在分析流程之前杉女，需要準(zhǔn)備好參考基因組和轉(zhuǎn)錄本注釋文件瞻讽，這些文件可以從多個(gè)資源中獲取，例如：
Ensembl（https://www.ensembl.org/index.html）和UCSC（https://genome.ucsc.edu/）等數(shù)據(jù)庫(kù)提供了大量的參考基因組和注釋文件宠纯，可以根據(jù)需要下載卸夕。
NCBI提供的 RefSeq 數(shù)據(jù)庫(kù)（https://www.ncbi.nlm.nih.gov/refseq/）也提供了參考基因組和注釋文件，可以在 NCBI 的 FTP 網(wǎng)站（ftp://ftp.ncbi.nlm.nih.gov/genomes/refseq/）上找到并下載婆瓜。
如果需要使用人類基因組的參考文件快集，可以考慮使用 GATK 團(tuán)隊(duì)提供的資源包（https://gatk.broadinstitute.org/hc/en-us/articles/360035890811-Resource-bundle），其中包含了多個(gè)版本的參考基因組和注釋文件廉白。
10X Genomics 也提供了已經(jīng)處理好的基因組和注釋文件个初。在其官網(wǎng)可以獲得。
一般來(lái)說(shuō)猴蹂，參考基因組文件是一個(gè) fasta 格式的文件院溺，包含了染色體序列和對(duì)應(yīng)的基因組坐標(biāo)信息。轉(zhuǎn)錄本注釋文件可以是 gtf 或 gff3 格式的文件磅轻，包含了基因和轉(zhuǎn)錄本的注釋信息珍逸，以及它們?cè)诨蚪M上的位置信息。
因?yàn)檫@里涉及到一些meta信息聋溜，我這里從Ensembl下載，下載sm_primary_assemble的fasta漱病，或者 primary_assemble.fa杨帽。從頭構(gòu)建一套用于dropseq 分析的文件嗤军。
mkdir reference # 隨便創(chuàng)建個(gè)文件夾放一下
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
在下載參考基因組和注釋文件后，可以將它們作為參數(shù)傳遞給 create_Drop-seq_reference_metadata.sh 腳本沿量，并根據(jù)腳本提示的要求提供相應(yīng)的文件路徑即可朴则。
生成一些 meta 文件
在 Drop-seq 流程中，需要使用參考基因組和轉(zhuǎn)錄本注釋文件（gtf）對(duì)測(cè)序數(shù)據(jù)進(jìn)行比對(duì)和定量汹想。這個(gè)腳本的作用就是幫助用戶生成這些必要的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 如果在環(huán)境變量中槽唾，就可以省略這三個(gè)參數(shù)庞萍。
create_Drop-seq_reference_metadata.sh 腳本通常會(huì)為參考基因組創(chuàng)建 .dict 文件忘闻，該文件包含參考基因組的字典信息，是 Drop-seq 流程所需的重要參考文件之一私恬。會(huì)重新生成一下fasta炼吴，gtf缺厉，一些intervals, 還會(huì)構(gòu)建STAR index提针。
在 create_Drop-seq_reference_metadata.sh 腳本中曹傀，picard CreateSequenceDictionary 命令用于創(chuàng)建 .dict 文件皆愉。該命令會(huì)讀取參考基因組的 .fasta 文件艇抠，并基于參考基因組序列的名稱家淤、長(zhǎng)度絮重、描述等信息生成 .dict 文件歹苦。需要注意的是，這里使用了 Picard 工具中的 CreateSequenceDictionary 命令來(lái)創(chuàng)建 .dict 文件狠角，而不是使用 STAR 工具創(chuàng)建丰歌。這是因?yàn)?Picard 工具在處理字典信息時(shí)更加嚴(yán)格和精確动遭，可以確保 Drop-seq 流程的正確性神得。
這一步生成的這些文件要用在后面alignment中作為參數(shù)，而不用剛剛下載的原始reference宵蕉。
生成表達(dá)矩陣
第一步：fqtobam
fastq 文件必須要轉(zhuǎn)換為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
# 上述寫法雖然可以運(yùn)行改执，但是picard更新了語(yǔ)法，
# 我們與時(shí)俱進(jìn)衬横，用下面的新語(yǔ)法更時(shí)尚：
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 是一個(gè)all-in-one的腳本冕香，它將 Drop-seq 實(shí)驗(yàn)中的多個(gè)處理步驟整合到一個(gè)腳本中后豫，簡(jiǎn)化了數(shù)據(jù)處理的流程挫酿。
具體來(lái)說(shuō)，Drop-seq_alignment.sh 包含了以下幾個(gè)步驟：
-
TagBamWithReadSequenceExtended：將每個(gè)read的UMI和barcode信息添加到bam文件的read名中惫霸。 -
FilterBAM：用于過(guò)濾低質(zhì)量的reads和reads mapping到基因組外的區(qū)域壹店。 -
TrimStartingSequence：用于去除 read 的起始序列芝加。 -
PolyATrimmer：用于去除read的3'端的polyA序列藏杖。 -
TagReadWithInterval:TagReadWithInterval函數(shù)將每個(gè)比對(duì)結(jié)果中的Read映射到參考基因組上的轉(zhuǎn)錄本進(jìn)行分類，并將分類結(jié)果添加到比對(duì)結(jié)果的SAM格式文件中点寥。 -
TagReadWithGeneFunction:根據(jù)所屬基因的注釋信息（如基因名敢辩、基因類型弟疆、外顯子位置等），將注釋信息添加到比對(duì)結(jié)果的SAM格式文件中的XF標(biāo)簽中。 -
DetectBeadSubstitutionErrors:檢測(cè)Drop-seq實(shí)驗(yàn)中可能存在的珠子替換錯(cuò)誤（Bead Substitution Error）嘀略，以保證單細(xì)胞RNA測(cè)序的準(zhǔn)確性。 -
DetectBeadSynthesisErrors：檢測(cè)每個(gè)UMI的錯(cuò)配率咒程，并將潛在的bead合成錯(cuò)誤的UMI標(biāo)記為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參數(shù)用的是create_Drop-seq_reference_metadata.sh生成的STAR 的index的路徑
-r參數(shù)用的是create_Drop-seq_reference_metadata.sh生成的新的fasta.gz
Drop-seq_alignment.sh 腳本的最終輸出是經(jīng)過(guò)處理后的 BAM 文件奶段，其中每條 read 都帶有 UMIs 和 barcodes 信息痹籍。但是這個(gè) BAM 文件并不是最終的基因表達(dá)矩陣文件。
第三步： 表達(dá)矩陣的生成
要想得到表達(dá)矩陣棺克，需要用 DigitalExpression 工具來(lái)從 BAM 文件中提取每個(gè)基因的表達(dá)量信息娜谊，生成一個(gè)基因表達(dá)矩陣文件斤讥。具體來(lái)說(shuō)，可以運(yùn)行以下命令來(lái)生成基因表達(dá)矩陣文件：
/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]` 是基因表達(dá)矩陣文件的輸出路徑钳踊，
`[output summary file]` 是 DigitalExpression 運(yùn)行過(guò)程的日志文件勿侯，
`[number of cores to use]` 是使用的 CPU 核數(shù)，
`[minimum read mapping quality]` 是過(guò)濾 BAM 文件中的低質(zhì)量 read 的參數(shù)祭埂，
`[cell barcode tag name]` 和 `[UMI tag name]` 是 BAM 文件中保存 cell barcode 和 UMI 信息的 tag 名稱蛆橡，
`[gene barcode tag name]` 是 Drop-seq 實(shí)驗(yàn)中加在 cDNA bead 上的 barcode tag 名稱泰演，
`[maximum edit distance to allow]` 是允許的最大編輯距離，
`[single or dual]` 指定是單端測(cè)序還是雙端測(cè)序藐握，
`[logging level]` 是日志輸出的詳細(xì)程度猾普。
在運(yùn)行 DigitalExpression 后缔御，會(huì)生成一個(gè)基因表達(dá)矩陣文件和一個(gè)日志文件，基因表達(dá)矩陣文件中每行對(duì)應(yīng)一個(gè)基因笤成，每列對(duì)應(yīng)一個(gè)細(xì)胞炕泳，記錄了每個(gè)細(xì)胞中每個(gè)基因的表達(dá)量上祈。
Drop-seq數(shù)據(jù)下游分析：
使用Seurat進(jìn)行后續(xù)分析
Seurat支持多種數(shù)據(jù)格式的導(dǎo)入登刺，包括10x Genomics的輸出文件、Drop-seq的輸出文件以及基因表達(dá)矩陣文件皇耗。
具體操作步驟如下：
- 安裝Seurat包郎楼。在R控制臺(tái)中輸入以下代碼：
install.packages("Seurat")
- 讀入基因表達(dá)矩陣文件窒悔。如果基因表達(dá)矩陣文件是以.tab或.csv格式存儲(chǔ)的，可以使用
read.table或read.csv函數(shù)讀入阶界。例如膘融，如果基因表達(dá)矩陣文件名為counts.tab，可以使用以下代碼將其讀入：
counts <- read.table("counts.tab", header = TRUE, row.names = 1, sep = "\t")
其中，header = TRUE表示第一行是列名屯耸，row.names = 1表示第一列是行名蹭劈，sep = "\t"表示使用制表符分隔铺韧。
- 將基因表達(dá)矩陣文件轉(zhuǎn)換為Seurat對(duì)象。使用
CreateSeuratObject函數(shù)將基因表達(dá)矩陣文件轉(zhuǎn)換為Seurat對(duì)象塔逃。例如湾盗，可以使用以下代碼將其轉(zhuǎn)換：
seuratObj <- CreateSeuratObject(counts)
到此立轧，將基因表達(dá)矩陣文件轉(zhuǎn)換為Seurat對(duì)象后，可以對(duì)其進(jìn)行下游分析帐萎，例如數(shù)據(jù)質(zhì)控疆导、細(xì)胞聚類瑰艘、細(xì)胞亞群分析紫新、基因差異表達(dá)分析等。
再后續(xù)的步驟就不在這里寫了囤耳。
報(bào)錯(cuò)收集
- 在使用create_Drop-seq_reference_metadata.sh時(shí)報(bào)錯(cuò)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是一個(gè)Java虛擬機(jī)參數(shù)充择，用于啟用并行的老年代垃圾收集器。然而宰僧，該參數(shù)并不是所有的Java虛擬機(jī)版本都支持琴儿。因此嘁捷，建議嘗試將該參數(shù)從命令中刪除雄嚣，使用默認(rèn)的垃圾收集器，或者使用其他可用的垃圾收集器鼓鲁。將命令中的-XX:+UseParallelOldGC選項(xiàng)更改為-XX:+UseParallelGC坐桩，并再次運(yùn)行命令即可封锉。
參考文獻(xiàn)：
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.
