流程

前置条件：

  • 在环境变量中可调用的 bedtools

本文参考：

Get intronic and intergenic sequences based on gff file from Biostars

https://www.biostars.org/p/112251/

本文可满足的需求：

  • 得到特定区域的 bed 文件（例如 exon / intron 等）。
  • 在得到区域信息的同时进行链信息的区分。

下载基因组注释文件（gtf）

以人类最新版本的 gencode 注释为例：

1
wget https://ftp.ebi.ac.uk/pub/databases/gencode/Gencode_human/release_46/gencode.v46.basic.annotation.gtf.gz

也可以选择已有的 gtf 文件进行。

得到 transcript 和 exon 区域信息

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
GTF="/path/to/gencode.v46.basic.annotation.gtf.gz" # 改为自己实际的注释文件路径，建议使用压缩格式以兼容以下命令
BASENAME=$(basename "$GTF" .gtf.gz)

TRANSCRIPT_BED="${BASENAME}.transcript.bed"
FORWARD_TRANS_BED="${BASENAME}.transcript.fors.bed"
BACKWARD_TRANS_BED="${BASENAME}.transcript.bacs.bed"
DOUBLE_TRANS_BED="${BASENAME}.transcript.doubletrans.bed"
EXON_BED="${BASENAME}.exon.bed"
INTRON_BED="${BASENAME}.intron.bed"
CDS_BED="${BASENAME}.cds.bed"

# 得到 transcript 区域
awk '$3=="transcript" {print $1, $4-1, $5, $7}' OFS='\t' <(zcat $GTF) > $TRANSCRIPT_BED

# 得到 exon 区域
awk '$3=="exon" {print $1, $4-1, $5, $7}' OFS='\t' <(zcat $GTF) > $EXON_BED

以上命令将得到 .transcript.bed 和 .exon.bed 结尾的文件，里面包含基因组上所有 transcript 和 exon 的信息。

不关注链信息时

1
2
3
4
5
6
7
8
9
10
# 合并排序 transcript 区域
bedtools sort -i $TRANSCRIPT_BED | bedtools merge -i - > $TRANSCRIPT_BED.tmp 
mv $TRANSCRIPT_BED.tmp $TRANSCRIPT_BED

# 合并排序 exon 区域
bedtools sort -i $EXON_BED | bedtools merge -i - > $EXON_BED.tmp 
mv $EXON_BED.tmp $EXON_BED

# 在 transcript 中排除 exon 得到 intron
bedtools subtract -a $TRANSCRIPT_BED -b $EXON_BED > $INTRON_BED

关注链信息时

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
# 得到正负链信息
awk '$4=="+"' $TRANSCRIPT_BED | bedtools sort -i - | bedtools merge -i - | awk 'OFS="\t"{print $0, "+"}' > $FORWARD_TRANS_BED
awk '$4=="-"' $TRANSCRIPT_BED | bedtools sort -i - | bedtools merge -i - | awk 'OFS="\t"{print $0, "-"}' > $BACKWARD_TRANS_BED

# 找出冲突区域（正负链上都存在 transcript）
bedtools intersect -a $FORWARD_TRANS_BED -b $BACKWARD_TRANS_BED > $DOUBLE_TRANS_BED

# 排除 transcript 冲突区域
bedtools subtract -a <(cat $FORWARD_TRANS_BED $BACKWARD_TRANS_BED | bedtools sort -i -) -b $DOUBLE_TRANS_BED > $TRANSCRIPT_BED

# 排除 exon 冲突区域
awk '$3=="exon" {print $1, $4-1, $5, $7}' OFS='\t' <(zcat $GTF) | bedtools sort -i - | bedtools merge -i - > $EXON_BED
bedtools intersect -a $TRANSCRIPT_BED -b $EXON_BED > $EXON_BED.tmp
mv $EXON_BED.tmp $EXON_BED

# 得到包含链信息的 intron
bedtools subtract -a $TRANSCRIPT_BED -b $EXON_BED > $INTRON_BED

相关的用途

通过上述命令，同理也可以获得例如 CDS 或 UTR 的 bed 文件。下游分析中，我们可以通过这些 bed 文件筛选在特定区域中的突变进行相关探索：

1
bedtools intersect -a input.vcf -b CDS.bed > input.cds.vcf

将 input.vcf 和 CDS.bed 请换为实际路径即可得到所有在 CDS 区域中的 variants。