玩转参考基因组

在研究基因序列层面变化的时候，对基因组序列有一个全面的认知及学会怎么对序列基本操作十分重要。

参考基因组简介

1.1 参考基因组基本格式

参考基因组是以 .fa 结尾的 FASTA 格式文件。以 human 基因组为例。

在参考基因组中，每一个染色体第一行都是以 “>” 开头的描述序列的头部行（header line），数据库的不同，这里的信息可能略有不同。

• ***>***：表示这是FASTA格式的头部行。每个序列在FASTA文件中都以这种符号开头的行来描述其信息。
• ***1***：通常表示序列的标识符（ID）。在这种情况下，可能指的是染色体1。
• ***dna:chromosome***：指示序列的类型，这里表示该序列是一个染色体的DNA序列。
• ***chromosome:GRCh38:1:1:248956422:1***：这是一个详细的描述，包含多个信息：
  • ***chromosome***：指示这是一个染色体。
  • ***GRCh38***：表示使用的是人类基因组的GRCh38版本（Genome Reference Consortium Human Build 38）。
  • ***1***：指染色体1。
  • ***1:248956422***：表示序列的范围，从位置1到248,956,422，这是染色体1的全长。
  • ***1***：可能是指序列的方向或版本信息。
• ***REF***：通常表示这是参考序列（Reference Sequence），用于参考基因组的标准序列。

第一行以后就是这个染色体的序列了，每行有固定的字符数。可以从头部行看到这个染色体有248956422 个碱基。

1.2 为什么有那么多的 N，而不是常见的 A T C G 呢？

在参考基因组序列中，除了常见的碱基符号A、T、C、G（分别代表腺嘌呤、胸腺嘧啶、胞嘧啶和鸟嘌呤）之外，还可能出现其他符号。这些符号用于表示不确定性或变异。

N：表示未知碱基。通常用于未测序或不确定的区域。

$ grep -v "^>" Homo_sapiens.GRCh38.dna.chromosome.1.fa | tr -d '\n' | fold -w1 | sort | uniq -c
67070277 A
48055043 C
48111528 G
18475410 N
67244164 T

可以看到有 18475410 个碱基都是没有被测定的， 在1号染色体中占比 7.4% 左右。

1.3 染色体 1 序列最后一行是什么样的？

可以看到 1 号染色体序列的最后一行是两个 N字符，接着就是 10 号染色体。

因为是按照字符排序的，"1"、"10"、"11"、"2"这样的顺序，所以1 号染色体后就是10 号染色体。

参考基因组注释文件

2.1 注释文件基本格式

参考基因组注释文件提供了关于基因组序列的详细信息，包括基因的位置、功能、转录本、外显子、内含子等。常见的基因组注释文件格式包括GFF（General Feature Format）、GTF（Gene Transfer Format）和BED（Browser Extensible Data）格式。

这里以human gtf文件为例。

GTF文件由9个字段组成，每个字段用制表符分隔：

1. seqname：序列名称，通常是染色体编号（如*chr1、chrX*）或其他序列标识符。
2. source：注释的来源或生成数据的程序（如*HAVANA、ENSEMBL*）。
3. feature：特征类型（如*gene、transcript、exon、CDS*）。
4. start：特征的起始位置。
5. end：特征的结束位置。
6. score：特征的分数或评价值，通常是*.*表示无评分。
7. strand：链信息，*+表示正链，-*表示负链。
8. frame：阅读框信息，对于编码序列（CDS）有用，可能为*0、1、2或.*表示不适用。
   • 0：表示CDS的起始位置正好在密码子边界上。
   • 1：表示CDS的起始位置在密码子第一个碱基之后。
   • 2：表示CDS的起始位置在密码子第二个碱基之后。
   • **.**：对于非编码特征（如基因、转录本、外显子）或不适用的情况下使用。
9. attribute：包含键值对的附加信息，通常以分号分隔，例如*gene_id "GENE_ID"; transcript_id "TRANSCRIPT_ID";*。

2.2 注释文件中有哪些特征

可以看到第三列为基因组特征。

$ awk '{print $3}' *gtf | sort | uniq

CDS
Selenocysteine
exon
five_prime_utr
gene
start_codon
stop_codon
three_prime_utr
transcript

1. **CDS (Coding Sequence)**：
   • 表示编码序列，是基因的一部分，包含从起始密码子到终止密码子的核苷酸序列，翻译成蛋白质。
2. Selenocysteine：
   • 是一种特殊的氨基酸，由UGA密码子编码，通常作为终止密码子，但在特定情况下编码硒半胱氨酸。
3. Exon：
   • 外显子是基因的一部分，在转录后保留在成熟的mRNA中。外显子可能包含编码序列和非编码序列。
4. **Five_prime_utr (5' UTR)**：
   • 5'非翻译区是mRNA的起始部分，在编码序列之前。它不编码蛋白质，但在调节翻译中发挥重要作用。
5. Gene：
   • 是遗传信息的基本单位，包含编码序列和调控序列，负责控制特定性状或功能。
6. Start_codon：
   • 起始密码子是mRNA翻译起始的信号，通常是AUG，编码氨基酸甲硫氨酸。
7. Stop_codon：
   • 终止密码子是翻译终止的信号，通常为UAA、UAG或UGA，不编码任何氨基酸。
8. **Three_prime_utr (3' UTR)**：
   • 3'非翻译区是mRNA的末端部分，在编码序列之后。它在调节mRNA的稳定性和翻译效率中起作用。
9. Transcript：
   • 转录本是基因转录的产物，包括所有外显子和内含子，在剪接后形成成熟mRNA。

2.3 一个基因的完整信息应该怎么解读？

举例

chr1  example  gene          1000  5000  .  +  .  gene_id "gene1";
chr1  example  transcript    1000  5000  .  +  .  gene_id "gene1"; transcript_id "transcript1";
chr1  example  exon          1000  1200  .  +  .  gene_id "gene1"; transcript_id "transcript1";
chr1  example  exon          1500  1700  .  +  .  gene_id "gene1"; transcript_id "transcript1";
chr1  example  exon          2000  2400  .  +  .  gene_id "gene1"; transcript_id "transcript1";
chr1  example  five_prime_utr 1000 1099  .  +  .  gene_id "gene1"; transcript_id "transcript1";
chr1  example  CDS           1100  1200  .  +  0  gene_id "gene1"; transcript_id "transcript1";
chr1  example  CDS           1500  1700  .  +  2  gene_id "gene1"; transcript_id "transcript1";
chr1  example  CDS           2000  2300  .  +  0  gene_id "gene1"; transcript_id "transcript1";
chr1  example  three_prime_utr 2301 2400  .  +  .  gene_id "gene1"; transcript_id "transcript1";

• Gene:
  • chr1 example gene 1000 5000 . + . gene_id "gene1";
  • 表示基因的范围，从位置1000到5000。
• Transcript:
  • chr1 example transcript 1000 5000 . + . gene_id "gene1"; transcript_id "transcript1";
  • 表示转录本的范围，与基因范围相同。
• Exon:
  • 1000-1200
  • 1500-1700
  • 2000-2400
  • 三个外显子分别位于：
• 5' UTR:
  • chr1 example five_prime_utr 1000 1099 . + . gene_id "gene1"; transcript_id "transcript1";
  • 在转录本的起始部分，不被翻译。
• CDS:
  • 1100-1200
  • 1500-1700
  • 2000-2300
  • 编码序列分布在：
  • 这些区域将被翻译成蛋白质。
• 3' UTR:
  • chr1 example three_prime_utr 2301 2400 . + . gene_id "gene1"; transcript_id "transcript1";
  • 在转录本的末端部分，不被翻译。

2.4 如何在不同参考基因组版本之间转换注释坐标

可以明白，这种基于起始位置和终止位置的注释文件，如果参考基因组的版本更迭，将不再准确，最方便的方法是下载最新的版本，但如果你有特殊的需求，也可以进行“坐标转换”，有网页和工具能实现这个功能。

• liftOver
• CrossMap

参考基因组索引

3.1 参考基因组索引简介

参考基因组索引被称为基因组目录，为 FAI 格式文件，通常由 samtools faidx 命令生成。FAI 文件用于快速访问大型 FASTA 文件中的特定序列。它的每一行对应于 FASTA 文件中的一个序列，包含以下列信息：

1. 序列名称：对应于 FASTA 文件中序列的名称（即 > 后的部分）。
2. 序列长度：该序列的碱基总数。
3. 偏移量：该序列在 FASTA 文件中的起始字节位置。这是该序列在文件中的起始位置，以字节为单位。
4. 行长度：FASTA 文件中每一行的字符数（不包括换行符）。
5. 行总长度：FASTA 文件中每一行的字符总数（包括换行符）。

对于偏移量，这里举例说明：

假设有一个 FASTA 文件 *example.fasta*，内容如下：

>chr1
AGCTTAGCTA
GCTAGCTAGC
>chr2
CGTAGCTAGC
TAGCTAGCTA
CCT

使用 samtools faidx 创建的索引文件 fasta.fai 可能会是这样的：

chr1    20      6       10      11
chr2    23      34      10      11

• chr1 的偏移量（6）：
  • > 和 chr1 占用 5 个字节。
  • 第一个换行符占用 1 个字节。
  • 因此，*chr1* 的序列数据开始于第 6 个字节。
• chr2 的偏移量（34）：
  • chr1 的序列数据占用 20 个字节（两行，每行 10 个字符）。
  • 两个换行符（每行末尾一个）占用 2 个字节。
  • >chr2 占用 5 个字节。
  • 换行符占用 1 个字节。
  • 因此，*chr2* 的序列数据开始于第 34 个字节（6 + 20 + 2 + 6 = 34）。

快速了解自己的参考基因组

4.1 查看参考基因组大小

这里强烈推荐SeqKit工具。

$ seqkit stat Homo_sapiens.GRCh38.dna.primary_assembly.fa
file                                         format  type  num_seqs        sum_len  min_len       avg_len      max_len
Homo_sapiens.GRCh38.dna.primary_assembly.fa  FASTA   DNA        194  3,099,750,718      970  15,978,096.5  248,956,422

4.2 有多少个染色体

# Ensemble GRCh38 为例，统计染色体数。（结果为22+X+Y+MT=25
$ cat *.fa | grep ">*chromosome*" | wc -l
25

4.3 每个染色体的长度

可通过samtools faidx命令获取每个染色体的长度。

samtools faidx *.fa
cat *.fa.fai

4.4 参考基因组注释类型查看

$ cat *.gtf | grep -v "#" | cut -f3 | sort | uniq -c
 899302 CDS
    130 Selenocysteine
1668828 exon
 175449 five_prime_utr
  63140 gene
  98967 start_codon
  92935 stop_codon
 211674 three_prime_utr
 254129 transcript

从参考基因组中提取关注序列

5.1 使用samtools 工具

# 提取整个染色体序列
samtools faidx Homo_sapiens.GRCh38.dna.primary_assembly.fa 1 > chr1.fa

# 提取染色体1上1000-2000的位置
samtools faidx Homo_sapiens.GRCh38.dna.primary_assembly.fa 1:1000-2000 > sequence.fa

5.2 使用 bedtools 工具

这里我们先有一个bed 格式文件命名为 foo.bed

1       2581559 2584533 gene
1       2581559 2584533 transcript
1       2581559 2581650 exon
1       2583369 2583495 exon
1       2584124 2584533 exon

bedtools getfasta -fi Homo_sapiens.GRCh38.dna.primary_assembly.fa -bed foo.bed -fo output.fa

拓展 ：bio

bio 是一个非常优秀的软件，在获取基因的序列等方面极大的简化了步骤，这里演示一下获取fa数据的能力。

6.1 从 GeneBank 获取 fa 数据

$ bio fetch NC_045512 --format fasta  | head -3
>NC_045512.2 Severe acute respiratory syndrome coronavirus 2 isolate Wuhan-Hu-1, complete genome
ATTAAAGGTTTATACCTTCCCAGGTAACAAACCAACCAACTTTCGATCTCTTGTAGATCTGTTCTCTAAA
CGAACTTTAAAATCTGTGTGGCTGTCACTCGGCTGCATGCTTAGTGCACTCACGCAGTATAATTAATAAC

6.2 从GeneBank 获取 gtf 数据

$ bio fetch NC_045512 --format gff  | head -7
##gff-version 3
#!gff-spec-version 1.21
#!processor NCBI annotwriter
##sequence-region NC_045512.2 1 29903
##species .cgi?id=2697049
NC_045512.2     RefSeq  region  1       29903   .       +       .       ID=NC_045512.2:1..29903;Dbxref=taxon:2697049;collection-date=Dec-2019;country=China;gb-acronym=SARS-CoV-2;gbkey=Src;genome=genomic;isolate=Wuhan-Hu-1;mol_type=genomic RNA;nat-host=Homo sapiens;old-name=Wuhan seafood market pneumonia virus
NC_045512.2     RefSeq  five_prime_UTR  1       265     .       +       .       ID=id-NC_045512.2:1..265;gbkey=5'UTR

6.3 从 Ensembl 获取数据

# 获取cdna数据
bio fetch ENST00000288602 --type cdna

# 获取cds数据
bio fetch ENST00000288602 --type cds

# 获取蛋白数据
bio fetch ENST00000288602 --type protein