多行fasta檔案分解成單個檔案

在使用kmer進行統計時，需要分別統計每條序列的kmer數目。如果所有樣本的fasta檔案均在乙個多行fasta檔案裡，如果把每一條序列提取出來？

有兩種方法，第一種方法先把序列id提取出來，然後採用grp + for迴圈的方法：

#獲得序列的id
grep
">" multiline.fa|
sed's/>//'
>fas.id
#分拆成單個檔案
for i in
`cat fas.id `
dogrep -a 1 -w $i multiline.fa >
$i.fa
done

當fasta序列行數不多時，這個**執行時間還可以接受。10萬條序列大約24小時可以處理完畢，1小時可以處理5000行左右。如果行數太大，時間就長，可以使用awk：

awk
'}' multiline.fa

速度可以實現1小時處理100萬條序列。

假如mutiline.fa為

$ cat mutiline.fa
>seq1
atcggggg
>seq2
gggctaaaaa

這個命令比較複雜，首先第乙個語句if($0~/^>/)a=$0，讀到標題行以">「開頭時（」^」）時，記錄標題在a中。其他鹼基行採用system函式執行shell中的echo命令。由於記錄的a中開頭含有">"，因此需要用轉義符號\消除原來的意義，但該轉義符號會導致」被轉義，因此需要再次被轉義，再次加乙個\，即為"echo \\"a。shell中即為命令echo \>seq1

這將列印》seq1,後面需要換行寫序列，因此需要跟乙個換行符號「\n」，由於引號「」在system命令中有意義，因此前後引號都需要轉義，即變成了\"\n\"，在awk 的system中shell命令兩邊需要用引號"「括起來，左右各加乙個引號」，即為"\"\n\"", 後面直接跟$0為序列，shell 中表示式為echo \>seq1"\n"tcgggg。這個地方用到兩層的轉義。

$ echo \>seq1"\n"attggggg
>seq1
attggggg

那麼下面就是把這個寫入檔案，檔案名字為seq2，即>seq2,正好是a，所以把a放在後面即可，所以在shell中的表示式即為echo \>seq1「\n」atcgggg>seq1`

$ echo \>seq2"\n"attggggg>seq1
$ ls
seq1

同理後面的seq2，seq3……也是類似方法

隨後對單個檔案計算kmer數目

~/bin/kmc3bx/kmc -k6 -ci1 -fa -cs2000 test.fa test_ind temp/

~/bin/kmc3bx/kmc_dump test_ind test_mer

多行fasta檔案分解成單個檔案

筆記 python對FASTA檔案的處理

使用python讀取和分析fasta檔案

fasta與fastq格式檔案解讀

多行fasta檔案分解成單個檔案

筆記 python對FASTA檔案的處理

使用python讀取和分析fasta檔案

fasta與fastq格式檔案解讀

相關推薦