NGS基礎(chǔ) - GTF/GFF文件格式解讀和轉(zhuǎn)換這篇文章有讀者留言想要提取外顯子，內(nèi)含子，啟動子，基因體，非編碼區(qū)，編碼區(qū)，TSS上游1500,TSS下游500的序列。下面我們就來示范如何提取這些序列。

NGS基礎(chǔ) - 參考基因組和基因注釋文件提到了如何下載對應(yīng)的基因組序列和基因注釋文件。

假如我們已經(jīng)拿到了基因組序列文件GRCh38.fa和基因注釋文件GRCh38.gtf，也可從文后鏈接獲取。

查看下文件內(nèi)容和格式

基因組序列文件為FASTA格式，查看命令和內(nèi)容如下（測試文件，只有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)，安裝方式如下 (若不理解，見這個為生信學(xué)習(xí)打造的開源Linux教程真香的軟件安裝部分)：

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

提取轉(zhuǎn)錄本序列、CDS和蛋白序列

gffread -h可以參考所有可用參數(shù)，如果有特殊情況需要考慮的，還需配合其它參數(shù)使用。

1.獲取轉(zhuǎn)錄本序列

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

內(nèi)容如下：

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

2.獲取CDS序列

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

內(nèi)容如下

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

內(nèi)容如下

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

解析GTF文件的結(jié)構(gòu)

針對本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 ;

這個查看信息在哪一列是很常用的檢查文件結(jié)構(gòu)提取對應(yīng)信息的方式，簡化為一個腳本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 ;

提取基因啟動子序列

首先確定啟動子區(qū)域，這里定義轉(zhuǎn)錄起始位點上游1000 bp和下游500 bp為啟動子區(qū)域。

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

啟動子區(qū)域如下 （這個bed文件也可以用于ChIP-seq類型的數(shù)據(jù)分析確定peak是否在啟動子區(qū)域）

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: 考慮到正反鏈（對于啟動子區(qū)域，是否考慮鏈的信息關(guān)系不太大） 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文件中有轉(zhuǎn)錄本類型的注釋，包含下面這些注釋類型

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 # 查看文件內(nèi)容 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: 考慮到正反鏈（對于啟動子區(qū)域，是否考慮鏈的信息關(guān)系不太大） 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

提取一個個內(nèi)含子序列

確定內(nèi)含子區(qū)域

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 # 查看文件內(nèi)容 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 +

提取序列同上。

往期精品(點擊圖片直達文字對應(yīng)教程)

機器學(xué)習(xí)

后臺回復(fù)“生信寶典福利第一波”或點擊閱讀原文獲取教程合集

總結(jié)

以上是生活随笔為你收集整理的如何快速从基因组中提取基因、转录本、蛋白、启动子、非编码序列？的全部內(nèi)容，希望文章能夠幫你解決所遇到的問題。

如果覺得生活随笔網(wǎng)站內(nèi)容還不錯，歡迎將生活随笔推薦給好友。