提取出一个组装基因组的gap（N）和重复序列区域，保存为bed格式-优快云博客

本文介绍了一种使用Perl脚本从FASTA文件中提取基因组非序列位置（如N碱基的位置）的方法，并展示了如何将提取到的数据转换为标准BED格式文件，便于进一步的生物信息学分析。

参见：

Question: How to extract all non-seqenced positions from a genome (Fasta file)?

test.fa

>chr1
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNtaaattgttt
taaattgtttctgtttgcagttgacatgatctNNNNNatagaaaacacca
ataactctgccaaaaaatttagaattcataaatgaatttagtaaagttgc
>chr2
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNN
NNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNNtaaattgttt
taaattgtttctgtttgcagttgacatgatcttatatatagaaaacacca
ataactctgccaaaaaatttagaattcataaatgaatttagtaaagttgc

perl一行命令

perl -ne 'chomp;if( />(.*)/){$head = $1; $i=0; next};@a=split("",$_); foreach(@a){$i++;if($_ eq "N" && $s ==0 ){print "$head\t$i"; $s =1}elsif($s==1 && $_ ne "N"){print "\t$i\n";$s=0}}' test.fa

转为规范化的bed

cat gap.bed | awk 'BEGIN{i=0}{i++;print $1,$5,$6,"Gap"i}' > gap.2.bed

转载于:https://www.cnblogs.com/leezx/p/8647699.html