搜索
查看: 1984|回复: 1

最方便的还是R语言呀,获取所有基因的上游100bp序列

[复制链接]

4

主题

51

帖子

482

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
482
发表于 2017-11-9 10:00:07 | 显示全部楼层 |阅读模式
不需要去找基因组序列,不需要去找基因结构坐标信息,也不需要写脚本了。

一切,都是在R里面完成,如下:
[AppleScript] 纯文本查看 复制代码
source("https://bioconductor.org/biocLite.R")
biocLite("BSgenome.Hsapiens.UCSC.hg19")
biocLite("TxDb.Hsapiens.UCSC.hg19.knownGene")

# BSgenome data package was made from the following source data files:
# chromFa.zip from [url]http://hgdownload.cse.ucsc.edu/goldenPath/hg19/bigZips/[/url]
library(BSgenome.Hsapiens.UCSC.hg19)
library(TxDb.Hsapiens.UCSC.hg19.knownGene)
BSgenome.Hsapiens.UCSC.hg19
genome <- BSgenome.Hsapiens.UCSC.hg19
seqlengths(genome)
genome$chr1  # same as genome[["chr1"]]

library(TxDb.Hsapiens.UCSC.hg19.knownGene)
txdb <- TxDb.Hsapiens.UCSC.hg19.knownGene
gn <- sort(genes(txdb))
up1000 <- flank(gn, width=1000)
up1000seqs <- getSeq(genome, up1000)



上一篇:Python小测007
下一篇:YouTube视频下载的23种方法大集合
回复

使用道具 举报

634

主题

1182

帖子

4030

积分

管理员

Rank: 9Rank: 9Rank: 9

积分
4030
发表于 2017-11-15 07:41:37 | 显示全部楼层
也有现成的轮子做这件事:Get a BED file of all regions not covered by the probes (+500 bp up/down)

bedtools complement -i probes.500bp.bed -g hg18.genome > probes.500bp.complement.bed
你这个问题很复杂,需要打赏,请点击 http://www.bio-info-trainee.com/donate 进行打赏,谢谢
回复 支持 反对

使用道具 举报

您需要登录后才可以回帖 登录 | 立即注册

本版积分规则

QQ|手机版|小黑屋|生信技能树    

GMT+8, 2019-6-20 20:17 , Processed in 0.029885 second(s), 25 queries .

Powered by Discuz! X3.2

© 2001-2013 Comsenz Inc.