hg38_intron_exton_bed文件

一、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