参考基因组与数据比对

约 1346 个字 7 行代码 预计阅读时间 5 分钟

参考基因组准备

FASTA与GFF/GTF文件格式

参考基因组包含两种文件: 基因组文件 fasta与注释文件 gff/gtf。

FASTA 格式

FASTA和FASTQ一样，都是一种记录序列的文本格式，前者同样可以记录氨基酸序列。

与FASTQ相比，FASTA的格式更加简单，其核心的序列表示块仅包含两部分：序列描述与序列内容，前者以">"开头，后者允许空格与空行的出现，直到下一个">"算作这一段序列的结束。

举例如下:

>MCHU - Calmodulin - Human, rabbit, bovine, rat, and chicken
MADQLTEEQIAEFKEAFSLFDKDGDGTITTKELGTVMRSLGQNPTEAELQDMINEVDADGNGTID
FPEFLTMMARKMKDTDSEEEIREAFRVFDKDGNGYISAAELRHVMTNLGEKLTDEEVDEMIREA
DIDGDGQVNYEEFVQMMTAK*

了解更多，看看wikipedia！

FASTA format: https://en.wikipedia.org/wiki/FASTA_format

gff/gtf

基因转移格式GTF (gene transfer format) 是一种用于保存基因结构信息的文件格式，它以固定的格式标记了基因组件在基因组上的坐标。

了解更多，看看Ensembl是怎么说的？

GFF/GTF File Format - Definition and supported options: https://asia.ensembl.org/info/website/upload/gff.html

FASTA文件负责回答“序列是什么”，GFF/GTF文件负责回答“基因组有什么组件，他们在哪里？”

GFF格式的总结如下表:

列数	关键字	描述
第1列	seqid	参考序列的id，通常是chromosome或者contig编号。该id的取名不能以'>'开头，不能包含空格。
第2列	source	注释的来源。如果未知，则用点（.）代替。一般指明产生此gff3文件的软件如havana或数据库如ensembl。
第3列	type	属性的类型。建议使用符合SO惯例的名称（sequence ontology），如gene, mRNA, start_codon, stop_codon, exon, CDS等。
第4列	start position	属性对应片段的起点，对应在第一列上的坐标。从1开始计数。
第5列	end position	属性对应片段的终点。一般比起点的数值要大。
第6列	score	得分，对于一些可以量化的属性，可以在此设置一个数值以表示程度的不同。如果为空，用点（.）代替。
第7列	strand	"+"表示正链，"-"表示负链，"."表示不需要指定正负链。
第8列	phase	步进。对于编码蛋白质的CDS来说，本列指定下一个密码子开始的位置。可以是0、1或2，表示到达下一个密码子需要跳过的碱基个数。对于其它属性，则用点（.）代替。
第9列	attributes	含众多属性feature的列表。格式为“标签=值”（tag=value）。不同属性之间用分号“;”隔开，可以存在空格。不过若有",=;"，则用URL（URL escaping rule）进行转义，同时TAB空格也需要转换为"%09"表示。

关于GFF格式第九列的详细描述如下:

ID	描述
ID	属性feature的唯一标识，一个GFF文件内ID具有唯一性
Name	属性feature所展示的名称
Alias	属性feature的第二个name，可以不具有唯一性
Parent	属性feature的上一级ID，可以将exons聚集成transcripts, transcripts聚集成genes, and so forth。一个feature可以有多个parents。
Target	表示比对的目标区域，格式为"target_id start end [strand]"，其中strand是可选（"+" or "-"）。如果target_id含有空格，必须转换为'%20'
Gap	比对结果的Gap信息
Note	文本描述
Ontology_term	A reference to an ontology term，对应的GO数据库ID
Is_circular	表明属性feature是否为环化

GTF文件格式总结如下:

数据库

常见三大基因组数据库: Ensembl、NCBI与UCSC。

三者链接如下:

Ensembel: https://asia.ensembl.org/

NCBI: https://www.ncbi.nlm.nih.gov/

UCSC: https://genome.ucsc.edu/

在Ensembel中可以下载到对应物种的基因组测序文件，如果是转录组文件则要注意从cdna目录下下载。

注释文件的release版本应当同测序文件一致。

对于不同数据库所提供的基因ID各有其特点，对于Ensembl来说，其基因ID格式为

ENS[species prefix][feature type prefix][a unique eleven digit number]

一般来说人类默认省略物种名。

有关于功能前缀Feature prefixes的说明:

功能前缀	含义
E	exon
FM	Ensembl protein family
G	gene
GT	gene tree
P	protein
R	regulatory feature（调控序列特征）
T	transcript

数据比对

数据比对的过程就好像是“拼图”，只是这个“拼图”是有重叠的那种，我们会将我们测序出来的读段比对到参考基因组上，这样就知道每个读段具体来自于哪个基因的那个片段上。

数据比对主要涉及两个工具: Hisat2和

Hisat2

Hisat2主要用于转录组数据比对，其运行前需要先构建对基因组的索引。

Hisat2 官方网站: https://daehwankimlab.github.io/hisat2/

Hisat2 官方仓库: https://github.com/DaehwanKimLab/hisat2

# 通过conda安装
❯ conda install bioconda::hisat2
# 构建索引
❯ hisat2-build <your_fasta_file> <prefix_for_index_file>
# 进行比对
# 输出文件名一般是.sam结尾
❯ hisat2 -x /path/to/your/index_file_prefix -1 /path/to/your/fastq/file -2 /path/to/your/fastq/file -p <thread_amount> -S /path/to/your/output/file

Hisat2常用参数如下:

参数	描述
-x	索引文件的前缀
-1	双末端测序结果的第一个文件。 若有多组数据，使用逗号将文件分隔。Reads的长度可以不一致。
-2	双末端测序结果的第二个文件。 若有多组数据，使用逗号将文件分隔，并且文件顺序要和-1参数的对应。Reads的长度可以不一致。
-U	单端数据文件。 若有多组数据，使用逗号将文件分隔。可以和-1,-2参数同时使用。Reads的长度可以不一致。
--rna-strandness	连特异性参数
-p	线程数

由Hisat2可以得到SAM文件，其可以通过samtools转化为BAM文件（体积更小）