macでインフォマティクス

macでインフォマティクス

NGS関連のインフォマティクス情報についてまとめています。

fastaのフォーマットを変換したり、指定サイズを取り出す seqmajic

 

Documentation

https://seqmagick.readthedocs.io/en/latest/

対応フォーマット

f:id:kazumaxneo:20181202093745p:plain

拡張子によってフォーマットが自動認識される。

 

インストール

mac os10.14の miniconda3-5.0環境でテストした。

依存

  • Python >= 3.4
  • biopython >= 1.70

本体 Github

#Anaconda環境でcondaを使い導入
conda install -c bioconda -y seqmagick

#python2環境では0.6系が最新
conda install -c bioconda -y seqmagick==0.6.1

#またはpipでも導入できる
pip install seqmagick

> seqmagick -h

$ seqmagick -h

usage: seqmagick [-h] [-V] [-v] [-q]

                 {help,convert,info,mogrify,primer-trim,quality-filter,extract-ids,backtrans-align}

                 ...

 

seqmagick - Manipulate sequence files.

 

positional arguments:

  {help,convert,info,mogrify,primer-trim,quality-filter,extract-ids,backtrans-align}

    help                Detailed help for actions using help <action>

    convert             Convert between sequence formats

    info                Info action

    mogrify             Modify sequence file(s) in place.

    primer-trim         Find a primer sequence in a gapped alignment, trim to

                        amplicon

    quality-filter      Filter reads based on quality scores

    extract-ids         Extract the sequence IDs from a file

    backtrans-align     Given a protein alignment and unaligned nucleotides,

                        align the nucleotides using the protein alignment.

                        Protein and nucleotide sequence files must contain the

                        same number of sequences, in the same order, with the

                        same IDs.

 

optional arguments:

  -h, --help            show this help message and exit

  -V, --version         Print the version number and exit

  -v, --verbose         Be more verbose. Specify -vv or -vvv for even more

  -q, --quiet           Suppress output

 

実行方法

1、convert: Biopythonサポートフォーマット間の変換(リンク

 

例えばfastaフォーマットをphyフォーマットに変換(アライメント実行後のfastaを使う)

seqmagick convert input.fa output.phy

stockholm format(wiki)に変換(アライメント実行後のfastaを使う)

seqmagick convert input.fa output.sto

multi-fastaから先頭5配列を出力

seqmagick convert --head 5 input.fa firs_five_seq.fa

multi-fastaからラスト5配列を出力。ただし、5配列のうち、100bp以下は出力されない(--head 5 => --min-length 100 )。

seqmagick convert --tail 5 --min-length 100 input.fa last_five_seq.fa

multi-fastaの100bp以上の配列を対象に、最後から5配列を出力(--min-length 100 => --head 5 )

seqmagick convert --min-length 100 --tail 5 input.fa last_five_seq.fa

multi-fastaの100bp以上の配列を対象に、それぞれの配列の1bpから5bpまでのポジションを抽出

seqmagick convert --min-length 100 --cut 1:5 input.fa pos1-5.fa

 DNA <=> RNA変換

#DNA => RNA
seqmagick convert --transcribe dna2rna input_dna.fa output_rna.fa 

#RNA => DNA
seqmagick convert --transcribe rna2dna input_rna.fa output_dna.fa
  • --transcribe {dna2rna,rna2dna}
  • --translate {dna2protein,rna2protein,dna2proteinstop,rna2proteinstop} 

 duplicationを除く(配列をチェックしている。ヘッダーが違っても配列が100%同一だと除かれる)

seqmagick convert --deduplicate-sequences input.fa output.fa

 

オプション一覧

$ seqmagick convert 

usage: seqmagick convert [-h] [--line-wrap N]

                         [--sort {length-asc,length-desc,name-asc,name-desc}]

                         [--apply-function /path/to/module.py:function_name[:parameter]]

                         [--cut start:end[,start2:end2]] [--relative-to ID]

                         [--drop start:end[,start2:end2]] [--dash-gap]

                         [--lower] [--mask start1:end1[,start2:end2]]

                         [--reverse] [--reverse-complement] [--squeeze]

                         [--squeeze-threshold PROP]

                         [--transcribe {dna2rna,rna2dna}]

                         [--translate {dna2protein,rna2protein,dna2proteinstop,rna2proteinstop}]

                         [--ungap] [--upper] [--deduplicate-sequences]

                         [--deduplicated-sequences-file FILE]

                         [--deduplicate-taxa] [--exclude-from-file FILE]

                         [--include-from-file FILE] [--head N]

                         [--max-length N] [--min-length N]

                         [--min-ungapped-length N] [--pattern-include REGEX]

                         [--pattern-exclude REGEX] [--prune-empty]

                         [--sample N] [--seq-pattern-include REGEX]

                         [--seq-pattern-exclude REGEX] [--tail N]

                         [--first-name] [--name-suffix SUFFIX]

                         [--name-prefix PREFIX]

                         [--pattern-replace search_pattern replace_pattern]

                         [--strip-range] [--input-format FORMAT]

                         [--output-format FORMAT]

                         [--alphabet {protein,dna,dna-ambiguous,rna,rna-ambiguous}]

                         source_file dest_file

ヘルプは"seqmagick convert --help"

 

2、mogrify: 修正コマンド(modifyは入力ファイルが上書きされるので注意)

ギャップを除く

seqmagick mogrify --ungap alignment.fa

 逆相補鎖にする

seqmagick mogrify --reverse-complement input.fasta

 複数の処理を同時実行

seqmagick mogrify --tail 5 --reverse-complement --ungap --deduplicate-sequences input.fasta

 

3、info: fastaやstoの要約統計

#カレントの全fastaとstoファイル
seqmagick info *.{fasta,sto}

$ seqmagick info *.{fasta,fna}

name                                             alignment    min_len   max_len   avg_len  num_seqs

adapter1.fasta                                   FALSE             41        41     41.00         1

adapter2.fasta                                   FALSE             47        47     47.00         1

alignment.fasta                                  TRUE            7608      7608   7608.00        59

escherichia_coli.fasta                           FALSE        5231428   52314285231428.00         1

ex2.fasta                                        FALSE           1182      2890   1491.07        76

rRNA.fasta                                       FALSE           2810      5183   3281.47        59

test.fasta                                       FALSE           2910      3018   2968.33         3

 

4、extract-ids: fastaのヘッダーID抽出(ヘッダー全てではない)

seqmagick extract-ids input.fasta > header

 

他にもクオリティフィルタリングを行うquality-filter(link)、アライメントファイルからプライマー部位をトリムするprimer-trim(link)などがあります。

引用

https://github.com/fhcrc/seqmagick