根据坐标在基因组上面拿到碱基序列来设计引物

做DNA测序的朋友们一般来说，都会拿到突变位点信息，不管是SNV还是INDEL，都是一个基因组上面的坐标而已。而高通量测序的结果通常是需要做一下实验验证，最常见的就是sanger测序啦，需要设计引物来捕获一下突变位点附近的序列信息，查看是否该位点真的具有突变信息。一般来说，sanger测序能验证过的高通量测序的结果就不会受到审稿人的质疑啦！

如果仅仅是一两个位点，我们可以很容易通过各种各样的网页工具去查询到它的序列信息，但是高通量测序的结果往往是成千上万的，就算是节省成本，一般来说也会挑选100个左右的位点拿去设计引物进行sanger测序，一个个网页查询工作量有点大，这个时候就可以使用代码实现批量查询。

首先我们使用R语言模拟22个突变位点：

很简单的代码，这里我们在22条染色体上面各随机挑选一个位点哈，仅仅是作为程序的演示而已：

> pos=data.frame(chr=paste0('chr',1:22),start=sample(1:10000000,22))
> pos
     chr   start
1   chr1 2022626
2   chr2  696733
3   chr3 3250387
4   chr4 7673854
5   chr5 5408537
6   chr6 9719502
7   chr7 6581990
8   chr8 9601594
9   chr9 4787975
10 chr10 3528978
11 chr11 5885445
12 chr12 4356111
13 chr13 9586571
14 chr14 5893113
15 chr15 2299890
16 chr16 5854945
17 chr17 3117896
18 chr18 1789465
19 chr19 7853784
20 chr20 6409488
21 chr21 3040456
22 chr22 8896738

然后使用BSgenome::getSeq函数根据位点进行序列摘取

其中参考基因组序列来自于 BSgenome.Hsapiens.UCSC.hg38 包，这个包非常大，大家下载安装的时候一定要切换好镜像高速下载哦！

options(BioC_mirror="https://mirrors.tuna.tsinghua.edu.cn/bioconductor/")
options("repos" = c(CRAN="http://mirrors.cloud.tencent.com/CRAN/")) 
options(download.file.method = 'libcurl')
options(url.method='libcurl')
BiocManager::install("BSgenome.Hsapiens.UCSC.hg38",ask = F,update = F)

如果已经安装好对应的包，就可以直接使用BSgenome::getSeq函数，全部代码如下：

pos=data.frame(chr=paste0('chr',1:22),start=sample(1:10000000,22))
pos
library("BSgenome.Hsapiens.UCSC.hg38")
library("GenomicRanges")
# 真实情况下其实是读取你的突变坐标文件：
# pos=read.table('pos.txt')
# head(pos)
# 突变位点前后400bp供引物设计
pos1=GRanges(seqnames=pos[,1], ranges=IRanges(start=pos[,2]-400,end=pos[,2]))
pos2=GRanges(seqnames=pos[,1], ranges=IRanges(start=pos[,2],end=pos[,2]))
pos3=GRanges(seqnames=pos[,1], ranges=IRanges(start=pos[,2]+1,end=pos[,2]+401))
seq1 = BSgenome::getSeq(BSgenome.Hsapiens.UCSC.hg38, pos1)
seq2 = BSgenome::getSeq(BSgenome.Hsapiens.UCSC.hg38, pos2)
seq3 = BSgenome::getSeq(BSgenome.Hsapiens.UCSC.hg38, pos3)

输出可以是fasta文件或者txt文件，通常不会选择fasta文件，因为绝大部分没有生物信息学背景的生物学家其实不懂它。

# 
# names(seq) = paste0("SEQUENCE_", seq_along(seq)) 
# Biostrings::writeXStringSet(seq, "my.fasta")

tmp=cbind(as.data.frame(pos),
      as.data.frame(as.character(seq1)),
      as.data.frame(as.character(seq2)),
      as.data.frame(as.character(seq3)))
write.table(tmp,file = 'myFastq.txt',
            row.names = F,quote = F,col.names = F)

前面的代码里面，提取突变位点前后400bp供引物设计，不是很方便展现成为教程，所以我修改成为前后4bp，如下所示：

chr1 2022626 CCTCA A CTAG
chr2 696733 TCCCT T AGGT
chr3 3250387 CTACT T ACAC
chr4 7673854 CCACC C ACCC
chr5 5408537 GTAAA A ACTA
chr6 9719502 ATATT T AATT
chr7 6581990 TGGTT T GGCC
chr8 9601594 AACAC C CTGA
chr9 4787975 AAAGC C AAAC
chr10 3528978 TCATA A TCAC
chr11 5885445 AGATT T AATG
chr12 4356111 GTGGA A GAGC
chr13 9586571 NNNNN N NNNN
chr14 5893113 NNNNN N NNNN
chr15 2299890 NNNNN N NNNN
chr16 5854945 ATTGT T GGTT
chr17 3117896 TCAAA A CCCC
chr18 1789465 TTCTT T TACA
chr19 7853784 GGGAC C CGCC
chr20 6409488 CCAGG G GCTT
chr21 3040456 NNNNN N NNNN
chr22 8896738 NNNNN N NNNN

可以看到，每个突变位点上下游的4bp碱基序列都提取出来啦，就可以根据这些序列去设计引物做sanger测序验证。

bioconductor超级值得学习

假如我组建一个学习小组，关于这个bioconductor的，大家会加入吗？欢迎文末留言讨论：

http://bioconductor.riken.jp/packages/3.7/workflows/vignettes/sequencing/inst/doc/sequencing.html
http://bioconductor.org/packages/devel/workflows/vignettes/sequencing/inst/doc/sequencing.html
http://bioconductor.org/packages/3.3/bioc/vignettes/GenomicRanges/inst/doc/GenomicRangesHOWTOs.pdf

尤其是 http://www.bioconductor.org/packages/release/workflows/html/cytofWorkflow.html

当然了，看再多的教程，也需要活学活用！