8. 参考基因组

1. 背景引入

  本小节开始讲述转录组测序的准备工作.因为做的是有参的基因组分析,所以首先是准备参考基因组、测序数据.当数据准备完成后,接下来是比对参考基因组,表达定量,合并成表达矩阵,差异表达分析.
  上面是转录组分析的大致步骤,这节我们介绍的是参考基因组.

2. 准备参考基因组

2.1 下载参考基因组

  我们准备参考基因组序列的基w本原则是下载物种最新、最优、最近的基因组.
  现有多个生物数据库存储他人的测序信息.但是有的数据库侧重不同的方向,比如GENCODE着重于人与小鼠的基因组.下图是GENCODE提供的人的基因组页面的fasta部分.我们解释一下\(GRCh38.p14\)都代表什么意思

  • GRC:The Genome Reference Consortium的缩写,即为基因组参考联盟.早期在发表文章的时候,发现参考基因组的序列不一致,会导致实验结果对不上,因此为了避免这种情况,GRC成立以统一参考基因组序列.
  • h38: 第38个版本.
  • p14: 第14版补丁.
      我们在测序人类的基因组做了很多工作,但仍有得到一些\(DNA\)片段,却不知道它属于哪个染色体的情况.\(Region\)的\(ALL\)表示包括这些不清楚情况的\(DNA\)片段.而\(pri\)则不包括这些,而是主要是染色体上的DNA.DNA不一定全在染色体上,叶绿体、线粒体上也有DNA.
    image
      GTF、GFF3文件是基因注释文件,下面两个是编码DNA的注释.PRI同样是在染色体上的DNA注释,而ALL是非染色体上的也包含.
    image
      Ensembl主要做基因组注释,并且侧重于脊椎动物方面.我们可以看一下小鼠的参考基因组序列.后面的13、14表示小鼠的13、14号染色体,nonchromosomal表示小鼠上未挂载在染色体上的DNA序列.primary_assembly表示主要挂载在染色体上的序列.toplevel也包含一些未挂载在染色体上的序列.rm.chromosome.2.fa.gz中的rm表示repeat mask,也就是将重复序列屏蔽掉.下面还有文件是将rm换成sm,详情可以看这个.

因此，在构建的目标物种 重复序列库中应排除掉这部分序列，即去除与已知物种基因相似性高的序列。 在获得重复 序列库后，可利用这部分序列将基因组中存在重复序列相似片段或区域“ 屏蔽”(mask) 。 所 谓屏蔽就是将原序列中的"A、T、C、G“ 用“ Ｎ”(hard mask)或小写的" ａ、t、ｃ、ｇ" (soft mask)表示，这样后续的基因预测软件将这部分序列按重复序列处理。对基因组中重复序列处理的 好坏将直接影响后续基因注释的质量。

image
image

  Ensembl同样提供基因组注释的下载,GTF是现在应用比较多的基因组注释文件,而GFF3是早期使用的文件格式(现大多不用).
image
  点进去后,文件目录如下所示,同样是一个只包含挂载在染色体(chr),一个包含所有的.
image
  Ensembl只做脊椎动物的基因组,但由于基因组注释做得很好,因此开创了网站做其他生物的基因组注释.
  \(NCBI\)即 National Center for Biotechnology Information,包含的基因组序列是最全的,但也是因为它的序列信息由用户上传,因此序列信息比较乱,如果有其他下载渠道不建议在此处下载.
  \(JGI\ PhytoZome\)是主要做植物的基因组网址.
  最后是物种数据库,很多物种在发表基因组的时候,会建立一个数据库.
image
  本小节要准备的基因组序列是柑橘的数据,对应就有柑橘的基因组数据库网址.
  打开后,点击Download进入下载页面,下面有对应的1.0、2.0,表示测序的基因组序列版本,这里我们下载最新的\(3.0\),有cDNA、genome(基因组)的文件,但我们需要下载的是基因组和基因组注释的文件(gff和genome).
image

2.2 基因组处理

  我们下载完成后,最后建立一个\(download.sh\)记录一下下载的指令,方便后面查阅下载了哪些文件.接下来我们可以先看一下gensome的文件,它的格式是第一行是染色体序号,第二行表示该染色体上的DNA序列,偶数行都很长,因此需要处理.
image
  我们按照fasta的规范处理.这里需要使用到seqtk包.
image
  安装好后,使用\(seqtk\ seq\ -l\)来划分行:
image

seqtk seq -l 60 raw/SWO.v3.0.genome.fa > genome.fa
# 注意要把60放在-l的邻近后面,因为给出的文档是-l INT说明后面接一个整数.

  我们处理好参考基因组文件后,接下来的步骤是对参考基因组构建索引文件,方便后续的的比对.这个操作有点类似于makeblastdb建立序列数据库,但这里我们使用的软件是seqkit.通过查阅文档可以知道构建index的指令是seqkit faidx

seqkit faidx genome.fa

之后会多出来文件genome.fa.fai.这个文件.这个文件与\(fasta\)文件一起,可以快速查看对应区域的序列.

  • 第一列 NAME : 序列的名称,只保留">"后,第一个空白之前的内容;
  • 第二列 LENGTH: 序列的长度,单位为bp;
  • 第三列 OFFSET: 第一个碱基的偏移量,从0开始计数，换行符也统计进行;
  • 第四列 LINEBASES : 行的碱基数,单位为bp；
  • 第五列 LINEWIDTH : 行宽,行的字节数,包括换行符,在windows系统中换行符为\r\n,要在序列长度的基础上加2;
    image

2.3 基因组注释处理

  之前我们说过,基因组注释文件现在大多是\(gtf\)格式,如果我们拿到\(gff\)格式文件,建议直接转为\(gtf\)格式.转换使用的软件是\(gffread\).

gffread -T raw/SWO.v3.0.gene.model.gff3 -o SWO_genes.gtf

  我们来解释一下\(gtf\)文件格式.gtf全称为gene transfer format,主要是用来对基因进行注释,当前所广泛使用的gtf格式为第二版(\(gtf2\)).以下均基于\(gtf2\)叙述.\(gtf\)同\(gff3\)很相似,也是9列内容,其内容如下:

  • seqname: 序列的名字.通常格式染色体ID或是contig ID.
  • source: 注释的来源.通常是预测软件名或是公共数据库.
  • feature: 基因结构.CDS,start_codon,stop_codon是一定要含有的类型.
  • start: 开始位点,从1开始计数.
  • end：结束位点.
  • score: 这一列的值表示对该类型存在性和其坐标的可信度,不是必须的,可以用点"."代替.
  • strand: 链的正向与负向,分别用加号+和减号-表示.
  • frame: 密码子偏移,可以是0、1或2.
  • attributes: 必须要有以下两个值.
    • gene_id value; 表示转录本在基因组上的基因座的唯一的ID。gene_id与value值用空格分开，如果值为空，则表示没有对应的基因。
    • transcript_id value; 预测的转录本的唯一ID。transcript_id与value值用空格分开，空表示没有转录本。

  我们转化完了,可以看一下文件.\(transcipt\)表示转录本,\(exon\)是外显子,\(cds\)是\(Coding\ sequence\)的缩写,是指编码一段蛋白产物的序列,是与蛋白质密码子一一对应的序列，注意其与mRNA序列的差异;ORF是open reading frame的缩写,翻译成开放阅读框,是指从一个起始密码子开始到一个终止密码子结束的一段序列,但并不是所有读码框都能表达出蛋白产物；CDS必定是一个ORF,但也可能包括多个ORF,相反，每个ORF不一定都是CDS.
  \(gtf\)文件包含比较底层的序列比如\(exon\)而不是像\(gff\)文件一层层罗列\(gene,mRNA\)等等.\(gtf\)比较便于程序处理.因此采用\(gtf\)更多.

image
  有时候我们获得的\(gtf\)文件没有\(transcript\)那一行,但是这样后续数据分析处理会不方便.所以我们要让\(gtf\)文件转为\(ensembl\)的\(gtf\)文件,这样就会带\(transcript\).

gtftk convert_ensembl -i temp.gtf > SWO_genes.gtf

  我们转换为gtf文件后,会发现\(ensembl\)的\(gtf\)带有\(gene\)那一行.这是\(ensembl\)的格式,也附有转录本的位置信息.
image

2.4 提取最长转录本

  转录本其实就是基因通过转录形成的一种或多种可供编码蛋白质的成熟的mRNA.一个基因有可能有多个转录本,原因是由于不同的剪接方式造成的,这也是我们之前说过的可变剪切.
  基因注释有不同水平.当我们对一个物种的基因组研究到了只知道外显子、内含子等序列位置,不知道外显子有多少种剪切形式,这就是注释到了基因水平.如果我们知道有多少种剪切形式,那就是注释到了转录本水平.
  当我们想提取最长转录本时,首先要确定基因组注释到了什么水平,一个很简单的方法就是确定\(gtf\)文件有多少\(gene\)以及多少\(transcript\),对比它们的数目,如果\(gene\)与\(transcript\)数目相同,就是注释到了基因水平.
image
  这里先用\(awk\)方法测试了一下,当然我们可以选择其他软件,比如\(gtftk\).
image
  先拉回正题,为什么我们要选择最长转录本.当相同基因得到不同转录本后,它们的功能是不一样的,这时我们要用哪一条转录本去代表这个基因呢?一种常见的做法就是使用最长转录本.这是我自己理解的,老师归纳的在这里:
image
  必须声明,现时期最佳选择其实是\(super\ transcript\),也就是在有多个剪切拼接时,对比不同的剪切,当到重合片段时,哪个长就选择哪个,这样一定是最佳的.以下面这张图举个例子.最佳\(transcipt\)实际是1(或3)、4、2.但提取\(supertranscript\)计算量大,操作麻烦,实际工作中还是使用最长转录本的多.
image

  确定好目的后,接下来我们需要将每一条基因的最长转录本提取出来.建议看文档理解.这里必须说明一下\(feature\)是什么,对于gtf文件来说,第三列的表头就是\(feature\),而值就是\(feature\)的值.也就是说,我们选择\(feature\)那一栏值为\(transcript\)的行.
image

# 提取最长转录本到文件中
gtftk short_long -l -i SWO_genes.gtf | gtftk select_by_key -k feature -v transcript | gtftk tabulate  -s "|" -k gene_id,transcript_id > longest_map.txt

  得到的文件形式如下:
image

  获取了每个基因对应的最长转录本后,我们安装基因的ID将每个最长转录本的蛋白质序列提取出来.我们查阅raw/SWO.v3.0.protein.fa的文件形式如下图,可以发现蛋白质序列正好对应着转录本ID.因此我们需要根据最长转录本的ID,将对应蛋白质序列提取出来.
image

sed '1d' longest_map.txt | cut -d '|' -f 2 | seqtk subseq raw/SWO.v3.0.protein.fa - > longest_transcript.proteins.fa
# 借此复习一下指令
# sed '1d' 是删除文件第一行
# cut  可以每行提取指定区域字符. -d 后面接分隔符 -f 每行取第2个区域 这里也可以使用awk
# seqtk subseq 根据输入的name.list提取in.fa里的序列
# - 据老师说,-可以理解为一个占位符,上个管道的输出结果给予了这个减号.相当于In.fa
# 最后得到最长转录本ID对应的蛋白质序列

  关于-占位符,看到一篇文章Click也是用了.这个到底啥时候能用???感觉是需要多个参数的时候可以使用.
  得到的文件格式如图所示:
image

  接下来就是将转录本ID置换成基因ID.从这里开始需要把前面的建立\(longest_map.txt\)改一下,改成下面的代码,不然替换不了基因ID(看官方的文档没看出啥来,可能是我水平太差吧)

#必须要让transcript_id在前,且两列用tab分割而不是其他字符.
gtftk short_long -l -i SWO_genes.gtf | gtftk select_by_key -k feature -v transcript | gtftk tabulate  -k transcript_id,gene_id > longest_map_modified.txt

  之后就能正常替换了.$表示匹配到尾巴.但我真的迷惑\(-p\)真的是接正则表达吗?用>(.+)$就换不了了.官方文档Go.

seqkit replace -K -p '(.+)$' -r '{kv}' -k longest_map_modified.txt longest_transcript.proteins.fa > proteins.fa
posted @ 2023-09-28 09:24  acmloser  阅读(368)  评论(0)  编辑  收藏  举报