ゲノム配列からアノテーションされた遺伝子を抽出する

[備忘録]ゲノム解析　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