kuroの覚え書き

96の個人的覚え書き

seqkitでfastaファイルから一部を取り出す

multi fastaファイルから一部の遺伝子だけ取り出したサブfastaファイルを作るには

samtools faidx TAIR10_cDNA.fasta AT1G01010.1 AT1G01020.1 AT1G01030.1 AT1G01040.1 AT1G01050.1 > subset.fasta

のようにsamtoolsを使えばいいのだけれど、fastaファイルがcDNAで遺伝子のリストがtranscript IDではなくてgene IDだったとしたら、これでは抽出できない。
そんなときは

seqkit faidx TAIR10_cDNA.fasta -r AT1G01010* AT1G01020* AT1G01030* AT1G01040* AT1G01050* > subset.fasta

という感じでseqkitを使って正規表現で抽出すると良い。

更に遺伝子のリストがlist.txtなどテキストファイルに列記されているなら。

list=`cat list.txt`
seqkit faidx TAIR10_cDNA.fasta -r $list > subset.fasta

みたいにテキストを読み込んで渡してやる。

さらにさらに、遺伝子のリストが10000遺伝子とかめちゃめちゃ膨大だと引数の文字制限にぶち当たって

Argument list too long

とか文句言われるので、

list=`cat list.txt`
echo "$list" | xargs seqkit faidx TAIR10_cDNA.fa -r  >> subset.fasta

とxargsを使うと良いみたい。

>>ではなくて>でも内容は同じだが、何故か順番は変わってくる模様。

`はバックスラッシュなので注意。