【佳学基因检测】如何从基因组序列文件中获取特定基因的全部序列、编码序列、启动子序列？

一、从基因组序列文件获取特定基因序列需要参照基因组序列和注释文件

1. 从NCBI数据库下载人类基因组参照基因组数据文件。https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_genomic.fna.gz，下载后的文件格式是FASTA文件格式。文件的存储为：GCF_000001405.35_GRCh38.p9_genomic.fna，查看前5行的内容：

head -5 /media/jiaxue/0B8B16F90B8B16F9/reference/GCF_000001405.35_GRCh38.p9_genomic.fna

显示结果为：

>NC_000001.11 Homo sapiens chromosome 1, GRCh38.p7 Primary Assembly

NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN

这里显示的是文件1号染色体的头文件、前4行序列文件。

2. 从NCBI下载人类基因组注释文件，下载地址为：https://ftp.ncbi.nlm.nih.gov/refseq/H_sapiens/annotation/GRCh38_latest/refseq_identifiers/GRCh38_latest_genomic.gff.gz，存储为：GRCh38_latest_genomic.gff.gz，将文件解压为GRCh38_latest_genomic.gff基因注释文件为GTF格式，共有9列，看前9列信息(第三列包含了不同的元件注释）

cut -f 1-9 /media/jiaxue/0B8B16F90B8B16F9/reference/GRCh38_latest_genomic.gff | head

##gff-version 3

#!gff-spec-version 1.21

#!processor NCBI annotwriter

#!genome-build GRCh38.p14

#!genome-build-accession NCBI_Assembly:GCF_000001405.40

#!annotation-date 03/15/2023

#!annotation-source NCBI RefSeq GCF_000001405.40-RS_2023_03

##sequence-region NC_000001.11 1 248956422

##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=9606

NC_000001.11 RefSeq region 1 248956422 . + . ID=NC_000001.11:1..248956422;Dbxref=taxon:9606;Name=1;chromosome=1;gbkey=Src;genome=chromosome;mol_type=genomic DNA

显示注释文件前15行内容：

head -15 /media/jiaxue/0B8B16F90B8B16F9/reference/GRCh38_latest_genomic.gff
显示内容为：

##gff-version 3

#!gff-spec-version 1.21

#!processor NCBI annotwriter

#!genome-build GRCh38.p14

#!genome-build-accession NCBI_Assembly:GCF_000001405.40

#!annotation-date 03/15/2023

#!annotation-source NCBI RefSeq GCF_000001405.40-RS_2023_03

##sequence-region NC_000001.11 1 248956422

##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=9606

NC_000001.11 RefSeq region 1 248956422 . + . ID=NC_000001.11:1..248956422;Dbxref=taxon:9606;Name=1;chromosome=1;gbkey=Src;genome=chromosome;mol_type=genomic DNA

NC_000001.11 BestRefSeq pseudogene 11874 14409 . + . ID=gene-DDX11L1;Dbxref=GeneID:100287102,HGNC:HGNC:37102;Name=DDX11L1;description=DEAD/H-box helicase 11 like 1 (pseudogene);gbkey=Gene;gene=DDX11L1;gene_biotype=transcribed_pseudogene;pseudo=true

NC_000001.11 BestRefSeq transcript 11874 14409 . + . ID=rna-NR_046018.2;Parent=gene-DDX11L1;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;Name=NR_046018.2;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2

NC_000001.11 BestRefSeq exon 11874 12227 . + . ID=exon-NR_046018.2-1;Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2

NC_000001.11 BestRefSeq exon 12613 12721 . + . ID=exon-NR_046018.2-2;Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2

NC_000001.11 BestRefSeq exon 13221 14409 . + . ID=exon-NR_046018.2-3;Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2

显示内容经过整理以说明不同的序列片段的注释内容的不同。

##gff-version 3

#!gff-spec-version 1.21

#!processor NCBI annotwriter

#!genome-build GRCh38.p14

#!genome-build-accession NCBI_Assembly:GCF_000001405.40

#!annotation-date 03/15/2023

#!annotation-source NCBI RefSeq GCF_000001405.40-RS_2023_03

##sequence-region NC_000001.11 1 248956422

##species https://www.ncbi.nlm.nih.gov/Taxonomy/Browser/wwwtax.cgi?id=9606

NC_000001.11 RefSeq region 1 248956422 . + . ID=NC_000001.11:1..248956422;Dbxref=taxon:9606; Name=1;chromosome=1;gbkey=Src;genome=chromosome;mol_type=genomic DNA

NC_000001.11 BestRefSeq pseudogene 11874 14409 . + . ID=gene-DDX11L1; Dbxref=GeneID:100287102,HGNC:HGNC:37102; Name=DDX11L1;description=DEAD/H-box helicase 11 like 1 (pseudogene);gbkey=Gene;gene=DDX11L1;gene_biotype=transcribed_pseudogene;pseudo=true

NC_000001.11 BestRefSeq transcript 11874 14409 . + . ID=rna-NR_046018.2; Parent=gene-DDX11L1;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;Name=NR_046018.2;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2

NC_000001.11 BestRefSeq exon 11874 12227 . + . ID=exon-NR_046018.2-1; Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2

NC_000001.11 BestRefSeq exon 12613 12721 . + . ID=exon-NR_046018.2-2; Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2

NC_000001.11 BestRefSeq exon 13221 14409 . + . ID=exon-NR_046018.2-3; Parent=rna-NR_046018.2;Dbxref=GeneID:100287102,GenBank:NR_046018.2,HGNC:HGNC:37102;gbkey=misc_RNA;gene=DDX11L1;product=DEAD/H-box helicase 11 like 1 (pseudogene);pseudo=true;transcript_id=NR_046018.2

3. 基因组注释文件信息内容的解释

什么是GFF文件?

GFF格式是Sanger研究所最先提出的，一种简单的、方便的对于DNA、RNA以及蛋白质序列的特征进行描述的一种数据格式，比如基因序列的起点和终点坐标。GFF格式是通过基因解码技术中用来注释基因序列的通用格式。

GFF文件包含了那些信息？

GFF文件由tab键隔开的9列组成，每一列代表不同的信息，下面是佳学基因对各列的说明：

第一列：参考序列的编号，是chromosome or scaffold的编号；

第二列：基因信息注释来源，一般为数据库例或者注释的机构，如果未知，用“."代替；

第三列：基因信息的类型，如gene、mRNA、exon、CDS、UTR等；

第四列：第三列的基因信息在参考序列上的起始位置；

第五列：第三列的基因信息在参考序列上的终止位置；

第六列：注释信息可信度得分，是注释信息可能性的说明，可以是序列相似性比对时的E-values值或者基因预测时的P-values值，“.”表示为空；

第七列：该基因信息在基因序列的DNA链的标识，是正链(+)还是负链(-)上；

第八列：当基因信息是CDS时，表示起始编码的位置，有效值为0、1、2，0表示该编码框的第一个密码子的第一个碱基位于其5'末端；1表示该编码框的第一个密码子的第一个碱基位于该编码区外；2表示该编码框的第一个密码子的第一、二个碱基位于该编码区外。

第九列：包含不同的注释信息，用多个不同的名称或者键值对来注释。不同的注释内容之间以分号相隔，佳学基因对常见信息进行一一解释说明：

ID--注释信息的编号，在一个GFF文件中必须唯一；

Name--注释信息的名称，可以重复；

Alias--别名；

Parent--指明该基因信息所从属的上一级ID。用于将exons聚集成transcript，将transripts聚集成gene；

Note--备注；

Dbxref--数据库索引

二、参照基因组基因信息提取软件介绍gffread

这里用到了gffread , 运行如下命令，安装gffread。
conda install -c bioconda gffread
运行 gffread -h 查看软件是否安装成功。
提取转录本序列、CDS和蛋白序列

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

1.获取转录本序列
转到注释文件所在的文件夹： cd /media/jiaxue/0B8B16F90B8B16F9/reference/

 gffread GRCh38_latest_genomic.gff -g GRCh38_latest_genomic.fna -w jiaxue.transcripts.fa 
输入基因组文件和注释文件需要匹配，否则会终止。输入匹配的文件后显示了如下记录：

FASTA index file GRCh38_latest_genomic.fna.fai created. 查看生成的转录本文件：

内容如下：

head GRCh38.transcripts.fa
>rna-NR_046018.2

CTTGCCGTCAGCCTTTTCTTTGACCTCTTCTTTCTGTTCATGTGTATTTGCTGTCTCTTAGCCCAGACTT

CCCGTGTCCTTTCCACCGGGCCTTTGAGAGGTCACAGGGTCTTGATGCTGTGGTCTTCATCTGCAGGTGT

CTGACTTCCAGCAACTGCTGGCCTGTGCCAGGGTGCAAGCTGAGCACTGGAGTGGAGTTTTCCTGTGGAG

AGGAGCCATGCCTAGAGTGGGATGGGCCATTGTTCATCTTCTGGCCCCTGTTGTCTGCATGTAACTTAAT

ACCACAACCAGGCATAGGGGAAAGATTGGAGGAAAGATGAGTGAGAGCATCAACTTCTCTCACAACCTAG

2.获取CDS序列

# 获取CDS序列

gffread GRCh38_latest_genomic.gff -g GRCh38_latest_genomic.fna -x jiaxue.CDS.fa

内容如下

head -150 jiaxue.CDS.fa
>rna-NM_001005484.2

ATGAAGAAGGTAACTGCAGAGGCTATTTCCTGGAATGAATCAACGAGTGAAACGAATAACTCTATGGTGA

CTGAATTCATTTTTCTGGGTCTCTCTGATTCTCAGGAACTCCAGACCTTcctatttatgttgttttttgt

aTTCTATGGAGGAATCGTGTTTGGAAACCTTCTTATTGTCATAACAGTGGTATCTGACTCCCACCTTCAC

TCTCCCATGTACTTCCTGCTAGCCAACCTCTCACTCATTGATCTGTCTCTGTCTTCAGTCACAGCCCCCA

AGATGATTACTGACTTTTTCAGCCAGCGCAAAGTCATCTCTTTCAAGGGCTGCCTTGTTCagatatttct

ccttcacttctttgGTGGGAGTGAGATGGTGATCCTCATAGCCATGGGCTTTGACAGATATATAGCAATA

TGCAAGCCCCTACACTACACTACAATTATGTGTGGCAACGCATGTGTCGGCATTATGGCTGTCACATGGG

GAATTGGCTTTCTCCATTCGGTGAGCCAGTTGGCGTTTGCCGTGCACTTACTCTTCTGTGGTCCCAATGA

GGTCGATAGTTTTTATTGTGACCTTCCTAGGGTAATCAAACTTGCCTGTACAGATACCTACAGGCTAGAT

ATTATGGTCATTGCTAACAGTGGTGTGCTCACTGTGTGTTCTTTTGTTCTTCTAATCATCTCATACACTA

TCATCCTAATGACCATCCAGCATCGCCCTTTAGATAAGTCGTCCAAAGCTCTGTCCACTTTGACTGCTCA

CATTACAGTAGTTCTTTTGTTCTTTGGACCATGTGTCTTTATTTATGCCTGGCCATTCCCCATCAAGTCA

TTAGATAAATTCCTTGCTGTATTTTATTCTGTGATCACCCCTCTCTTGAACCCAATTATATACACACTGA

GGAACAAAGACATGAAGACGGCAATAAGACAGCTGAGAAAATGGGATGCACATTCTAGTGTAAAGTTTTA

>rna-XM_047436352.1

3.获取蛋白序列

# 获取蛋白序列

gffread GRCh38_latest_genomic.gff -g GRCh38_latest_genomic.fna -y jiaxue.protein.fa

采用如下命令显示内容蛋白质序列： head -150 jiaxue.protein.fa

>rna-NM_001005484.2

MKKVTAEAISWNESTSETNNSMVTEFIFLGLSDSQELQTFLFMLFFVFYGGIVFGNLLIVITVVSDSHLH

SPMYFLLANLSLIDLSLSSVTAPKMITDFFSQRKVISFKGCLVQIFLLHFFGGSEMVILIAMGFDRYIAI

CKPLHYTTIMCGNACVGIMAVTWGIGFLHSVSQLAFAVHLLFCGPNEVDSFYCDLPRVIKLACTDTYRLD

IMVIANSGVLTVCSFVLLIISYTIILMTIQHRPLDKSSKALSTLTAHITVVLLFFGPCVFIYAWPFPIKS

LDKFLAVFYSVITPLLNPIIYTLRNKDMKTAIRQLRKWDAHSSVKF

>rna-XM_047436352.1

解析GTF文件的结构

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

head -n 1 GRCh38.gtf | sed 's/"/	/g' | tr '	' '
' | sed = | sed 'N;s/
/	/'
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/"/	/g' | tr '	' '
' | sed = | sed 'N;s/
/	/'
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/"/	/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/"/	/g' GRCh38.gtf | awk 'BEGIN{OFS=FS="	"}{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/"/	/g' GRCh38.gtf | awk -v type="${type}" 'BEGIN{OFS=FS="	"}{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.fa

head GRCh38.lincRNA.fa | cut -c 1-60
>ENST00000608495
GTCGCACGCGCTGGCCAAACGGGCGCACCAGACACTTTTCAGGGCCCTGCCAAAGACCTC
CTGGCGTCCCAGACACAAGAGATCCAGGCCAAGACTCACACTTCACAAGATACACAGACA
GGAACAGGAAATTCCATGAAACTTCCATTTACCCAATTAGCCGGACTCACTGAGCCCCAG
TCAACCAACTCCTACTAAAATTAAAAAGTAATGTGTGGTATAGATTGGAATAATAGACAT
AAACGATGGGAGGCGGAGAGGGGTGAGGGTTGAAAAATTACCTATTGGGTGCAACATTCA
AATGGGGCACTAGAAGCCCACTCCACCACTATGCAATATATGTATTTGTACCCCGTAAAT

提取一个个外显子序列

获取外显子的坐标

type="exon"
sed 's/"/	/g' GRCh38.gtf | awk -v type="${type}" 'BEGIN{OFS=FS="	"}{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/"/	/g' GRCh38.gtf | awk 'BEGIN{OFS=FS="	";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    +

提取序列同上。

(责任编辑：佳学基因)