Documentation
https://seqmagick.readthedocs.io/en/latest/
対応フォーマット
拡張子によってフォーマットが自動認識される。
インストール
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
実行方法
１、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配列を出力。ただし、５配列のうち、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
３、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
４、extract-ids: fastaのヘッダーID抽出（ヘッダー全てではない）
seqmagick extract-ids input.fasta > header
他にもクオリティフィルタリングを行うquality-filter（link）、アライメントファイルからプライマー部位をトリムするprimer-trim（link）などがあります。
引用
https://github.com/fhcrc/seqmagick
