NGS基础 - GTF/GFF文件格式解读和转换这篇文章有读者留言想要提取外显子，内含子，启动子，基因体，非编码区，编码区，TSS上游1500,TSS下游500得序列。下面我们就来示范如何提取这些序列。
NGS基础 - 参考基因组和基因注释文件提到了如何下载对应得基因组序列和基因注释文件。
假如我们已经拿到了基因组序列文件GRCh38.fa和基因注释文件GRCh38.gtf，也可从文后链接获取。
查看下文件内容和格式基因组序列文件为FASTA格式，查看命令和内容如下（测试文件，只有1条染色体）：
# 查看前10行，每行查看前40个字符# FASTA序列一般比较长，查看前面一部分字符是一个常用得方式head GRCh38.fa | cut -c 1-40>chr20NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
基因注释文件为GTF格式，只看前6列信息(第三列包含了不同得元件注释）
cut -f 1-6 GRCh38.gtf | headchr20 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 (github/gpertea/gffread)，安装方式如下 (若不理解，见这个为生信学习打造得开源Linux教程真香得软件安装部分)：
git clone github/gpertea/gffreadcd gffreadmake release提取转录本序列、CDS和蛋白序列
gffread -h可以参考所有可用参数，如果有特殊情况需要考虑得，还需配合其它参数使用。
1.获取转录本序列
gffread GRCh38.gtf -g GRCh38.fa -w GRCh38.transcripts.fa
内容如下：
head GRCh38.transcripts.fa>ENST00000608838ACAGGAATTCATATCGGGGTGATCACTCAGAAGAAAAGGTGAATACCGGATGTTGTAAGCTATTGAACTGCCACAAGTGATATCTTTACACACCATTCTGCTGTCATTGGGTAGCTTTGAACCCCAAAAATGTTGGAAGAATAATGTAGGACATTGCAGAAGACGATGTTTAGATACTGAAAGGTACATACTTCTTTGTAGGAACAAGCTATCATGCTGCATTTCTATAATATCACATGAATATACTCGACGACCAGCATTTCCTGTGATTCACCTAGAG
2.获取CDS序列
# 获取CDS序列gffread GRCh38.gtf -g GRCh38.fa -x GRCh38.cds.fa
内容如下
head GRCh38.cds.fa>ENST00000382410ATGAATATCCTGATGCTGACCTTCATTATCTGTGGGTTGCTAACTCGGGTGACCAAAGGTAGCTTTGAACCCCAAAAATGTTGGAAGAATAATGTAGGACATTGCAGAAGACGATGTTTAGATACTGAAAGGTACATACTTCTTTGTAGGAACAAGCTATCATGCTGCATTTCTATAATATCACATGAATATACTCGACGACCAGCATTTCCTGTGATTCACCTAGAGGATATAACATTGGATTATAGTGATGTGGACTCTTTTACTGGTTCCCCAGTATCTATGTTGAATGATCTGATAACATTTGACACAACTAAATTTGGAGAAACCATGACACCTGAGACCAATACTCCTGAGACTACTATGCCACCATCTGAGGCCACTACTCCCGAGACTACTATGCCACCATCTGAGACTGCTACTTCCGAGACTATGCCACCACCTTCTCAGACAGCTCTTACTCATAATTAA>ENST00000382398ATGAAGTCCCTACTGTTCACCCTTGCAGTTTTTATGCTCCTGGCCCAATTGGTCTCAGGTAATTGGTATG
3.获取蛋白序列
# 获取蛋白序列gffread GRCh38.gtf -g GRCh38.fa -y GRCh38.protein.fa
内容如下
head GRCh38.protein.fa>ENST00000382410MNILMLTFIICGLLTRVTKGSFEPQKCWKNNVGHCRRRCLDTERYILLCRNKLSCCISIISHEYTRRPAFPVIHLEDITLDYSDVDSFTGSPVSMLNDLITFDTTKFGETMTPETNTPETTMPPSEATTPETTMPPSETATSETMPPPSQTALTHN>ENST00000382398MKSLLFTLAVFMLLAQLVSGNWYVKKCLNDVGICKKKCKPEEMHVKNGWAMCGKQRDCCVPADRRANYPVFCVQTKTTRISTVTATTATTTLMMTTASMSSMAPTPVSPTG>ENST00000382388MGLFMIIAILLFQKPTVTEQLKKCWNNYVQGHCRKICRVNEVPEALCENGRYCCLNIKELEACKKITKPPRPKPATLALTLQDYVTIIENFPSLKTQST解析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 chr202 ensembl_havana3 gene4 872505 970946 .7 +8 .9 gene_id 10 ENSG0000017859111 ; gene_version 12 613 ; gene_name 14 DEFB12515 ; gene_source 16 ensembl_havana17 ; gene_biotype 18 protein_coding19 ;
针对本GTF，对于transcript元件，基因名字 (Gene symbol)在第18列。
sed -n '2p' GRCh38.gtf | sed 's/"/\t/g' | tr '\t' '\n' | sed = | sed 'N;s/\n/\t/'1 chr202 havana3 transcript4 872505 970946 .7 +8 .9 gene_id 10 ENSG0000017859111 ; gene_version 12 613 ; transcript_id 14 ENST0000060883815 ; transcript_version 16 117 ; gene_name 18 DEFB12519 ; gene_source 20 ensembl_havana21 ; gene_biotype 22 protein_coding23 ; transcript_name 24 DEFB125-20225 ; transcript_source 26 havana27 ; transcript_biotype 28 processed_transcript29 ; transcript_support_level 30 231 ;
这个查看信息在哪一列是很常用得检查文件结构提取对应信息得方式，简化为一个脚本checkCol.sh
检查某个文件得指定行（默认为第壹行）
checkCol.sh -f GRCh38.gtf1 chr202 ensembl_havana3 gene4 872505 970946 .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 chr202 ensembl_havana3 gene4 872505 970946 .7 +8 .9 gene_id 10 ENSG0000017859111 ; gene_version 12 613 ; gene_name 14 DEFB12515 ; gene_source 16 ensembl_havana17 ; gene_biotype 18 protein_coding19 ;提取基因启动子序列
首先确定启动子区域，这里定义转录起始位点上游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.bedchr20 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.fahead GRCh38.promoter.simplename.fa | cut -c 1-60>DEFB125ATAATTTGAAGTGAGGTAATGTGATTCCTCTAGTTTTGTTCTTTTTGCTTAGGATGGCTT>DEFB126AATATTCAAGAGAATGCCAAGAAAGCTACAAGAACAAATAGCAGGTCAGTCGTTGCCTGG>DEFB127ATATCCGTCACCTCAAACATTTATCATTTGTATTGGGAACATTCAAAATCCTCTCTTCTA>DEFB128AAAAAAGAAAAAGAACTCCAAGTCTAATAAGACCAGAGACCTGCCCTTTATGGGTCTGCA>DEFB129GAGTGGAAGGTGGGAGGAGGGAGAGGATGAGGAAAAATAACTAATGGACACTAGGCTTAA提取基因序列
提取基因序列得操作也类似于提取启动子序列。这里要注意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.bedhead GRCh38.gene.bedchr20 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_RNAlincRNAmiRNAmisc_RNAprocessed_pseudogeneprocessed_transcriptprotein_codingrRNAscaRNAsense_intronicsense_overlappingsnoRNAsnRNATECtranscribed_processed_pseudogenetranscribed_unitary_pseudogenetranscribed_unprocessed_pseudogeneunitary_pseudogeneunprocessed_pseudogene
我们只筛选lincRNA
grep 'transcript_biotype "lincRNA"' GRCh38.gtf >GRCh38.lincRNA.gtfgffread GRCh38.lincRNA.gtf -g GRCh38.fa -w GRCh38.lincRNA.fahead GRCh38.lincRNA.fa | cut -c 1-60>ENST00000608495GTCGCACGCGCTGGCCAAACGGGCGCACCAGACACTTTTCAGGGCCCTGCCAAAGACCTCCTGGCGTCCCAGACACAAGAGATCCAGGCCAAGACTCACACTTCACAAGATACACAGACAGGAACAGGAAATTCCATGAAACTTCCATTTACCCAATTAGCCGGACTCACTGAGCCCCAGTCAACCAACTCCTACTAAAATTAAAAAGTAATGTGTGGTATAGATTGGAATAATAGACATAAACGATGGGAGGCGGAGAGGGGTGAGGGTTGAAAAATTACCTATTGGGTGCAACATTCAAATGGGGCACTAGAAGCCCACTCCACCACTATGCAATATATGTATTTGTACCCCGTAAAT提取一个个外显子序列
获取外显子得坐标
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.bedchr20 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.bedchr20 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 +
提取序列同上。
