特殊物种cellranger基因组质量评估

男，

一个长大了才会遇到的帅哥，

稳健，潇洒，大方，靠谱。

一段生信缘，一棵技能树，

一枚大型测序工厂的螺丝钉，

一个随机森林中提灯觅食的津门旅客。

---我们实验室是研究雪莲的，可以做10X单细胞转录组吗？
···可以

--- 我们实验室前几年做了雪莲的基因组，没有发表，师兄做的，不知道质量怎么样，可以做10X单细胞转录组吗？
···可以

---我们实验室做的雪莲三代转录组，有一个基因组，基于这个可以做10X单细胞转录组吗？
···可以

所以说，基因组是生命科学实验室基础建设的一部分，在不远的将来，单细胞也会是。

要回答上述问题，首先要明白的一点就是：基因组是什么？

基因组主要有两个文件：

• fa序列文件

1>15 dna:chromosome chromosome:GRCh38:15:1:101991189:1 REF
2NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
3NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
4NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
5NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
6NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN


• gtf注释文件

1!genome-build GRCh38.p12
2#!genome-version GRCh38
3#!genome-date 2013-12
4#!genome-build-accession NCBI:GCA_000001405.27
5#!genebuild-last-updated 2018-01
61       havana  gene    29554   31109   .       +       .       gene_id "ENSG00000243485"; gene_version "5"; gene_name "MIR1302-2HG"; gene_source "havana"; gene_biotype "lincRNA"
71       havana  transcript      29554   31097   .       +       .       gene_id "ENSG00000243485"; gene_version "5"; transcript_id "ENST00000473358"; transcript_version "1"; gene_name "MIR1302-2HG"; gene_source "havana"; gene_biotype "lincRNA"; transcript_name "MIR1302-
81       havana  exon    29554   30039   .       +       .       gene_id "ENSG00000243485"; gene_version "5"; transcript_id "ENST00000473358"; transcript_version "1"; exon_number "1"; gene_name "MIR1302-2HG"; gene_source "havana"; gene_biotype "lincRNA"; transcript_name 
91       havana  exon    30564   30667   .       +       .       gene_id "ENSG00000243485"; gene_version "5"; transcript_id "ENST00000473358"; transcript_version "1"; exon_number "2"; gene_name "MIR1302-2HG"; gene_source "havana"; gene_biotype "lincRNA"; transcript_name 
101       havana  exon    30976   31097   .       +       .       gene_id "ENSG00000243485"; gene_version "5"; transcript_id "ENST00000473358"; transcript_version "1"; exon_number "3"; gene_name "MIR1302-2HG"; gene_source "havana"; gene_biotype "lincRNA"; transcript_name 
111       havana  transcript      30267   31109   .       +       .       gene_id "ENSG00000243485"; gene_version "5"; transcript_id "ENST00000469289"; transcript_ve
12


组装

序列文件就是基因组的序列以fa格式存储，这里我们看到在GRCh38版本中染色体两端加了很多N。
从序列文件我们可以得到什么？

• 组装水平：染色体，contig，还是scaffold水平？

• 组装质量评估：

动植物基因组de novo工作，其组装指标的好坏直接影响着整个基因组的质量。而评估基因组组装结果，contigN50和scaffoldN50是第一指标，即contig/ scaffoldN50：将contig/scaffold长度从长到短进行排序并累加，当累加和达到contig/scaffold总长度的50%的时候，最后参与加和的那一条contig/scaffold长度即为contig/ scaffoldN50的长度。一般来说，contig/scaffoldN50越长，表示组装结果越好。

但是，N50指标高就意味着组装结果就一定可靠吗？

不一定！将一些不相关的reads或者contig错误的连接为scaffold，一样可以达到很高的scaffoldN50。

目前高水平文章发表，组装指标固然是一方面，但真正决定文章发表档次的，是生物学故事是否足够完美，有亮点。我们知道，后续分析依赖的基础便是组装得到的基因组，因此，不可靠的组装结果，对基因组后续分析会造成很大的困扰，甚至会得出错误的生物学结论。

那么，如何才能检验一个基因组组装结果的可靠性呢？

1、 序列一致性评估：

基因组是通过reads组装得到，这一步，是将reads比到基因组上，验证reads对基因组的覆盖情况，用于评估组装的完整性以及测序的均匀性。较高的mapping rate（90%以上）以及coverage（95%以上）认为组装结果和reads有比较好的一致性。

2、 序列完整性评估：

所谓完整性评估，即评估组装得到的基因组对基因区的覆盖程度，一般需要借助RNA方面的证据进行评估，如EST数据或RNA reads。由于用来评估的RNA方面证据不同，得到的比例也会有差别。一般来说，50%的scaffold覆盖基因的95%以上，85%的scaffold覆盖基因的90%以上，认为组装较完整。

3、 准确性评估：

通过全长BAC序列，可以通过与组装结果的比对，对组装结果的正确性进行验证，从BAC序列和scaffold是否具有较好的一致性来判断组装质量。

4、 保守性基因评估：

即根据广泛存在于大量真核生物中的保守蛋白家族集合（248个core gene库），对组装得到基因组进行评估，评估组装基因组中的core gene的准确性和完整性。可以通过该物种和同源物种cegma的比例，判断保守基因组装情况。

有没有现成的方法来评估呢？

有的，LAI: 评估基因组质量一个标准

得到的LAI值按照如下评估标准进行分类：

Category	LAI	Examples
Draft	0 ≤ LAI < 10	Apple (v1.0), Cacao (v1.0)
Reference	10 ≤ LAI < 20	Arabidopsis (TAIR10), Grape (12X)
Gold	20 ≤ LAI	Rice (MSUv7), Maize (B73 v4)

注释

注释就是以位置信息来注明基因组的序列每一段都是什么功能（一种描述）。

那么，如何对基因组序列进行注释呢？基因组组装完成后，或者是完成了草图，就不可避免遇到一个问题，需要对基因组序列进行注释。注释之前首先得构建基因模型，有三种策略：

• 从头注释(de novo prediction)：通过已有的概率模型来预测基因结构，在预测剪切位点和UTR区准确性较低

• 同源预测(homology-based prediction)：有一些基因蛋白在相近物种间的保守型高，所以可以使用已有的高质量近缘物种注释信息通过序列联配的方式确定外显子边界和剪切位点

• 基于转录组预测(transcriptome-based prediction)：通过物种的RNA-seq数据辅助注释，能够较为准确的确定剪切位点和外显子区域。

在高通量测序的时代，基因组序列的获得已经不是难题了，但是每段序列的注释依然需要也是值得花一些精力的。

我的基因组可以做10X单细胞转录组了吗？

在对基因组有了基本的认识之后，我们来回答这个问题。

Cell Ranger uses an aligner called STAR, which peforms splicing-aware alignment of reads to the genome. Cell Ranger then uses the transcript annotation GTF to bucket the reads into exonic, intronic, and intergenic, and by whether the reads align (confidently) to the genome. A read is exonic if at least 50% of it intersects an exon, intronic if it is non-exonic and intersects an intron, and intergenic otherwise.

基本的注释信息：

Column	Name	Description
1	Chromosome	Must refer to a chromosome/contig in the genome fasta.
2	Source	Unused.
3	Feature	cellranger count only uses rows where this line is exon.
4	Start	Start position on the reference (1-based inclusive).
5	End	End position on the reference (1-based inclusive).
6	Score	Unused.
7	Strand	Strandedness of this feature on the reference: + or -.
8	Frame	Unused.
9	Attributes	A semicolon-delimited list of key-value pairs of the form key "value". The attribute keys transcript_id and gene_idare required; gene_name is optional and may be non-unique, but if present will be preferentially displayed in reports.

也就是注释信息中必须要有exon，transcript_id，gene_id ，这个是做10X单细胞转录组对一个基因组最基本的要求。能组装到染色体水平当然更好，组装不到的话也可以。

有了fa以及gtf文件之后，我们就可以用cellrang的mkerf流程来构建10X专用的参考基因组了：

1cellranger mkref --genome=output_genome --fasta=input.fa --genes=input.gtf


构建好之后，是这样的：

1genome_output/
2├── fasta
3│   └── genome.fa
4├── genes
5│   └── genes.gtf
6├── pickle
7│   └── genes.pickle
8├── reference.json
9└── star # STAR genome index folder


• For the genome sequence, include all major chromosomes, unplaced and unlocalized scaffolds, but do not include patches and alternative haplotypes.

• In Ensembl, the recommended genome file to download is annotated as "primary assembly." - In NCBI, it is "no alternative - analysis set."

• For the GTF file, genes must be annotated with feature type 'exon' (column 3). - Prior to mkref, GTF annotation files from Ensembl and NCBI are typically filtered with `mkgtf` to include only a subset of the annotated gene biotypes.

Creating a Reference Package with cellranger mkref

关于特殊物种细胞类型的注释

这很大程度上取决于我们的基因组注释情况，如果是斑马鱼这种模式生物，一般的研究者是很多的，文献检索是可以获得有益的背景（marker基因或者表达谱）知识用于细胞类型鉴定的：

有了marker基因或者表达谱细胞的定义就和人鼠的没有什么方法学上的区别了。

第二种，新鲜的基因组，用于定量转录组的基因功能还不清楚，基因名只是自定义的编号。这个有两个方法来做：

• 对基因做功能富集，看某群的差异基因在功能上富集到哪，根据功能结合生物学知识来做。

R包clusterProfiler的纯无参自定义物种注释的GO、KEGG富集分析及GSEA

• 同源基因。用同源基因来将特殊物种的基因与已知基因构建联系，自然界基因并不是每个物种一套特意的基因，有许多是同源的。

基本思路也是构建从已知到未知的证据链。

---

References

[1] LAI: 评估基因组质量一个标准: https://www.jianshu.com/p/7d794d22e0a0
[2] 如何对基因组序列进行注释: https://www.jianshu.com/p/931e9821c45a
[3] STAR: https://links.jianshu.com/go?to=https%3A%2F%2Fgithub.com%2Falexdobin%2FSTAR
[4] feature type 'exon' (column 3): https://links.jianshu.com/go?to=https%3A%2F%2Fsupport.10xgenomics.com%2Fsingle-cell-gene-expression%2Fsoftware%2Fpipelines%2Flatest%2Fadvanced%2Freferences%23gtf
[5] filtered with : https://links.jianshu.com/go?to=https%3A%2F%2Fsupport.10xgenomics.com%2Fsingle-cell-gene-expression%2Fsoftware%2Fpipelines%2Flatest%2Fadvanced%2Freferences%23mkgtf
[6] Creating a Reference Package with cellranger mkref: https://links.jianshu.com/go?to=https%3A%2F%2Fsupport.10xgenomics.com%2Fsingle-cell-gene-expression%2Fsoftware%2Fpipelines%2Flatest%2Fadvanced%2Freferences%23header
[7] R包clusterProfiler的纯无参自定义物种注释的GO、KEGG富集分析及GSEA: https://links.jianshu.com/go?to=http%3A%2F%2Fblog.sciencenet.cn%2Fblog-3406804-1213409.html
[8] 关于人类参考基因组的一些认识: https://www.jianshu.com/p/3806afaf0c8c
[9] https://www.cnblogs.com/leezx/p/5710819.html
[10] Why Use Zebrafish to Study Human Diseases?: https://links.jianshu.com/go?to=https%3A%2F%2Firp.nih.gov%2Fblog%2Fpost%2F2016%2F08%2Fwhy-use-zebrafish-to-study-human-diseases