ゲノム配列からアノテーションされた遺伝子を抽出する
[備忘録]ゲノム解析 GFFREADの使い方
はじめに
こんにちは。
軟体動物のオミクス解析を行っている大学院生です。
今回はゲノムを扱う初学者が知っておくと便利なツール=GFFREADを紹介したいと思います。
背景
全生物のゲノム解読プロジェクトが進み、ゲノム解析の研究結果が増え続けています。そして、医学、農学、生物学などを中心に、あらゆる分野でゲノム配列を扱う場面が増えました。そんな中で私のようにゲノム解析初学者がゲノム配列をダウンロードして必要な配列を取得しなければいけない場面に遭遇することもあるかと思っています。そんな時に知っておくと心強いバイオインフォマティクスツールがGFFREADになります。
今回の概要
論文発表されたアノテーションファイル(GFF3 (GTF)ファイル)の情報を基に、ゲノムFASTAファイルから配列情報を抽出する。
ツール
・GFFREAD
必要なファイル
ゲノムFASTAファイルとGFF3 (GTF)ファイルの2つのデータを揃えるだけで良い。
*この2つのデータファイルは各論文内にアクセッション番号が開示されている。もしくはNCBIやEnsemblなどゲノムデータベースにアクセスするのも良い。
インストール
Conda環境で簡単にインストールすることができる。
$ conda install gffread
↓参考↓
https://bioconda.github.io/recipes/gffread/README.html
コマンド
$ gffread <input_gff> [-g <genomic_seqs_fasta> [-w <exons.fa>] [-x <cds.fa>] [-y <tr_cds.fa>]]
ex)↓エキソン領域の配列をs取得したい時↓
$ gffread genome.gff3 -w genome_exon_output.fa -g genome.fa
・-w exons.fa ：エキソンの配列を取得
・-x cds.fa ：cdsの塩基配列を取得
・-y tr_cds.fa ：cdsのアミノ酸配列を取得
抜きたいタグ(属性)がある場合
productの部分も同時に抜き出したい時は次のコマンドを追加で打つと良い。
ex) --table "product"
実用例 (ミトコンドリアゲノムから配列を入手)
ミトコンドリアゲノムから配列情報を抜き出してみる。例えば、tRNA配列の場合、ID情報だけではどの遺伝子か分かりにくいので、productも同時に抜き出すと良い(text:gff3の一部を参照)。
(例)gff3の一部
....(省略)....
NC_001276.1 RefSeq region 1 18224 . + . ID=NC_001276.1:1..18224;Dbxref=taxon:29159;Is_circular=true;Name=MT;gbkey=Src;genome=mitochondrion;mol_type=genomic DNA
NC_001276.1 RefSeq tRNA 736 803 . + . ID=rna-NC_001276.1:736..803;anticodon=(pos:768..770);gbkey=tRNA;product=tRNA-Ile
NC_001276.1 RefSeq exon 736 803 . + . ID=exon-NC_001276.1:736..803-1;Parent=rna-NC_001276.1:736..803;anticodon=(pos:768..770);gbkey=tRNA;product=tRNA-Ile
NC_001276.1 RefSeq tRNA 804 871 . + . ID=rna-NC_001276.1:804..871;anticodon=(pos:836..838);gbkey=tRNA;product=tRNA-Thr
NC_001276.1 RefSeq exon 804 871 . + . ID=exon-NC_001276.1:804..871-1;Parent=rna-NC_001276.1:804..871;anticodon=(pos:836..838);gbkey=tRNA;product=tRNA-Thr
NC_001276.1 RefSeq tRNA 892 960 . + . ID=rna-NC_001276.1:892..960;anticodon=(pos:923..925);gbkey=tRNA;product=tRNA-Glu
NC_001276.1 RefSeq exon 892 960 . + . ID=exon-NC_001276.1:892..960-1;Parent=rna-NC_001276.1:892..960;anticodon=(pos:923..925);gbkey=tRNA;product=tRNA-Glu
....(省略)....
そこで、次のコマンドを実行。
gffread sequence.gff3 -w sequence_output.fasta --table product -g sequence.fasta
すると次のfastaファイルが得られる。
sequence_output.fasta
>rna-NC_001276.1:736..803 tRNA-Ile
AGTATTGTGCCAGAGTTTTAATGGGCTTTGTTGATGTCAAAGAATACGAGAAATTTTCTCCAATACTT
>rna-NC_001276.1:804..871 tRNA-Thr
TCAGTTTTAGTATAAGTTTATTACGAGGGTCTTGTAAACCTTAGATCTGTATGTTTATGGGAACTGAA
>rna-NC_001276.1:892..960 tRNA-Glu
TCTGAGCTGGTGTAACTAAGCATGTAAGACTTTCAATCTTAAGGAGGGCAATCTTTCTCTCGTTCAGAT
....(省略)....
他
・GTFとGFFの変換 https://qiita.com/MaedaTaro_Umiushi/items/4b9d37d23614e53a3cdf
・目的の遺伝子の前後の配列を抜き出す
など、さまざまなオプションがある。
→→→ヘルプ参照
引用
ヘルプ
$ gffread -h
gffread v0.11.6. Usage:
gffread <input_gff> [-g <genomic_seqs_fasta> | <dir>][-s <seq_info.fsize>]
[-o <outfile>] [-t <trackname>] [-r [[<strand>]<chr>:]<start>..<end> [-R]]
[-CTVNJMKQAFPGUBHZWTOLE] [-w <exons.fa>] [-x <cds.fa>] [-y <tr_cds.fa>]
[-i <maxintron>] [--bed] [--table <attrlist>] [--sort-by <refseq_list.txt>]
Filter, convert or cluster GFF/GTF/BED records, extract the sequence of
transcripts (exon or CDS) and more.
By default (i.e. without -O) only transcripts are processed, discarding any
other non-transcript features. Default output is a simplified GFF3 with only
the basic attributes.
<input_gff> is a GFF file, use '-' for stdin
Options:
-i discard transcripts having an intron larger than <maxintron>
-l discard transcripts shorter than <minlen> bases
-r only show transcripts overlapping coordinate range <start>..<end>
(on chromosome/contig <chr>, strand <strand> if provided)
-R for -r option, discard all transcripts that are not fully
contained within the given range
-U discard single-exon transcripts
-C coding only: discard mRNAs that have no CDS features
--nc non-coding only: discard mRNAs that have CDS features
--ignore-locus : discard locus features and attributes found in the input
-A use the description field from <seq_info.fsize> and add it
as the value for a 'descr' attribute to the GFF record
-s <seq_info.fsize> is a tab-delimited file providing this info
for each of the mapped sequences:
<seq-name> <seq-length> <seq-description>
(useful for -A option with mRNA/EST/protein mappings)
Sorting: (by default, chromosomes are kept in the order they were found)
--sort-alpha : chromosomes (reference sequences) are sorted alphabetically
--sort-by : sort the reference sequences by the order in which their
names are given in the <refseq.lst> file
Misc options:
-F preserve all GFF attributes (for non-exon features)
--keep-exon-attrs : for -F option, do not attempt to reduce redundant
exon/CDS attributes
-G do not keep exon attributes, move them to the transcript feature
(for GFF3 output)
--keep-genes : in transcript-only mode (default), also preserve gene records
--keep-comments: for GFF3 input/output, try to preserve comments
-O process other non-transcript GFF records (by default non-transcript
records are ignored)
-V discard any mRNAs with CDS having in-frame stop codons (requires -g)
-H for -V option, check and adjust the starting CDS phase
if the original phase leads to a translation with an
in-frame stop codon
-B for -V option, single-exon transcripts are also checked on the
opposite strand (requires -g)
-P add transcript level GFF attributes about the coding status of each
transcript, including partialness or in-frame stop codons (requires -g)
--add-hasCDS : add a "hasCDS" attribute with value "true" for transcripts
that have CDS features
--adj-stop stop codon adjustment: enables -P and performs automatic
adjustment of the CDS stop coordinate if premature or downstream
-N discard multi-exon mRNAs that have any intron with a non-canonical
splice site consensus (i.e. not GT-AG, GC-AG or AT-AC)
-J discard any mRNAs that either lack initial START codon
or the terminal STOP codon, or have an in-frame stop codon
(i.e. only print mRNAs with a complete CDS)
--no-pseudo: filter out records matching the 'pseudo' keyword
--in-bed: input should be parsed as BED format (automatic if the input
filename ends with .bed*)
--in-tlf: input GFF-like one-line-per-transcript format without exon/CDS
features (see --tlf option below); automatic if the input
filename ends with .tlf)
Clustering:
-M/--merge : cluster the input transcripts into loci, discarding
"duplicated" transcripts (those with the same exact introns
and fully contained or equal boundaries)
-d <dupinfo> : for -M option, write duplication info to file <dupinfo>
--cluster-only: same as -M/--merge but without discarding any of the
"duplicate" transcripts, only create "locus" features
-K for -M option: also discard as redundant the shorter, fully contained
transcripts (intron chains matching a part of the container)
-Q for -M option, no longer require boundary containment when assessing
redundancy (can be combined with -K); only introns have to match for
multi-exon transcripts, and >=80% overlap for single-exon transcripts
-Y for -M option, enforce -Q but also discard overlapping single-exon
transcripts, even on the opposite strand (can be combined with -K)
Output options:
--force-exons: make sure that the lowest level GFF features are considered
"exon" features
--gene2exon: for single-line genes not parenting any transcripts, add an
exon feature spanning the entire gene (treat it as a transcript)
--t-adopt: try to find a parent gene overlapping/containing a transcript
that does not have any explicit gene Parent
-D decode url encoded characters within attributes
-Z merge very close exons into a single exon (when intron size<4)
-g full path to a multi-fasta file with the genomic sequences
for all input mappings, OR a directory with single-fasta files
(one per genomic sequence, with file names matching sequence names)
-w write a fasta file with spliced exons for each transcript
--w-add <N> for the -w option, extract additional <N> bases
both upstream and downstream of the transcript boundaries
-x write a fasta file with spliced CDS for each GFF transcript
-y write a protein fasta file with the translation of CDS for each record
-W for -w and -x options, write in the FASTA defline the exon
coordinates projected onto the spliced sequence;
for -y option, write transcript attributes in the FASTA defline
-S for -y option, use '*' instead of '.' as stop codon translation
-L Ensembl GTF to GFF3 conversion (implies -F; should be used with -m)
-m <chr_replace> is a name mapping table for converting reference
sequence names, having this 2-column format:
<original_ref_ID> <new_ref_ID>
-t use <trackname> in the 2nd column of each GFF/GTF output line
-o write the records into <outfile> instead of stdout
-T main output will be GTF instead of GFF3
--bed output records in BED format instead of default GFF3
--tlf output "transcript line format" which is like GFF
but exons, CDS features and related data are stored as GFF
attributes in the transcript feature line, like this:
exoncount=N;exons=<exons>;CDSphase=<N>;CDS=<CDScoords>
<exons> is a comma-delimited list of exon_start-exon_end coordinates;
<CDScoords> is CDS_start:CDS_end coordinates or a list like <exons>
--table output a simple tab delimited format instead of GFF, with columns
having the values of GFF attributes given in <attrlist>; special
pseudo-attributes (prefixed by @) are recognized:
@chr, @start, @end, @strand, @numexons, @exons, @cds, @covlen, @cdslen
-v,-E expose (warn about) duplicate transcript IDs and other potential
problems with the given GFF/GTF records