一、bed文件介绍
bed文件是一种记录基因组不同(功能)区域在基因组上的位置以及其它注释信息的文本文件。它包含了由空格或者tab分隔的不同列,以记录不同的信息,每一行对应一个区域。它最早出现于人类基因组计划中,后被广泛应用。因为它不直接在基因组上进行标记和修改,在使用上更具效率。
bed文件最开始并没有一个标准的格式,因此 UCSC Genome Browser 对它的描述逐渐成为了大家的参考标准。它最少为3列,最多可为12列。bed文件辅助 UCSC Genome Browser 对不同片段进行可视化展示,因此第三行以后的信息多和如何展示这一区域有关,我就不在这里赘述了。以下是前三行的内容
chrom:染色体或者scaffold的名字;
chromStart:在染色体或者scaffold上的起始位置(包含起始位置对应的碱基对),如果是染色体的话,第一个碱基对的位置被标记为0;
chromEnd:在染色体或者scaffold上的终止位置(该终止位置的碱基对不被包含在内);
* scaffold:我们在二代测序中,获得的片段是reads,由reads组装成的为contigs,而contigs进一步组装,就能得到scaffold。
一个最简单的bed文件如下图所示:
因此,如果我们想要获得hg38参考基因组相对应的、包含外显子注释信息的参考基因组,最重要的,就是获得每一个染色体上外显子的起始位置和终止位置。
二、gtf文件介绍
相对应的gtf文件能够帮助我们获得想要的注释信息。gtf文件是由tab分隔的、包含不同列的文本文件。和bed文件相似,它同样包含了一定会有的列(9列),以及可有可没有的列(3列),定义如下:
seqname:序列名称;
source:这行信息的来源,产生这行信息的软件或者数据库;
feature:这行信息包含的特性类型,如基因、重复序列;
start:起始位置,第一个碱基对从1开始排序;
end:终止位置
score:评分
strand:正义链(+)反义链(-)
frame:为0,1,2三个数中的一个。如果为0,表明第一个密码子从第一个碱基开始。如果为1,表明有一个多出来的碱基,第一个密码子从第二个碱基开始。如果果为2,表明有两个多出来的碱基,第一个密码子从第三个碱基开始。
attribute:包含和这行信息有关的其他信息
一个简单的gtf文件示例如下图所示:
三、由gtf文件获得bed文件
因此,我们想要通过gtf文件获得含有外显子信息的bed文件,只需要找到第三列为exon的列,提取出它的第一列(序列名称)、第四列(起始位置)、第五列(终止位置),再把起始位置减一就可以了。下面是在Linux系统中的操作步骤。
其他作者用的是gencode(https://www.gencodegenes.org/human/)中下载的参考基因组,就用了gencode中提供的gtf文件。
我们首先要把gtf文件下载下来
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/gencode.v40.annotation.gtf.gz
接着,使用grep、awk、sort、bedtools等对gtf文件进行处理。为了节省空间,我们可以使用“|”将上一步处理得到的结果直接作为下一步的输入,而不进行保存。首先是解压gtf文件,并筛选出包含'transcript_type "protein_coding"'的行
gunzip -c gencode.v40.annotation.gtf.gz | \
grep 'transcript_type "protein_coding"' | \
解压后的gtf文件如下图所示
筛选出包含'transcript_type "protein_coding"'的行如下图所示
我们可以发现,包含'transcript_type "protein_coding"'的行里不仅有exon,还有star_codon等信息,而我们只需要exon的部分。接着,我们在这些行中找到第三列为exon的行,选取这些行里的第一列、第四列、第五列。然后,我们把第四列先变成数字,再减一。最后,我们在第一列和第四列后面加上分隔符“\t”(tab),在第五列后面加上分隔符“\n”以达到换行的目的。
awk '($3=="exon") {printf("%s\t%s\t%s\n",$1,int($4)-1,$5);}' | \
这一步处理后,我们得到的结果部分如下
最后我们再对最终的结果进行排序,先按照第一列进行排序(-k1,1)再按照第二列1以数值的方式进行排序(-k2,2n)。最后,我们再用bedtools把可能有重叠的部分整合之后,就能得到我们想要的结果。
sort -T . -t $'\t' -k1,1 -k2,2n | \
bedtools merge > hg38.bed
完整代码如下:
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_40/gencode.v40.annotation.gtf.gz
gunzip -c gencode.v40.annotation.gtf.gz | \
grep 'transcript_type "protein_coding"' | \
awk '($3=="exon") {printf("%s\t%s\t%s\n",$1,int($4)-1,$5);}' | \
sort -T . -t $'\t' -k1,1 -k2,2n | \
bedtools merge > hg38.bed
4、通过gtf文件获得bed12文件
###
在infer_experiment.py(RSeQC包)判断数据是否是链特异性
需要bed12文件,含有所有12列信息的bed文件
#安装gtfToGenePre
conda install -c bioconda ucsc-gtftogenepred
#准备好基因组gtf文件,从gtf转换为GenePred格式
gtfToGenePred -genePredExt -geneNameAsName2 genes.gtf gene.tmp
#从GenePred文件提取信息得到bed文件
awk '{print $2"\t"$4"\t"$5"\t"$1"\t0\t"$3"\t"$6"\t"$7"\t0\t"$8"\t"$9"\t"$10}' gene.tmp > genes_refseq.bed12