佳學基因遺傳病基因檢測機構排名,三甲醫(yī)院的選擇

基因檢測就找佳學基因!

熱門搜索
  • 癲癇
  • 精神分裂癥
  • 魚鱗病
  • 白癜風
  • 唇腭裂
  • 多指并指
  • 特發(fā)性震顫
  • 白化病
  • 色素失禁癥
  • 狐臭
  • 斜視
  • 視網膜色素變性
  • 脊髓小腦萎縮
  • 軟骨發(fā)育不全
  • 血友病

客服電話

4001601189

在線咨詢

CONSULTATION

一鍵分享

CLICK SHARING

返回頂部

BACK TO TOP

分享基因科技,實現人人健康!
×
查病因,阻遺傳,哪里干?佳學基因準確有效服務好! 靶向用藥怎么搞,佳學基因測基因,優(yōu)化療效 風險基因哪里測,佳學基因
當前位置:????致電4001601189! > 關于佳學 > 技術優(yōu)勢 >

【佳學基因檢測】如何從基因組序列文件中獲取特定基因的全部序列、編碼序列、啟動子序列?

假如我們已經拿到了基因組序列文件 GRCh38.fa 和基因注釋文件 GRCh38.gtf ,也可從文后鏈接獲取。 查看下文件內容和格式 基因組序列文件為FASTA格式,查看命令和內容如下(測試文件,只有1條染


佳學基因檢測】如何從基因組序列文件中獲取特定基因的全部序列、編碼序列、啟動子序列?


一、從基因組序列文件獲取特定基因序列需要參照基因組序列和注釋文件


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
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基因注釋文件為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表示該編碼框的先進個密碼子的先進個堿基位于該編碼區(qū)外;2表示該編碼框的先進個密碼子的先進、二個堿基位于該編碼區(qū)外。
 
第九列:包含不同的注釋信息,用多個不同的名稱或者鍵值對來注釋。不同的注釋內容之間以分號相隔,佳學基因對常見信息進行一一解釋說明:
 
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
G
>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    ;

提取基因啟動子序列

首先確定啟動子區(qū)域,這里定義轉錄起始位點上游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ū)域如下 (這個bed文件也可以用于ChIP-seq類型的數據分析確定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ū)域,是否考慮鏈的信息關系不太大)
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: 考慮到正反鏈(對于啟動子區(qū)域,是否考慮鏈的信息關系不太大)
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

提取一個個內含子序列

確定內含子區(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
# 查看文件內容
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    +

提取序列同上。

(責任編輯:佳學基因)
頂一下
(4)
100%
踩一下
(0)
0%
推薦內容:
來了,就說兩句!
請自覺遵守互聯網相關的政策法規(guī),嚴禁發(fā)布色情、暴力、反動的言論。
評價:
表情:
用戶名: 驗證碼: 點擊我更換圖片

Copyright © 2013-2033 網站由佳學基因醫(yī)學技術(北京)有限公司,湖北佳學基因醫(yī)學檢驗實驗室有限公司所有 京ICP備16057506號-1;鄂ICP備2021017120號-1

設計制作 基因解碼基因檢測信息技術部