前言

已经有一个月没有更新文章，之前在知乎看到过一段话，大体是说写博客最重要的不是怎么把博客做的好看，而是如何持之以恒地产出，这句话让我感触颇深，正值国庆期间放假，就记录一下这段时间在分析数据时候发现的东西，以后万一有其他朋友发现并疑惑也好快速找到答案。

这里进行一个简要的概括：

  1. 长度为 1 并不代表该构成该基因的所有片段组合长度为 1，例如 1bp 的 exon 可能是微外显子，在剪切后与其他 exon 组合共同构成成熟 mRNA。
  2. CDS 长度不为 3 的整数倍，说明该基因的 CDS 注释不完善，例如缺失 3’ 或 5’ 端。
  3. 出于第 2 点现象，建议在下载注释时，选择 basic 版本，因为该版本中所有的注释都是完善的。

基因组注释中长度为 1 的现象解释

在使用基因组注释（版本：https://www.gencodegenes.org/human/release_19.html ）时，我发现了两个问题：

  1. 有些 CDS 的长度为 1。
  2. CDS 的长度不一定为 3 的倍数。

由于我是先发现的第一个问题，此后才发现的第二个问题，因此顺序如上，首先看一看第一个问题：

1
2
3
4
5
6
7
8
9
10
11
12
13
# First, download the GENCODE v19 annotation file
$ wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_19/gencode.v19.annotation.gtf.gz

# Next, extract and count the annotations with a length of 1bp
$ zcat gencode.v19.annotation.gtf.gz | awk '$4==$5' | cut -f 3 | sort | uniq -c | awk '$1 > 1'
  516 CDS
   19 exon
  249 start_codon
  377 stop_codon
 1058 UTR

# Optionally, view the full lines of annotations with a length of 1bp for closer examination
$ zcat gencode.v19.annotation.gtf.gz | awk '$4==$5' | less

能看到，不仅仅是 CDS，有许多其他的 annotation type 也仅有 1bp 的长度，甚至包括起始密码子和终止密码子，这个问题之后再讨论，首先说一说对于 CDS 而言为什么会有上面两个问题：

问题 ① 可能原因：

  • 长度为 1 有可能是微外显子（microexon），在剪切后它参与到蛋白翻译中。
  • 注释信息并不完全，见下文。

问题 ② 可能原因：

  • 对于一个基因，它可能具有多种剪切方式，对于每一种剪切方式，其涉及到的所有 CDS 的总长度为 3 的倍数，而一些外显子可能在多种剪切中出现，但在注释中它仅标明了一种。

    • 通过以下命令可以排除该点：

      1
      2
      3
      4
      5
      6
      7
      8
      9
      10
      11
      12
      $ zcat gencode.v19.annotation.gtf.gz | awk '$3=="CDS" && $1=="chr6"' | grep "99999716" | cut -f 1-5
      chr6    HAVANA  CDS     99999716        99999771
      chr6    HAVANA  CDS     99999716        99999771
      chr6    HAVANA  CDS     99999716        99999771
      chr6    HAVANA  CDS     99999716        99999771
      chr6    HAVANA  CDS     99999716        99999771
      chr6    HAVANA  CDS     99999716        99999771
      chr6    HAVANA  CDS     99999716        99999771
      chr6    HAVANA  CDS     99999716        99999771
      chr6    HAVANA  CDS     99999716        99999771
      chr6    HAVANA  CDS     99999716        99999771
      # 说明在不同剪切中重复出现的 CDS 会被记录

  • 注释信息并不完全，例如 tag "cds_start_NF" tag "cds_end_NF" 表明对于该转录本而言，起始密码子和终止密码子没有在注释的 CDS 区域中找到。

此外，另一些情况例如 alternative_5_UTR nonsense_mediated_decay 也可能有所影响。

如果你对具体的情况感兴趣，可以运行以下命令：

1
2
$ zcat gencode.v19.annotation.gtf.gz | awk '$3=="CDS"' | awk -F '\t' '{split($9, a, "gene_name \""); split(a[2], b, "\""); print b[1], ($5-$4+1)/3}' | awk '{ a[$1] += $2 } END { for (i in a) print i, a[i] }' | sort >> cdsdevide3.txt
$ zcat gencode.v19.annotation.gtf.gz | awk '$3=="CDS"' | awk -F '\t' '{split($9, a, "transcript_name \""); split(a[2], b, "\""); print b[1], ($5-$4+1)/3}' | awk '{ a[$1] += $2 } END { for (i in a) print i, a[i] }' >> transcriptdevide3.txt

例子：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 长度/3
A1BG-001 495
A1BG-004 305.667
A2ML1-001 1454
A2ML1-003 1004
A2ML1-004 963
A2ML1-006 181.333
A2ML1-008 81.3333
A2ML1-009 183.667

# 对应注释
transcript_name "A1BG-004"; exon_number 1; exon_id "ENSE00003185372.1"; level 2; protein_id "ENSP00000470909.1"; tag "mRNA_start_NF"; tag "mRNA_end_NF"; tag "cds_start_NF"; tag "cds_end_NF"; havana_gene "OTTHUMG00000183507.1";
transcript_name "A2ML1-008"; exon_number 1; exon_id "ENSE00002207222.1"; level 2; protein_id "ENSP00000440662.1"; tag "mRNA_end_NF"; tag "cds_end_NF"; havana_gene "OTTHUMG00000128499.7"; havana_transcript "OTTHUMT00000395958.1";
transcript_name "A2ML1-009"; exon_number 2; exon_id "ENSE00001798845.1"; level 2; protein_id "ENSP00000440057.1"; tag "mRNA_end_NF"; tag "cds_end_NF"; havana_gene "OTTHUMG00000128499.7"; havana_transcript "OTTHUMT00000395959.1";
transcript_name "A2ML1-006"; exon_number 6; exon_id "ENSE00001348621.1"; level 2; protein_id "ENSP00000445674.1"; tag "mRNA_start_NF"; tag "cds_start_NF"; havana_gene "OTTHUMG00000128499.7"; havana_transcript "OTTHUMT00000395963.1";

上述例子中，其他能被整除的皆具有 tag "basic"

讲完了 CDS，再来看看起始密码子和终止密码子是怎么回事。

当然，真相很简单，直接打印出所有的长度-数量表就能知道是为什么了：

1
2
3
4
$ zcat gencode.v19.annotation.gtf.gz | awk '($3=="start_codon" || $3=="stop_codon") {print ($5-$4+1)}' | sort | uniq -c
    626 1
    626 2
 159088 3

如果想要更详细地了解具体情况可以输入以下命令：

1
$ zcat gencode.v19.annotation.gtf.gz | awk '($3=="start_codon" || $3=="stop_codon") && ($5-$4+1)!=3' | less

在输入了上述命令以后可以看到，有些密码子是跨越了两个外显子的，也就是说，它们在基因组上的位置并不连续，但在转录和剪切后它们则共同构成终止密码子。

后记

事实上探索这些问题时，我突然意识到它的情况如何对我的分析并没有影响，但是我依然狠狠地将时间荒废在了这问题上。

或许有些时候重要的不是结果而是过程，探索这些问题的时候注释文件的结构和各种 shell 命令在我脑子里刻的都更深了些，这也算是一种进步。

这次探索也让我意识到基因注释是一个非常具有挑战性的操作，虽然用软件可以轻轻松松得到结果，但是这些结果也需要仔细的审查和验证才能知道它到底是不是真实可靠的。