在处理全基因组序列时，有时候需要一段一段的遍历基因组序列，分析各个区域序列的特征。比如我有一个细菌基因组5M，我需要截断成5000bp的片段，或者每隔1000bp截取5000bp，一直到把整个基因组截取完毕，然后对这些所有的片段做分析。
【通过awk字符串处理函数substr遍历字符串（如基因组）】当然最简单的是用python的字符串加索引功能，详情见本人文章使用countvectorizer 和tf-idf来编码文字/DNA序列 中的第一部分“1. 对dna 打断成kmer并处理成word格式”。这儿推荐一条更好的使用awk的命令，只需要一条命令，速度更快；
目标是将全基因组序列的fasta文件，每隔1000bp取5000bp的片段，直至取到末尾，将这些获得片段可以用于进一步的分析
$less 525362.3.fna
>NZ_GG693015Lactobacillus ruminis ATCC 25644 genomic scaffold SCAFFOLD1, whole genome shotgun sequence.[Lactobacillus ruminis ATCC 25644 | 525362.3]
tttaaaagttattttcgattcgatttaacaaaatcaatttgccgaccaaattattatcgcacctttcaaaagaaatgtcaacaataattttatttaattttgcgttgctcaattacaacaactttaatatcttaccacgttcaaacatatgcgtcaacatctttttgaaccaaatattttaaaagttgattgctcaatcagcgacctttttagtataccacgttaagaaaaaataacaaggcatttcttcaaaataatttatctacttgaatatttatttggaattctccttaattcctaaaacagttcgcctttttagcactgcattacaaaaaaaagaagccaagaatatcattctttagcttccttccaagttattttccctcagcgatatcctctaaggcacttttgacccatccgtcaattctgcttgcaatcggtgcaaaaccttcatgaacgtgacgattc...............# 使用awk 每隔1000bp打断一次
$ awk '{if(NR%2==1)a=$1;
else{len=length($0)-5000;
i=1;
while(i<=len){print a"_"i;
print (substr($0,i,5000));
i=i+6000;
}}}'525362.3.fna >525362.3_frag.fa$ grep ">" 525362.3_frag.fa
>NZ_GG693015_1
>NZ_GG693015_6001
>NZ_GG693015_12001
>NZ_GG693015_18001
>NZ_GG693015_24001
>NZ_GG693015_30001
>NZ_GG693015_36001
>NZ_GG693015_42001
>NZ_GG693015_48001
>NZ_GG693015_54001
>NZ_GG693015_60001
>NZ_GG693015_66001
>NZ_GG693015_72001
>NZ_GG693015_78001
>NZ_GG693015_84001
>NZ_GG693015_90001
可以看出已经全部打断完毕了
这个代码中
if(NR%2==1)a=$1记录奇数行，即序列抬头，记下序列的名字。对于偶数行（else命令段），首先获得序列的长度，由于要截取5000bp，因此序列起始位置应该不高于序列全长-5000，即len=length($0)-5000,随后使用while循环，第一次从1处取5000bp,使用awk内置函数substr（）字符串截取函数：(substr($0,i,5000))。第二次取得时候跳过1000bp，即从5000+1000=6000的位置取第二个5000bp，开始位置为i=i+6000=6001，即(substr($0,6001,5000))，从6001处取到第二个5000bp后，再次空出1000bp，再取5000，依次循环，即可取到所有序列。如果中间不想空1000bp，只需要把步长6000改成5000即可，即awk '{if(NR%2==1)a=$1;
else{len=length($0)-5000;
i=1;
while(i<=len){print a"_"i;
print (substr($0,i,5000));
i=i+5000;
}}}' genome.fa >genome_frag.fa
