NGS基础 - GTF/GFF文件格式解读和转换这篇文章有读者留言想要提取外显子，内含子，启动子，基因体，非编码区，编码区，TSS上游1500,TSS下游500的序列。下面我们就来示范如何提取这些序列。
NGS基础 - 参考基因组和基因注释文件提到了如何下载对应的基因组序列和基因注释文件。
假如我们已经拿到了基因组序列文件GRCh38.fa和基因注释文件GRCh38.gtf，也可从文后链接获取。
查看下文件内容和格式
基因组序列文件为FASTA格式，查看命令和内容如下（测试文件，只有1条染色体）：
代码语言：javascript复制# 查看前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.获取转录本序列
代码语言：javascript复制gffread GRCh38.gtf -g GRCh38.fa -w GRCh38.transcripts.fa内容如下：
代码语言：javascript复制head GRCh38.transcripts.fa
>ENST00000608838
ACAGGAATTCATATCGGGGTGATCACTCAGAAGAAAAGGTGAATACCGGATGTTGTAAGCTATTGAACTG
CCACAAGTGATATCTTTACACACCATTCTGCTGTCATTGGGTAGCTTTGAACCCCAAAAATGTTGGAAGA
ATAATGTAGGACATTGCAGAAGACGATGTTTAGATACTGAAAGGTACATACTTCTTTGTAGGAACAAGCT
ATCATGCTGCATTTCTATAATATCACATGAATATACTCGACGACCAGCATTTCCTGTGATTCACCTAGAG2.获取CDS序列
代码语言：javascript复制# 获取CDS序列
gffread GRCh38.gtf -g GRCh38.fa -x GRCh38.cds.fa内容如下
代码语言：javascript复制head GRCh38.cds.fa
>ENST00000382410
ATGAATATCCTGATGCTGACCTTCATTATCTGTGGGTTGCTAACTCGGGTGACCAAAGGTAGCTTTGAAC
CCCAAAAATGTTGGAAGAATAATGTAGGACATTGCAGAAGACGATGTTTAGATACTGAAAGGTACATACT
TCTTTGTAGGAACAAGCTATCATGCTGCATTTCTATAATATCACATGAATATACTCGACGACCAGCATTT
CCTGTGATTCACCTAGAGGATATAACATTGGATTATAGTGATGTGGACTCTTTTACTGGTTCCCCAGTAT
CTATGTTGAATGATCTGATAACATTTGACACAACTAAATTTGGAGAAACCATGACACCTGAGACCAATAC
TCCTGAGACTACTATGCCACCATCTGAGGCCACTACTCCCGAGACTACTATGCCACCATCTGAGACTGCT
ACTTCCGAGACTATGCCACCACCTTCTCAGACAGCTCTTACTCATAATTAA
>ENST00000382398
ATGAAGTCCCTACTGTTCACCCTTGCAGTTTTTATGCTCCTGGCCCAATTGGTCTCAGGTAATTGGTATG3.获取蛋白序列
代码语言：javascript复制# 获取蛋白序列
gffread GRCh38.gtf -g GRCh38.fa -y GRCh38.protein.fa内容如下
代码语言：javascript复制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
检查某个文件的指定行（默认为第一行）
代码语言：javascript复制checkCol.sh -f GRCh38.gtf
1 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";检查标准输入的第一行
代码语言：javascript复制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是否在启动子区域）
代码语言：javascript复制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序列信息如下：
代码语言：javascript复制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如果不想要坐标信息，可对序列名字做一下简化
代码语言：javascript复制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 . 提取基因序列
代码语言：javascript复制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文件中有转录本类型的注释，包含下面这些注释类型
代码语言：javascript复制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
代码语言：javascript复制grep 'transcript_biotype "lincRNA"' GRCh38.gtf >GRCh38.lincRNA.gtf
gffread GRCh38.lincRNA.gtf -g GRCh38.fa -w GRCh38.lincRNA.fa
head GRCh38.lincRNA.fa | cut -c 1-60
>ENST00000608495
GTCGCACGCGCTGGCCAAACGGGCGCACCAGACACTTTTCAGGGCCCTGCCAAAGACCTC
CTGGCGTCCCAGACACAAGAGATCCAGGCCAAGACTCACACTTCACAAGATACACAGACA
GGAACAGGAAATTCCATGAAACTTCCATTTACCCAATTAGCCGGACTCACTGAGCCCCAG
TCAACCAACTCCTACTAAAATTAAAAAGTAATGTGTGGTATAGATTGGAATAATAGACAT
AAACGATGGGAGGCGGAGAGGGGTGAGGGTTGAAAAATTACCTATTGGGTGCAACATTCA
AATGGGGCACTAGAAGCCCACTCCACCACTATGCAATATATGTATTTGTACCCCGTAAAT提取一个个外显子序列
获取外显子的坐标
代码语言：javascript复制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 提取序列
代码语言：javascript复制# -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提取一个个内含子序列
确定内含子区域
代码语言：javascript复制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 提取序列同上。
