怎样用r程序分割fasta文件

前几天在柳城博客中看到《Perl处理Fasta序列的又一实例》一文。我以前就听说在生物医药领域perl使用极其广泛。我不懂生物学，是凑热闹的。

文章中给出了一段Perl脚本用来将一个fasta文件中的序列每两条一组分割成若干个文件。用Perl处理文本文件当然是非常合适的。我也经常处理文本文件，例如LabWindows程序采集的数据，不过一直使用Python，近期也没有学习Perl的计划，语言学得太多了反而容易把自己搞糊涂。Python非常好用，但是在进行一些简单的文本处理时不够方便，所以有时使用sed+awk作为补充。

对于上面提到的文章所述的例子，我试着用awk来手兆做（在Windows XP下用GnuWin32提毕激租供的awk）：

C:\LAB>awk "/^>/ &&i++%2==0{j++} {print >\铅滚"out_file\" j \".txt\"}" in_fasta.txt

似乎达到了要求。

1. 提取fasta文件abc.fas中序列>LG02的第164~202碱基之间序列，另存为abc_LG02_164_202.fas：

  $ more abc.fas|awk 'BEGIN{ORS=""}{if(/>/)print "\n"$0"@"else print $0}'|grep LG02|sed 's/@/\n/g'|grep -v '>'|awk '{print substr ($0,164+1,202-164-1)}' >abc_LG02_164_202.fas

在文件abc_LG02_164_202.fas第一行插入“>abc_LG02_164_202”

  $ sed  -i '1i >abc_LG02_164_202' abc_LG02_164_202.fas

2. 截取fasta文枣滑档件中序列XXXX及其前后n个碱基的序列：

  $ egrep  -o ‘.{n}XXXX.{n}’  abc.fas 

3. 从fasta文件中提取列表list中的序列（list格式：>a\n>b\n）：

  $ awk 'NR==FNR{a[$0]next} /^>/{b=($1 in a)} b'  list abc.fas >list.fas

4. 在fasta格式的氨基酸序列文件中，将序列中字符“.”去掉：

  $ more abc.fas|awk '{if(!/>/) gsub (/\./,"")print}' > abc.fasta

5. 去除空行,去除wei行

  $  sed '/^\s*$/d'

  $  sed -i '$d'  a.txt 

6. 批量重命名

  $  ls |grep fastq.gz|awk -F "." '{print$1}'|xargs -I {} mv {}.fastq.gz {}.fq.gz

7. 拆分文件，可将sh脚本拆分多个脚本同时后台运让芹行，避免运行任务过多内存不足

   $  凳乱split -l 10 abc.sh -d -a 4 abc_

8. 提取rsem分析结果的FPKM

ls|grep genes.results|awk '{print "awk \047{print $NF }\047 "$0 " >"$0".list"}'|sh

ls |grep results.list|awk 'BEGIN{ORS="\t"print "paste\tID"}{print$0}END{print" >Expression.FPKM\n"}'|sh

ls |grep results.list|sed 's/.genes.results.list//g'|awk 'BEGIN{ORS="\t"print "sed -i \0471i ID"}{print $0}END{print"\047 Expression.FPKM\n"}'|sh

more abc.list|awk '{print"grep \047\\b"$0"\\b\047 abc_Expression.FPKM >>abc.FPKM"}'|sh

---