【佳學(xué)基因檢測】如何從基因組序列文件中獲取特定基因的全部序列、編碼序列、啟動子序列？
一、從基因組序列文件獲取特定基因序列需要參照基因組序列和注釋文件
1. 從NCBI數(shù)據(jù)庫下載人類基因組參照基因組數(shù)據(jù)文件。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行的內(nèi)容：
head -5 /media/jiaxue/0B8B16F90B8B16F9/reference/GCF_000001405.35_GRCh38.p9_genomic.fna
顯示結(jié)果為：
>NC_000001.11 Homo sapiens chromosome 1, GRCh38.p7 Primary Assembly
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
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基因注釋文件為
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行內(nèi)容：
head -15 /media/jiaxue/0B8B16F90B8B16F9/reference/GRCh38_latest_genomic.gff 顯示內(nèi)容為：
##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
顯示內(nèi)容經(jīng)過整理以說明不同的序列片段的注釋內(nèi)容的不同。
##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. 基因組注釋文件信息內(nèi)容的解釋
什么是GFF文件?
GFF格式是Sanger研究所賊先提出的，一種簡單的、方便的對于DNA、RNA以及蛋白質(zhì)序列的特征進(jìn)行描述的一種數(shù)據(jù)格式，比如基因序列的起點(diǎn)和終點(diǎn)坐標(biāo)。GFF格式是通過基因解碼技術(shù)中用來注釋基因序列的通用格式。
GFF文件包含了那些信息？
GFF文件由tab鍵隔開的9列組成，每一列代表不同的信息，下面是佳學(xué)基因?qū)Ω髁械恼f明：
先進(jìn)列：參考序列的編號，是chromosome or scaffold的編號；
第二列：基因信息注釋來源，一般為數(shù)據(jù)庫例或者注釋的機(jī)構(gòu)，如果未知，用“."代替；
第三列：基因信息的類型，如gene、mRNA、exon、CDS、UTR等；
第四列：第三列的基因信息在參考序列上的起始位置；
第五列：第三列的基因信息在參考序列上的終止位置；
第六列：注釋信息可信度得分，是注釋信息可能性的說明，可以是序列相似性比對時(shí)的E-values值或者基因預(yù)測時(shí)的P-values值，“.”表示為空；
第七列：該基因信息在基因序列的DNA鏈的標(biāo)識，是正鏈(+)還是負(fù)鏈(-)上；
第八列：當(dāng)基因信息是CDS時(shí)，表示起始編碼的位置，有效值為0、1、2，0表示該編碼框的先進(jìn)個(gè)密碼子的先進(jìn)個(gè)堿基位于其5'末端；1表示該編碼框的先進(jìn)個(gè)密碼子的先進(jìn)個(gè)堿基位于該編碼區(qū)外；2表示該編碼框的先進(jìn)個(gè)密碼子的先進(jìn)、二個(gè)堿基位于該編碼區(qū)外。
第九列：包含不同的注釋信息，用多個(gè)不同的名稱或者鍵值對來注釋。不同的注釋內(nèi)容之間以分號相隔，佳學(xué)基因?qū)ΤＲ娦畔⑦M(jìn)行一一解釋說明：
ID--注釋信息的編號，在一個(gè)GFF文件中必須少有；
Name--注釋信息的名稱，可以重復(fù)；
Alias--別名；
Parent--指明該基因信息所從屬的上一級ID。用于將exons聚集成transcript，將transripts聚集成gene；
Note--備注；
Dbxref--數(shù)據(jù)庫索引
二、參照基因組基因信息提取軟件介紹gffread
這里用到了gffread , 運(yùn)行如下命令，安裝gffread。
conda install -c bioconda gffread
運(yùn)行 gffread -h 查看軟件是否安裝成功。
提取轉(zhuǎn)錄本序列、CDS和蛋白序列
gffread -h可以參考所有可用參數(shù)，如果有特殊情況需要考慮的，還需配合其它參數(shù)使用。
1.獲取轉(zhuǎn)錄本序列
轉(zhuǎn)到注釋文件所在的文件夾： 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. 查看生成的轉(zhuǎn)錄本文件：
內(nèi)容如下：
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內(nèi)容如下
head -150 jiaxue.CDS.fa
>rna-NM_001005484.2
ATGAAGAAGGTAACTGCAGAGGCTATTTCCTGGAATGAATCAACGAGTGAAACGAATAACTCTATGGTGA
CTGAATTCATTTTTCTGGGTCTCTCTGATTCTCAGGAACTCCAGACCTTcctatttatgttgttttttgt
aTTCTATGGAGGAATCGTGTTTGGAAACCTTCTTATTGTCATAACAGTGGTATCTGACTCCCACCTTCAC
TCTCCCATGTACTTCCTGCTAGCCAACCTCTCACTCATTGATCTGTCTCTGTCTTCAGTCACAGCCCCCA
AGATGATTACTGACTTTTTCAGCCAGCGCAAAGTCATCTCTTTCAAGGGCTGCCTTGTTCagatatttct
ccttcacttctttgGTGGGAGTGAGATGGTGATCCTCATAGCCATGGGCTTTGACAGATATATAGCAATA
TGCAAGCCCCTACACTACACTACAATTATGTGTGGCAACGCATGTGTCGGCATTATGGCTGTCACATGGG
GAATTGGCTTTCTCCATTCGGTGAGCCAGTTGGCGTTTGCCGTGCACTTACTCTTCTGTGGTCCCAATGA
GGTCGATAGTTTTTATTGTGACCTTCCTAGGGTAATCAAACTTGCCTGTACAGATACCTACAGGCTAGAT
ATTATGGTCATTGCTAACAGTGGTGTGCTCACTGTGTGTTCTTTTGTTCTTCTAATCATCTCATACACTA
TCATCCTAATGACCATCCAGCATCGCCCTTTAGATAAGTCGTCCAAAGCTCTGTCCACTTTGACTGCTCA
CATTACAGTAGTTCTTTTGTTCTTTGGACCATGTGTCTTTATTTATGCCTGGCCATTCCCCATCAAGTCA
TTAGATAAATTCCTTGCTGTATTTTATTCTGTGATCACCCCTCTCTTGAACCCAATTATATACACACTGA
GGAACAAAGACATGAAGACGGCAATAAGACAGCTGAGAAAATGGGATGCACATTCTAGTGTAAAGTTTTA
G
>rna-XM_047436352.1
3.獲取蛋白序列
# 獲取蛋白序列
gffread GRCh38_latest_genomic.gff -g GRCh38_latest_genomic.fna -y jiaxue.protein.fa采用如下命令顯示內(nèi)容蛋白質(zhì)序列： head -150 jiaxue.protein.fa
>rna-NM_001005484.2MKKVTAEAISWNESTSETNNSMVTEFIFLGLSDSQELQTFLFMLFFVFYGGIVFGNLLIVITVVSDSHLHSPMYFLLANLSLIDLSLSSVTAPKMITDFFSQRKVISFKGCLVQIFLLHFFGGSEMVILIAMGFDRYIAICKPLHYTTIMCGNACVGIMAVTWGIGFLHSVSQLAFAVHLLFCGPNEVDSFYCDLPRVIKLACTDTYRLDIMVIANSGVLTVCSFVLLIISYTIILMTIQHRPLDKSSKALSTLTAHITVVLLFFGPCVFIYAWPFPIKSLDKFLAVFYSVITPLLNPIIYTLRNKDMKTAIRQLRKWDAHSSVKF>rna-XM_047436352.1解析GTF文件的結(jié)構(gòu)
針對本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 ;
這個(gè)查看信息在哪一列是很常用的檢查文件結(jié)構(gòu)提取對應(yīng)信息的方式，簡化為一個(gè)腳本checkCol.sh
檢查某個(gè)文件的指定行（默認(rèn)為先進(jìn)行）
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";
檢查標(biāo)準(zhǔn)輸入的先進(jìn)行
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 ;
提取基因啟動子序列
首先確定啟動子區(qū)域，這里定義轉(zhuǎn)錄起始位點(diǎn)上游1000 bp和下游500 bp為啟動子區(qū)域。
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
啟動子區(qū)域如下 （這個(gè)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工具，官方有提供編譯好的二進(jìn)制文件，下載下來即可使用。
# -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
如果不想要坐標(biāo)信息，可對序列名字做一下簡化
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開始，前閉后開，所以要對序列的起始位置進(jìn)行-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文件中有轉(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.fa
head GRCh38.lincRNA.fa | cut -c 1-60
>ENST00000608495
GTCGCACGCGCTGGCCAAACGGGCGCACCAGACACTTTTCAGGGCCCTGCCAAAGACCTC
CTGGCGTCCCAGACACAAGAGATCCAGGCCAAGACTCACACTTCACAAGATACACAGACA
GGAACAGGAAATTCCATGAAACTTCCATTTACCCAATTAGCCGGACTCACTGAGCCCCAG
TCAACCAACTCCTACTAAAATTAAAAAGTAATGTGTGGTATAGATTGGAATAATAGACAT
AAACGATGGGAGGCGGAGAGGGGTGAGGGTTGAAAAATTACCTATTGGGTGCAACATTCA
AATGGGGCACTAGAAGCCCACTCCACCACTATGCAATATATGTATTTGTACCCCGTAAAT
提取一個(gè)個(gè)外顯子序列
獲取外顯子的坐標(biāo)
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
# 查看文件內(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
提取一個(gè)個(gè)內(nèi)含子序列
確定內(nèi)含子區(qū)域
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
# 查看文件內(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 +
提取序列同上。
(責(zé)任編輯：佳學(xué)基因)
