如何快速从基因组中提取基因、转录本、蛋白、启动子、非编码序列？

NGS基础 - GTF/GFF文件格式解读和转换这篇文章有读者留言想要提取外显子，内含子，启动子，基因体，非编码区，编码区，TSS上游1500,TSS下游500的序列。下面我们就来示范如何提取这些序列。

NGS基础 - 参考基因组和基因注释文件提到了如何下载对应的基因组序列和基因注释文件。

假如我们已经拿到了基因组序列文件GRCh38.fa和基因注释文件GRCh38.gtf，也可从文后链接获取。

查看下文件内容和格式

基因组序列文件为FASTA格式，查看命令和内容如下（测试文件，只有1条染色体）：

# 查看前10行，每行查看前40个字符
# FASTA序列一般比较长，查看前面一部分字符是一个常用的方式
head GRCh38.fa | cut -c 1-40
>chr20
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

基因注释文件为GTF格式，只看前6列信息(第三列包含了不同的元件注释）

cut -f 1-6 GRCh38.gtf | head
chr20    ensembl_havana    gene    87250    97094    .
chr20    havana    transcript    87250    97094    .
chr20    havana    exon    87250    87359    .
chr20    havana    exon    96005    97094    .
chr20    ensembl_havana    transcript    87710    96533    .
chr20    ensembl_havana    exon    87710    87767    .
chr20    ensembl_havana    CDS    87710    87767    .
chr20    ensembl_havana    start_codon    87710    87712    .
chr20    ensembl_havana    exon    96005    96533    .
chr20    ensembl_havana    CDS    96005    96414    .

安装提取工具gffread

这里用到了gffread (https://github.com/gpertea/gffread)，安装方式如下 (若不理解，见这个为生信学习打造的开源Linux教程真香的软件安装部分)：

git clone https://github.com/gpertea/gffread
cd gffread
make release

提取转录本序列、CDS和蛋白序列

gffread -h可以参考所有可用参数，如果有特殊情况需要考虑的，还需配合其它参数使用。

1.获取转录本序列

gffread GRCh38.gtf -g GRCh38.fa -w GRCh38.transcripts.fa

内容如下：

head GRCh38.transcripts.fa
>ENST00000608838
ACAGGAATTCATATCGGGGTGATCACTCAGAAGAAAAGGTGAATACCGGATGTTGTAAGCTATTGAACTG
CCACAAGTGATATCTTTACACACCATTCTGCTGTCATTGGGTAGCTTTGAACCCCAAAAATGTTGGAAGA
ATAATGTAGGACATTGCAGAAGACGATGTTTAGATACTGAAAGGTACATACTTCTTTGTAGGAACAAGCT
ATCATGCTGCATTTCTATAATATCACATGAATATACTCGACGACCAGCATTTCCTGTGATTCACCTAGAG

2.获取CDS序列

# 获取CDS序列
gffread GRCh38.gtf -g GRCh38.fa -x GRCh38.cds.fa

内容如下

head GRCh38.cds.fa
>ENST00000382410
ATGAATATCCTGATGCTGACCTTCATTATCTGTGGGTTGCTAACTCGGGTGACCAAAGGTAGCTTTGAAC
CCCAAAAATGTTGGAAGAATAATGTAGGACATTGCAGAAGACGATGTTTAGATACTGAAAGGTACATACT
TCTTTGTAGGAACAAGCTATCATGCTGCATTTCTATAATATCACATGAATATACTCGACGACCAGCATTT
CCTGTGATTCACCTAGAGGATATAACATTGGATTATAGTGATGTGGACTCTTTTACTGGTTCCCCAGTAT
CTATGTTGAATGATCTGATAACATTTGACACAACTAAATTTGGAGAAACCATGACACCTGAGACCAATAC
TCCTGAGACTACTATGCCACCATCTGAGGCCACTACTCCCGAGACTACTATGCCACCATCTGAGACTGCT
ACTTCCGAGACTATGCCACCACCTTCTCAGACAGCTCTTACTCATAATTAA
>ENST00000382398
ATGAAGTCCCTACTGTTCACCCTTGCAGTTTTTATGCTCCTGGCCCAATTGGTCTCAGGTAATTGGTATG

3.获取蛋白序列

# 获取蛋白序列
gffread GRCh38.gtf -g GRCh38.fa -y GRCh38.protein.fa

内容如下

head GRCh38.protein.fa
>ENST00000382410
MNILMLTFIICGLLTRVTKGSFEPQKCWKNNVGHCRRRCLDTERYILLCRNKLSCCISIISHEYTRRPAF
PVIHLEDITLDYSDVDSFTGSPVSMLNDLITFDTTKFGETMTPETNTPETTMPPSEATTPETTMPPSETA
TSETMPPPSQTALTHN
>ENST00000382398
MKSLLFTLAVFMLLAQLVSGNWYVKKCLNDVGICKKKCKPEEMHVKNGWAMCGKQRDCCVPADRRANYPV
FCVQTKTTRISTVTATTATTTLMMTTASMSSMAPTPVSPTG
>ENST00000382388
MGLFMIIAILLFQKPTVTEQLKKCWNNYVQGHCRKICRVNEVPEALCENGRYCCLNIKELEACKKITKPP
RPKPATLALTLQDYVTIIENFPSLKTQST

解析GTF文件的结构

针对本GTF，对于gene元件，基因名字 (Gene symbol)在第14列。

head -n 1 GRCh38.gtf | sed 's/"/\t/g' | tr '\t' '\n' | sed = | sed 'N;s/\n/\t/'
1    chr20
2    ensembl_havana
3    gene
4    87250
5    97094
6    .
7    +
8    .
9    gene_id
10    ENSG00000178591
11    ; gene_version
12    6
13    ; gene_name
14    DEFB125
15    ; gene_source
16    ensembl_havana
17    ; gene_biotype
18    protein_coding
19    ;

针对本GTF，对于transcript元件，基因名字 (Gene symbol)在第18列。

sed -n '2p' GRCh38.gtf | sed 's/"/\t/g' | tr '\t' '\n' | sed = | sed 'N;s/\n/\t/'
1    chr20
2    havana
3    transcript
4    87250
5    97094
6    .
7    +
8    .
9    gene_id
10    ENSG00000178591
11    ; gene_version
12    6
13    ; transcript_id
14    ENST00000608838
15    ; transcript_version
16    1
17    ; gene_name
18    DEFB125
19    ; gene_source
20    ensembl_havana
21    ; gene_biotype
22    protein_coding
23    ; transcript_name
24    DEFB125-202
25    ; transcript_source
26    havana
27    ; transcript_biotype
28    processed_transcript
29    ; transcript_support_level
30    2
31    ;

这个查看信息在哪一列是很常用的检查文件结构提取对应信息的方式，简化为一个脚本checkCol.sh

检查某个文件的指定行（默认为第一行）

checkCol.sh -f GRCh38.gtf1    chr20
2    ensembl_havana
3    gene
4    87250
5    97094
6    .
7    +
8    .
9    gene_id "ENSG00000178591"; gene_version "6"; gene_name "DEFB125"; gene_source "ensembl_havana"; gene_biotype "protein_coding";

检查标准输入的第一行

sed 's/"/\t/g' GRCh38.gtf | checkCol.sh -f -
1    chr20
2    ensembl_havana
3    gene
4    87250
5    97094
6    .
7    +
8    .
9    gene_id
10    ENSG00000178591
11    ; gene_version
12    6
13    ; gene_name
14    DEFB125
15    ; gene_source
16    ensembl_havana
17    ; gene_biotype
18    protein_coding
19    ;

提取基因启动子序列

首先确定启动子区域，这里定义转录起始位点上游1000 bp和下游500 bp为启动子区域。

sed 's/"/\t/g' GRCh38.gtf | awk 'BEGIN{OFS=FS="\t"}{if($3=="gene") {if($7=="+") {start=$4-1000; end=$4+500;} else {if($7=="-") start=$5-500; end=$5+1000; } if(start<0) start=0; print $1,start,end,$14,$10,$7;}}' >GRCh38.promoter.bed

启动子区域如下（这个bed文件也可以用于ChIP-seq类型的数据分析确定peak是否在启动子区域）

head GRCh38.promoter.bed
chr20    86250    87750    DEFB125    ENSG00000178591    +
chr20    141369    142869    DEFB126    ENSG00000125788    +
chr20    156470    157970    DEFB127    ENSG00000088782    +
chr20    189181    190681    DEFB128    ENSG00000185982    -
chr20    226258    227758    DEFB129    ENSG00000125903    +
chr20    256736    258236    DEFB132    ENSG00000186458    +
chr20    266186    267686    AL034548.1    ENSG00000272874    +
chr20    290278    291778    C20orf96    ENSG00000196476    -
chr20    295968    297468    ZCCHC3    ENSG00000247315    +
chr20    347724    349224    NRSN2-AS1    ENSG00000225377    -

然后提取序列。这里用到了bedtools工具，官方有提供编译好的二进制文件，下载下来即可使用。

# -name: 输出基因名字（bed文件的第四列）
# -s: 考虑到正反链（对于启动子区域，是否考虑链的信息关系不太大）
bedtools getfasta -name -s -fi GRCh38.fa -bed GRCh38.promoter.bed >GRCh38.promoter.fa

序列信息如下：

head GRCh38.promoter.fa | cut -c 1-60
>DEFB125::chr20:86250-87750(+)
ATAATTTGAAGTGAGGTAATGTGATTCCTCTAGTTTTGTTCTTTTTGCTTAGGATGGCTT
>DEFB126::chr20:141369-142869(+)
AATATTCAAGAGAATGCCAAGAAAGCTACAAGAACAAATAGCAGGTCAGTCGTTGCCTGG
>DEFB127::chr20:156470-157970(+)
ATATCCGTCACCTCAAACATTTATCATTTGTATTGGGAACATTCAAAATCCTCTCTTCTA
>DEFB128::chr20:189181-190681(-)
AAAAAAGAAAAAGAACTCCAAGTCTAATAAGACCAGAGACCTGCCCTTTATGGGTCTGCA
>DEFB129::chr20:226258-227758(+)
GAGTGGAAGGTGGGAGGAGGGAGAGGATGAGGAAAAATAACTAATGGACACTAGGCTTAA

如果不想要坐标信息，可对序列名字做一下简化

cut -d ':' -f 1 GRCh38.promoter.fa >GRCh38.promoter.simplename.fa
head GRCh38.promoter.simplename.fa | cut -c 1-60
>DEFB125
ATAATTTGAAGTGAGGTAATGTGATTCCTCTAGTTTTGTTCTTTTTGCTTAGGATGGCTT
>DEFB126
AATATTCAAGAGAATGCCAAGAAAGCTACAAGAACAAATAGCAGGTCAGTCGTTGCCTGG
>DEFB127
ATATCCGTCACCTCAAACATTTATCATTTGTATTGGGAACATTCAAAATCCTCTCTTCTA
>DEFB128
AAAAAAGAAAAAGAACTCCAAGTCTAATAAGACCAGAGACCTGCCCTTTATGGGTCTGCA
>DEFB129
GAGTGGAAGGTGGGAGGAGGGAGAGGATGAGGAAAAATAACTAATGGACACTAGGCTTAA

提取基因序列

提取基因序列的操作也类似于提取启动子序列。这里要注意GFF文件的序列位置是从1开始，而bed文件的位置是从0开始，前闭后开，所以要对序列的起始位置进行-1的操作。

type="gene"
sed 's/"/\t/g' GRCh38.gtf | awk -v type="${type}" 'BEGIN{OFS=FS="\t"}{if($3==type) {print $1,$4-1,$5,$14,".",$7}}' >GRCh38.gene.bed
head GRCh38.gene.bed
chr20    87249    97094    DEFB125    .    +
chr20    142368    145751    DEFB126    .    +
chr20    157469    159163    DEFB127    .    +
chr20    187852    189681    DEFB128    .    -
chr20    227257    229886    DEFB129    .    +
chr20    257735    261096    DEFB132    .    +

提取基因序列

bedtools getfasta -name -s -fi GRCh38.fa -bed GRCh38.gene.bed >GRCh38.gene.fa
# 查看序列
head GRCh38.gene.fa | cut -c 1-60
>DEFB125::chr20:87249-97094(+)
ACAGGAATTCATATCGGGGTGATCACTCAGAAGAAAAGGTGAATACCGGATGTTGTAAGC
>DEFB126::chr20:142368-145751(+)
GCCATACACTTCAGCAGAGTTTGCAACTTCTCTTCTAAGTCTTTATCCTTCCCCCAAGGC
>DEFB127::chr20:157469-159163(+)
CTCTGAGGAAGGTAGCATAGTGTGCAGTTCACTGGACCAAAAGCTTTGGCTGCACCTCTT
>DEFB128::chr20:187852-189681(-)
GGCACACAGACCACTGGACAAAGTTCTGCTGCCTCTTTCTCTTGGGAAGTCTGTAAATAT

提取非编码RNA的序列

在GTF文件中有转录本类型的注释，包含下面这些注释类型

ntisense_RNA
lincRNA
miRNA
misc_RNA
processed_pseudogene
processed_transcript
protein_coding
rRNA
scaRNA
sense_intronic
sense_overlapping
snoRNA
snRNA
TEC
transcribed_processed_pseudogene
transcribed_unitary_pseudogene
transcribed_unprocessed_pseudogene
unitary_pseudogene
unprocessed_pseudogene

我们只筛选lincRNA

grep 'transcript_biotype "lincRNA"' GRCh38.gtf >GRCh38.lincRNA.gtf
gffread GRCh38.lincRNA.gtf -g GRCh38.fa -w GRCh38.lincRNA.fahead GRCh38.lincRNA.fa | cut -c 1-60
>ENST00000608495
GTCGCACGCGCTGGCCAAACGGGCGCACCAGACACTTTTCAGGGCCCTGCCAAAGACCTC
CTGGCGTCCCAGACACAAGAGATCCAGGCCAAGACTCACACTTCACAAGATACACAGACA
GGAACAGGAAATTCCATGAAACTTCCATTTACCCAATTAGCCGGACTCACTGAGCCCCAG
TCAACCAACTCCTACTAAAATTAAAAAGTAATGTGTGGTATAGATTGGAATAATAGACAT
AAACGATGGGAGGCGGAGAGGGGTGAGGGTTGAAAAATTACCTATTGGGTGCAACATTCA
AATGGGGCACTAGAAGCCCACTCCACCACTATGCAATATATGTATTTGTACCCCGTAAAT

提取一个个外显子序列

获取外显子的坐标

type="exon"
sed 's/"/\t/g' GRCh38.gtf | awk -v type="${type}" 'BEGIN{OFS=FS="\t"}{if($3==type) {print $1,$4-1,$5,$14,$20,$7}}' >GRCh38.exon.bed
# 查看文件内容
head GRCh38.exon.bed
chr20    87249    87359    ENST00000608838    DEFB125    +
chr20    96004    97094    ENST00000608838    DEFB125    +
chr20    87709    87767    ENST00000382410    DEFB125    +
chr20    96004    96533    ENST00000382410    DEFB125    +
chr20    142368    142686    ENST00000382398    DEFB126    +
chr20    145414    145751    ENST00000382398    DEFB126    +
chr20    142633    142686    ENST00000542572    DEFB126    +
chr20    145414    145488    ENST00000542572    DEFB126    +
chr20    145578    145749    ENST00000542572    DEFB126    +
chr20    157469    157593    ENST00000382388    DEFB127    +

提取序列

# -name: 输出基因名字（bed文件的第四列）
# -s: 考虑到正反链（对于启动子区域，是否考虑链的信息关系不太大）
bedtools getfasta -name -s -fi GRCh38.fa -bed GRCh38.exon.bed >GRCh38.exon.fa# 查看序列信息
head GRCh38.exon.fa | cut -c 1-60
>ENST00000608838::chr20:87249-87359(+)
ACAGGAATTCATATCGGGGTGATCACTCAGAAGAAAAGGTGAATACCGGATGTTGTAAGC
>ENST00000608838::chr20:96004-97094(+)
GTAGCTTTGAACCCCAAAAATGTTGGAAGAATAATGTAGGACATTGCAGAAGACGATGTT
>ENST00000382410::chr20:87709-87767(+)
ATGAATATCCTGATGCTGACCTTCATTATCTGTGGGTTGCTAACTCGGGTGACCAAAG
>ENST00000382410::chr20:96004-96533(+)
GTAGCTTTGAACCCCAAAAATGTTGGAAGAATAATGTAGGACATTGCAGAAGACGATGTT

提取一个个内含子序列

确定内含子区域

sed 's/"/\t/g' GRCh38.gtf | awk 'BEGIN{OFS=FS="\t";oldtr="";}{if($3=="exon") {tr=$14; if(oldtr!=tr) {start=$5; oldtr=tr;} else {print $1,start,$4-1,tr,$20,$7; start=$5;} } }' >GRCh38.intron.bed
# 查看文件内容
head GRCh38.intron.bed
chr20    87359    96004    ENST00000608838    DEFB125    +
chr20    87767    96004    ENST00000382410    DEFB125    +
chr20    142686    145414    ENST00000382398    DEFB126    +
chr20    142686    145414    ENST00000542572    DEFB126    +
chr20    145488    145578    ENST00000542572    DEFB126    +
chr20    157593    158773    ENST00000382388    DEFB127    +
chr20    189681    187852    ENST00000334391    DEFB128    -
chr20    227346    229277    ENST00000246105    DEFB129    +

提取序列同上。

往期精品(点击图片直达文字对应教程)

机器学习

后台回复“生信宝典福利第一波”或点击阅读原文获取教程合集