【原】使用bedtools的getfasta功能来获取指定坐标上下游的序列

健明 2021-07-14

展开全文

前些天给学徒演示了猪狗的参考基因组构建索引就顺便布置了作业，有意思的是她下载的时候，在两个参考基因组文件里面犹豫不决：

<species>:   The systematic name of the species.
<assembly>:  The assembly build name.
<sequence type>:
 * 'dna' - unmasked genomic DNA sequences.
  * 'dna_rm' - masked genomic DNA.  Interspersed repeats and low
     complexity regions are detected with the RepeatMasker tool and masked
     by replacing repeats with 'N's.
  * 'dna_sm' - soft-masked genomic DNA. All repeats and low complexity regions

就是参考基因组是否带上rm后缀，我让她试一下，找一个fastq测序数据比对到两个参考基因组，结果她告诉我居然比对率很不一样，在rm后缀的参考基因组比对率比没有rm的参考基因组低10%。我仔细想了想，因为rm后缀的参考基因组意味着里面很多序列实际上是被NNNN占用了，所以一些在正常参考基因组里面比对成功的reads在rm后缀参考基因组比对失败很正常。

所以我让她提前了其中一个序列的比对坐标，然后去两个参考基因组里面看这个坐标里面的序列，是不是rm后缀的，被NNNN了。就发现她不会，所以提示了她getfasta可以根据BED/GFF/VCF文件提供的feature在染色体上的位置信息，从fasta中提取feature的碱基序列！

比如我想验证一些NGS得到的突变位点，需要获取位点上下游序列这样可以去设计引物做一代测序，位点坐标如下：

chr17    43045748
chr17    43045761
chr17    43057069
chr17    43057116
chr17    43094853

简单脚本做出bed格式文件：

$awk '{print $1,$2-300,$2+300}' tmp.pos
chr17 43045448 43046048
chr17 43045461 43046061
chr17 43056769 43057369
chr17 43056816 43057416
chr17 43094553 43095153
jianmingzeng 10:20:57 ~/tmp
# awk '{print $1"\t"$2-300"\t"$2+300}' tmp.pos > tmp.bed

标准用法是：bedtools getfasta -fi test.fa -bed test.bed -fo test.fa.out

bedtools getfasta -fi  ~/reference/genome/hg38/hg38.fa  -bed  tmp.bed -fo tmp.fa

写在最后，还有她下载的不是我曾下载过toplevel版本的基因组，比如人类的Homo_sapiens.GRCh38.dna.toplevel.fa.gz，文件大小1G，解压后54G！！！实际上用它对应的primary版本就够了：这个文件Homo_sapiens.GRCh38.dna.primary_assembly.fa.gz是正常的， primary的版本中是不包括haplotype info的，而top level中会包含大量的变异信息，而这部分是很冗余并且一般也用不太到。