新冠病毒、混合物种如何构建单细胞参考基因组？

在前面三篇文章中，我们着重介绍了 Cell Ranger 进行基因比对的理论模型【Cell Ranger 知多少（上）】，如何识别高质量细胞算法【Cell Ranger 知多少（中）】以及实操对细胞和基因的定量分析的代码以及对应的结果解读【Cell Ranger 知多少（下）】。系统地学习了这些后，大家在工作中是否更得心应手了呢？然而 10X Genomics 只为 Cell Ranger 提供了预先构建的人类和小鼠基因组参考，其他的物种则需要研究人员自己构建，因此，参考基因组的构建部分值得一讲。那么，话题来了：特殊物种的参考基因组如何构建？含有病毒序列的参考基因组如何构建？混合物种的参考基因组又当如何处理? 接下来就跟着小编来学习学习吧~

使用 cellranger mkref 构建参考基因组，需要 3 步：

1）获得输入文件

在分析之前需要找到对应的物种的参考基因组 FASTA 和 GTF 文件，10X Genomics 官方建议从 Ensembl 数据库中获得。如果感兴趣的物种无法从 Ensembl 中获得，那么来自其他来源的 GTF 和 FASTA 文件也可以，但是需要注意的是，构建基因组仅支持使用 GTF 文件，不支持 GFF 文件。如果物种的参考基因组中没有 GTF 文件，也不用担心，可以通过 GFF 文件转成 GTF 文件，可使用命令 gffread *gff -T -o *gtf。

以斑马鱼为例：

# 获取 FASTA 文件
wget ftp://ftp.ensembl.org/pub/release-98/fasta/danio_rerio/dna/Danio_rerio.GRCz11.dna.primary_assembly.fa.gz
gunzip Danio_rerio.GRCz11.dna.primary_assembly.fa.gz

# 获取 GTF 文件
wget ftp://ftp.ensembl.org/pub/release-99/gtf/danio_rerio/Danio_rerio.GRCz11.99.chr.gtf.gz
gunzip  Danio_rerio.GRCz11.99.chr.gtf.gz

2）过滤GTF

GTF 文件包含与蛋白质编码基因模型重叠的 non-polyA 转录本条目。由于注释的重叠，这些条目可能导致 reads 被标记为映射到多个基因(多映射)。如果 reads 被标记为多映射，它们将不被计算【 Cell Ranger 知多少（上）】。要从 GTF 中删除这些条目，向 cellranger mkgtf 命令添加这个过滤参数: --attribute=gene_biotype:protein_coding。如果使用的是一个不包含 gene_biotype 属性或缺少其他条目的 GTF 文件，也不用担心，可能仍然有足够的信息来建立参考基因组。一个最小的 GTF 文件只需要包含蛋白编码基因的外显子特征，即 GTF 文件第 3 列需存在 "exon" 特征类型。

cellranger mkgtf \\
    Danio_rerio.GRCz11.99.chr.gtf \\
    Danio_rerio.GRCz11.99.chr.filtered.gtf \\
    --attribute=gene_biotype:protein_coding

3）构建参考基因组

可通过 cellranger mkref --help 查看参数使用方法。

cellranger mkref 的参数作用

• genome 唯一的基因组名称，用于命名输出文件夹，应仅包含字母、数字字符和英文句点，连字符和下划线字符 [a-zA-Z0-9 _-] +。可通过多次指定该参数指定多个基因组。

• fasta 基因组参考 FASTA 文件的路径，可通过多次指定该参数指定多个基因组。

• genes 基因组包含基因注释的 GTF 文件的路径，可通过多次指定该参数指定多个基因组注释文件。

• nthreads （可选）使用 STAR 生成基因组索引使用的线程数，cellranger v5 此参数有 bug, 暂时不可用。

• memgb （可选）使用 STAR 比对时使用的最大内存（GB），默认值为 16。请注意，比对期间所使用的内存量必须大于输入 FASTA 文件的大小。

• ref-version （可选）参考基因组版本号。

  
  

执行代码:

cellranger mkref \\
    --genome=Danio.rerio_genome \\ # 构建好的参考基因组目录名称
    --fasta=Danio_rerio.GRCz11.dna.primary_assembly.fa \\ # FASTA 文件
    --genes=Danio_rerio.GRCz11.99.chr.filtered.gtf # GTF 文件

成功构建的标志：

# >>> Reference successfully created! <<<

# You can now specify this reference on the command line:
# cellranger --transcriptome=Danio.rerio_genome ...

此时，一个特殊物种的参考基因组就构建好了。

当研究对象是 COVID19+ 病人，进行 10X Genomics 转录组数据分析前需要构建带有  SARS-CoV-2 的人类参考基因组，那么如何将新冠病毒 SARS-CoV-2 添加到待构建的人类 GRCh38 参考基因组中呢？

首先需要下载到与研究对象匹配的 SARS-CoV-2 的序列，该序列可以在 NCBI 的病毒数据库（点击“阅读原文”即可获取）中选择。当找到与之匹配的 SARS-CoV-2 基因组后，首先下载基因组  FASTA 文件，将其添加到人类参考基因组中；然后需要为 SARS-CoV-2 制作一个 GTF 文件，建议将整个病毒基因组注释为一个单一的"基因"，特征类型标记为 "exon"；最后将制作好的 GTF 文件添加到人类 GTF 文件中。

1）从 NCBI 上下载 SARS-CoV-2  NC_045512 的 FASTA 文件中制作 GTF 文件。

## 查看 FASTA 文件数据
$ head -n 3 NC_045512.fa
# >NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome
# ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
# CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC

## 获取序列长度，该长度值在下一步制作 GTF 文件中使用
$ cat NC_045512.fa | grep -v "^>" | tr -d "\\n" | wc -c
# 29903

## 制作 GTF 文件
$ echo -e ''NC_045512.2\\tRefSeq\\texon\\t1\\t29903\\t.\\t+\\t.\\tgene_id "NC_045512.2"; transcript_id "NC_045512.2"; gene_name "NC_045512.2";'' > NC_045512.gtf

## 查看 GTF 文件
$ cat NC_045512.gtf
# NC_045512.2    RefSeq  exon    1   29903   .   +   .   gene_id "NC_045512.2"; transcript_id "NC_045512.2"; gene_name "NC_045512.2";

tips: 整理病毒 GTF 文件时要求 GTF 文件至少在第 3 列含有 exon 可识别类型，对应的第 9 列至少需要包含 transcript_id 和 gene_id 关键词，以便 cellranger mkref pipeline 识别。这两点是该步骤所需 GTF 文件的最小格式。

2）将 SARS-CoV-2 的 FASTA 文件和 GTF 文件分别添加到人的 GCRh38 FASTA 文件和 GTF 文件中，构建新的参考基因组文件。

## 处理 FASTA 文件
$ cat GRCh38-2020-A.fa NC_045512.fa > NC_045512_GRCh38-2020-A.fa

## 处理 GTF 文件
$ cat GRCh38-2020-A-genes.gtf  NC_045512.gtf > NC_045512_GRCh38-2020-A-genes.gtf

## 通过检测各文件行数进行 check 
$ grep -c ">" GRCh38-2020-A.fa
# 194
$ grep -c ">" NC_045512_GRCh38-2020-A.fa
# 195

$ wc -l GRCh38-2020-A-genes.gtf
# 2765974 GRCh38-2020-A-genes.gtf
$ wc -l NC_045512_GRCh38-2020-A-genes.gtf
# 2765975 NC_045512_GRCh38-2020-A-genes.gtf

当新的 FASTA 文件和 GTF 文件检查没问题后执行 cellranger mkref pipeline。

3）构建含有病毒的参考基因组

cellranger mkref \\
    --genome=GRCh38-2020-A_NC_045512 \\      
    --fasta=NC_045512_GRCh38-2020-A.fa \\
    --genes=NC_045512_GRCh38-2020-A-genes.gtf 

这样，含有病毒的参考基因组就构建好啦，是不是也不难哟~

如果我们需要构建人和小鼠的混合参考基因组，怎么处理呢？由于人和小鼠的染色体号和基因名大部分是重叠的，因此不能将两者的 FASTA 和 GTF 文件直接合并，好像有点麻烦。。。经过一番探索小编发现，cellranger mkref pipeline 可以直接完成这一需求，不需要进行任何文件处理。简直太惊喜了，一起来看看吧~

执行代码:

cellranger mkref \\
    --genome=hg19 --fasta=hg19.fa --genes=hg19.gtf \\
    --genome=mm10 --fasta=mm10.fa --genes=mm10.gtf

这里虽然和单个物种是类似的，但需要注意的是 --genome, --fasta, --genes 这 3 个参数对某一物种来说，是连续的。例如第一个 --genome 需要对应于第一个 --fasta 以及第一个 --genes，以此类推。当任务完成后，输出结果的目录为 hg19_and_mm10。

以上就是特殊情况（特殊物种，含有病毒，混合物种）的参考基因组构建的介绍，您掌握了吗？如果有疑问，欢迎在下方留言，小编会及时为您答疑解惑~

END

单细胞生信分析部  撰文

本文系欧易生物原创

欢迎转发到朋友圈

转载请注明文本转自欧易生物

欧易生物 | oebiotech

www.oebiotech.com/cn

☏ 400-808-5350 / 021-34781616