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を使うと良いみたい。
>>ではなくて>でも内容は同じだが、何故か順番は変わってくる模様。
`はバックスラッシュなので注意。