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

2023-08-28 06:06| 来源: 网络整理| 查看: 265

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.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";

检查标准输入的第一行

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

【本文地址】

公司简介

联系我们

今日新闻

推荐新闻

专题文章