前些天给学徒演示了猪狗的参考基因组构建索引 就顺便布置了作业，有意思的是她下载的时候，在两个参考基因组文件里面犹豫不决：
代码语言：javascript复制<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得到的突变位点，需要获取位点上下游序列这样可以去设计引物做一代测序，位点坐标如下：
代码语言：javascript复制chr17 43045748
chr17 43045761
chr17 43057069
chr17 43057116
chr17 43094853
简单脚本做出bed格式文件：
代码语言：javascript复制$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
代码语言：javascript复制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中会包含大量的变异信息，而这部分是很冗余并且一般也用不太到。
还有更多例子
所以说，bedtools 是能取代普通生信工程师的，更多实用小例子：
- Coverage analysis for targeted DNA capture. Thanks to Stephen Turner.
- Measuring similarity of DNase hypersensitivity among many cell types
- Extracting promoter sequences from a genome
- Comparing intersections among many genome interval files
- RNA-seq coverage analysis. Thanks to Erik Minikel.
- Identifying targeted regions that lack coverage. Thanks to Brent Pedersen.
- Calculating GC content for CCDS exons.
- Making a master table of ChromHMM tracks for multiple cell types.
