根据坐标在基因组上面拿到碱基序列来设计引物
做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
当然了,看再多的教程,也需要活学活用!
- 如何在Redhat中安装R的包及搭建R的私有源
- 什么是sparklyr
- 如何利用Dnsmasq构建小型集群的本地DNS服务器
- Cloudera Labs中的Phoenix
- 如何在CDH中使用Phoenix
- Java 8 时间 API 快速入门
- 如何在CDH中使用HPLSQL实现存储过程
- 如何掌握所有的编程语言
- 如何使用Sentry管理Hive外部表(补充)
- WebLogic XMLDecoder反序列化漏洞(CVE-2017-10271)漏洞复现&修复方案
- 如何在CDSW中使用R绘制直方图
- CTF学习交流群 第一期入群题writeup大放送
- 如何使用Hue创建Spark1和Spark2的Oozie工作流
- 【译】深入研究 Laravel 的依赖注入容器
- JavaScript 教程
- JavaScript 编辑工具
- JavaScript 与HTML
- JavaScript 与Java
- JavaScript 数据结构
- JavaScript 基本数据类型
- JavaScript 特殊数据类型
- JavaScript 运算符
- JavaScript typeof 运算符
- JavaScript 表达式
- JavaScript 类型转换
- JavaScript 基本语法
- JavaScript 注释
- Javascript 基本处理流程
- Javascript 选择结构
- Javascript if 语句
- Javascript if 语句的嵌套
- Javascript switch 语句
- Javascript 循环结构
- Javascript 循环结构实例
- Javascript 跳转语句
- Javascript 控制语句总结
- Javascript 函数介绍
- Javascript 函数的定义
- Javascript 函数调用
- Javascript 几种特殊的函数
- JavaScript 内置函数简介
- Javascript eval() 函数
- Javascript isFinite() 函数
- Javascript isNaN() 函数
- parseInt() 与 parseFloat()
- escape() 与 unescape()
- Javascript 字符串介绍
- Javascript length属性
- javascript 字符串函数
- Javascript 日期对象简介
- Javascript 日期对象用途
- Date 对象属性和方法
- Javascript 数组是什么
- Javascript 创建数组
- Javascript 数组赋值与取值
- Javascript 数组属性和方法