cat *.fa合并fasta文件的坑
平常我們?cè)诤喜asta文件的時(shí)候翩腐，常常會(huì)使用到一句命令
cat *.fa > total.fa
cat *.fasta > total.fasta
其實(shí)原理就是將所有的.fa或者.fasta文件按順序輸出到標(biāo)準(zhǔn)輸出之中，之后使用>重定向到新的文件中余蟹。
但是在有的時(shí)候，會(huì)發(fā)生一些奇怪的事情诗充。
正常情況
比如我有三個(gè)fasta文件
1.fasta
>1
ATGACTGATGCTGCAGT
CGTAGCTAGCTCGATGA
ATGTC
2.fasta
>2
GTTGTACGTGCTAGTGC
GCTAGCTAGCTGTAGCT
ACAC
3.fasta
>3
CGTAGCTAGCTGACTCC
CGAGC
將這三個(gè)文件放到一個(gè)文件中溺健，注意要完全復(fù)制里面的內(nèi)容。
# 新建文件夾惰蜜，將新生成的fasta文件放在里面昂拂，避免循環(huán)往復(fù)的cat而導(dǎo)致錯(cuò)誤
mkdir temp
cat *.fasta > ./temp/total.fasta
然后查看序列名稱
cat ./temp/total.fasta | grep "^>"
得到結(jié)果
>1
>2
>3
沒(méi)有問(wèn)題。
異常
1.fasta
>1
ATGACTGATGCTGCAGT
CGTAGCTAGCTCGATGA
ATGTC
2.fasta
>2
GTTGTACGTGCTAGTGC
GCTAGCTAGCTGTAGCT
ACAC
3.fasta
>3
CGTAGCTAGCTGACTCC
CGAGC
同樣的還是這三個(gè)抛猖，將里面內(nèi)容復(fù)制之前的fasta文件中
mkdir temp
cat *.fasta > ./temp/total.fasta
cat ./temp/total.fasta | grep "^>"
結(jié)果
>1
這次>打頭的行只有一行呢格侯，其他兩個(gè)跑哪兒去了？這里不加^位置錨定標(biāo)志來(lái)找一找它們藏在哪兒
cat ./temp/total.fasta | grep ">"
輸出為
>1
ATGTC>2
ACAC>3
結(jié)果>2和>3竟然是跟在前一個(gè)fasta文件的屁股后面财著！
原因
原因就在于fasta文件的結(jié)尾沒(méi)有一行空行联四，也就是在文件末尾不是由空行結(jié)尾，而是序列作為結(jié)尾撑教，在序列的尾巴那兒沒(méi)有換行符朝墩，所以導(dǎo)致了前面的fasta文件的序列與后面的fasta文件的序列名的>連到了一塊。在數(shù)據(jù)庫(kù)中下載的fasta文件一般都是以空行結(jié)尾伟姐，但是我們自已寫(xiě)腳本或者手動(dòng)生成的fasta文件合并的時(shí)候就可能出現(xiàn)異常情況收苏，這種情況會(huì)導(dǎo)致在后續(xù)分析的時(shí)候報(bào)錯(cuò)。
如果說(shuō)報(bào)錯(cuò)還好愤兵，就怕那種不報(bào)錯(cuò)的鹿霸，在這里三條序列被當(dāng)作一條序列對(duì)待了，那生成的結(jié)果自然就有問(wèn)題恐似！
-
正常情況的合并
正常情況.png -
異常情況下的合并
異常情況.png
解決辦法
如果不確定fasta文件是否是以空行結(jié)尾杜跷，那么保險(xiǎn)起見(jiàn)就這樣做：
for i in $(ls *.fasta);
do
cat ${i}
# 在每一個(gè)fasta輸出之后再輸出一個(gè)空行
echo
done > ./temp/total.fasta
cat ./temp/total.fasta | grep ">"
結(jié)果正常了
>1
>2
>3
