美文网首页生信笔记单细胞组学
单细胞 Drop-seq数据分析教程

单细胞 Drop-seq数据分析教程

作者: 11的雾 | 来源:发表于2023-03-16 17:43 被阅读0次
    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.2/Drop-seq_tools-2.5.1.zip
    unzip Drop-seq_tools-2.5.1.zip
    

    安装 bgzip

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

    以下是安装 bgzip 的步骤:

    1. 安装 zlib 库:bgzip 依赖于 zlib 库,需要先安装该库。在 Ubuntu 系统下可以使用以下命令安装:
    sudo apt-get update
    sudo apt-get install zlib1g-dev
    

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

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

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

    1. 安装 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的faata,或者 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")
    
    1. 读入基因表达矩阵文件。如果基因表达矩阵文件是以.tab或.csv格式存储的,可以使用read.tableread.csv函数读入。例如,如果基因表达矩阵文件名为counts.tab,可以使用以下代码将其读入:
    counts <- read.table("counts.tab", header = TRUE, row.names = 1, sep = "\t")
    

    其中,header = TRUE表示第一行是列名,row.names = 1表示第一列是行名,sep = "\t"表示使用制表符分隔。

    1. 将基因表达矩阵文件转换为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.

    相关文章

      网友评论

        本文标题:单细胞 Drop-seq数据分析教程

        本文链接:https://www.haomeiwen.com/subject/irddrdtx.html