非模式生物构建10x单细胞转录组CellRanger参考文件

生信技能树

发布于 2024-04-19 17:05:29

990

发布于 2024-04-19 17:05:29

10X单细胞上游定量标准流程运行Cellranger定量需要对应的参考基因组文件以及其配套的基因组注释信息文件，如果是人类和小鼠，官网即可下载构建好的文件压缩包，详见：https://www.10xgenomics.com/support/software/cell-ranger/downloads#reference-downloads

正常走cellranger的定量流程即可，代码我已经是多次分享了。参考：

差不多几个小时就可以完成全部的样品的cellranger的定量流程。

但是对于其它物种，就需要我们自己单独去构建cellranger 所需的参考基因组及注释文件（详见：https://www.10xgenomics.com/support/software/cell-ranger/latest/tutorials/cr-tutorial-mr）。

通常我们会去Ensembl网站独立下载，然后构建，Ensembl 的 GTF 文件包含可选标签，这使过滤变得容易。如果 Ensembl 无法提供您感兴趣的物种，则其他来源的 GTF 和 FASTA 文件也可以使用。但请注意，注释文件需要时GTF 格式，而 GFF 格式不受支持。

首先下载基因组注释信息

基因组注释是对参考基因组序列上的功能区域进行标识和描述的过程，它包括基因的位置、结构、功能和其他生物学特征。

GTF文件：

GTF（General Transfer Format）文件是一种基于纯文本的基因组注释文件格式，广泛用于存储基因、转录本、外显子等的注释信息。
GTF文件中的每一行代表一个基因组上的注释特征，通常包括以下字段：
- seqname：染色体或扫描序列的名称。
- source：注释信息的来源，例如基因预测软件或数据库。
- feature：特征类型，如基因（gene）、转录本（transcript）、外显子（exon）等。
- start：特征在序列上的起始位置（1-based）。
- end：特征在序列上的结束位置。
- score：一个可选的分数，通常用于表示序列比对的匹配程度。
- strand：特征在序列上的链，通常是+（正链）或-（负链）。
- frame：在某些情况下，指定编码区的阅读框架（0、1或2）。
- attribute：一系列键值对，提供额外的信息，如基因ID、外显子编号、注释来源等。

GTF文件为研究人员提供了一种标准化的方式来描述和共享基因组注释数据，它支持基因表达分析、变异分析和调控元件的识别等多种应用。

我们首先去Ensembl （ https://www.ensembl.org/index.html ）官网找到待处理物种的基因组注释信息（gtf文件），然后下载

##下载注释文件
wget -c https://ftp.ensembl.org/pub/release-111/gtf/gallus_gallus_gca000002315v5/Gallus_gallus_gca000002315v5.GRCg6a.111.gtf.gz
##文件大小
18M 10月  9 16:44 Gallus_gallus_gca000002315v5.GRCg6a.111.gtf.gz 
##解压后大小
488M 10月  9 16:44 Gallus_gallus_gca000002315v5.GRCg6a.111.gtf

为了避免在进行单细胞RNA测序数据分析时出现的多重映射（multi-mapped）问题，需要过滤一下GTF文件

详细信息可参考

Check for multi-mapped reads
Which reads are considered for UMI counting by Cell Ranger?

##过滤GTF
./cellranger-7.1.0/bin/cellranger mkgtf \
  Gallus_gallus_gca000002315v5.GRCg6a.111.gtf \
  Gallus_gallus_gca000002315v5.GRCg6a.111.filtered.gtf \
  --attribute=gene_biotype:protein_coding \
  --attribute=gene_biotype:lncRNA \
  --attribute=gene_biotype:antisense \
  --attribute=gene_biotype:IG_LV_gene \
  --attribute=gene_biotype:IG_V_gene \
  --attribute=gene_biotype:IG_V_pseudogene \
  --attribute=gene_biotype:IG_D_gene \
  --attribute=gene_biotype:IG_J_gene \
  --attribute=gene_biotype:IG_J_pseudogene \
  --attribute=gene_biotype:IG_C_gene \
  --attribute=gene_biotype:IG_C_pseudogene \
  --attribute=gene_biotype:TR_V_gene \
  --attribute=gene_biotype:TR_V_pseudogene \
  --attribute=gene_biotype:TR_D_gene \
  --attribute=gene_biotype:TR_J_gene \
  --attribute=gene_biotype:TR_J_pseudogene \
  --attribute=gene_biotype:TR_C_gene


## 按最新的教程其实也可以只保留编码蛋白质的基因的信息
cellranger mkgtf \
    Danio_rerio.GRCz11.105.gtf \
    Danio_rerio.GRCz11.105.filtered.gtf \
    --attribute=gene_biotype:protein_coding

上面的mkgtf操作到底是过滤掉了什么？

那么这一步做了什么，参考基因组注释文件有哪些变化，我们可以来详细看一下，首先对它参考基因组注释下载的原始文件进行基础的统计

 cat  Gallus_gallus_gca000002315v5.GRCg6a.111.gtf| awk '$3 == "gene" {print}'|cut -f 9 |perl -alne '{/gene_biotype "([^"]+)";/;print $1}' |sort |uniq  -c |sort -k1n

结果如下

     2 Mt_rRNA
      2 ribozyme
      2 vault_RNA
      2 Y_RNA
      5 misc_RNA
     13 processed_pseudogene
     15 scaRNA
     22 Mt_tRNA
     43 pseudogene
     74 snRNA
    108 rRNA
    200 snoRNA
    850 miRNA
  12447 lncRNA
  17077 protein_coding

然后统计一下过滤后的文件：

cat  Gallus_gallus_gca000002315v5.GRCg6a.111.filtered.gtf| awk '$3 == "gene" {print}'| cut -f 9 |perl -alne '{/gene_biotype "([^"]+)";/;print $1}' |sort |uniq  -c |sort -k1n

过滤后的结果统计：

  12447 lncRNA
  17077 protein_coding

可以看到，只剩下lncRNA和protein_coding。而参数中的antisense 、 IG_LV_gene 等，是参考基因组注释文件中本身就没有的。

其次是下载参考基因组FA文件

参考基因组是一个物种的代表性DNA序列，它通常是通过整合来自多个个体的遗传信息而构建的，旨在反映该物种的遗传多样性。参考基因组被用作比较基因组分析、基因发现、变异检测和基因功能研究的基础。

FA文件：

FA文件是FASTA格式的一种变体，通常用于存储参考基因组序列。
FASTA格式是一种生物序列文件格式，其中每个序列以一个以大于号（>）开头的标题行开始，后面跟着序列本身的一行或多行。
标题行通常包含序列的名称、来源和其他相关信息。
FA文件中的每个序列代表基因组中的一个染色体或片段，文件中可能包含整个基因组的所有染色体。

## 参考基因组下载
wget -c https://ftp.ensembl.org/pub/release-111/fasta/gallus_gallus_gca000002315v5/dna/Gallus_gallus_gca000002315v5.GRCg6a.dna.toplevel.fa.gz

##文件大小
306M 10月  4 19:46 Gallus_gallus_gca000002315v5.GRCg6a.dna.toplevel.fa.gz

##解压后大小
1.1G 10月  4 19:46 Gallus_gallus_gca000002315v5.GRCg6a.dna.toplevel.fa

构建CellRanger自定义参考

有了上面的两个文件，剩下的就很简单，执行cellranger mkref 即可

cellranger-7.1.0/bin/cellranger mkref \
  --genome=Gallus \
  --fasta=Gallus_gallus_gca000002315v5.GRCg6a.dna.toplevel.fa \
  --genes=Gallus_gallus_gca000002315v5.GRCg6a.111.filtered.gtf \
  --ref-version=1.0.0