如何使用 R 从 FASTA 文件中为 BED 文件中定义的每个间隔提取序列？答案

【问题标题】：How can I extract sequences from a FASTA file for each of the intervals defined in a BED file using R?如何使用 R 从 FASTA 文件中为 BED 文件中定义的每个间隔提取序列？
【发布时间】：2016-02-01 13:46:20
【问题描述】：

如何从 FASTA 文件中提取使用 R 在 BED 文件中定义的每个间隔的序列？使用的参考基因组是“Gallus gallus”，可以通过以下方式获得：

source("http://bioconductor.org/biocLite.R")
biocLite("BSgenome.Ggallus.UCSC.galGal4")
    library(BSgenome.Ggallus.UCSC.galGal4)

我的数据文件是 gRanges 包的结果

library("GenomicRanges")

> olaps
GRanges object with 2141 ranges and 0 metadata columns:
         seqnames               ranges strand
            <Rle>            <IRanges>  <Rle>
     [1]    chr14 [ 1665929,  1673673]      *
     [2]    chr14 [ 2587465,  2595209]      *
     [3]    chr14 [ 8143785,  8151529]      *
     [4]    chr14 [ 9779705,  9787449]      *
     [5]    chr14 [10281129, 10288873]      *
     ...      ...                  ...    ...
  [2137]    chr24   [3280553, 3288297]      *
  [2138]    chr24   [3330889, 3338633]      *
  [2139]    chr24   [3005641, 3015321]      *
  [2140]    chr24   [3319273, 3327017]      *
  [2141]    chr24   [5549545, 5557289]      *
  -------
  seqinfo: 31 sequences from an unspecified genome; no seqlengths

我可以在 data.table 中进行转换

olaps<- as.data.table(olaps)

要使用的示例：

olaps<-"seqnames    start      end width strand
chr1  1665929  1673673  7745      *
chr1  2587465  2595209  7745      *
chr1  8143785  8151529  7745      *
chr2  9779705  9787449  7745      *
chr2 10281129 10288873  7745      *"
olaps<-read.table(text=olaps,header=T)

预期结果：像这样（fasta 格式）：

>SEQUENCE_1
ACTGACTAGCATCGCAT...
>SEQUENCE_2
ACGTAGAGAGGGACATA...
>SEQUENCE_3...

我尝试使用这个包到现在都不成功：

source("http://bioconductor.org/biocLite.R")
biocLite("rtracklayer")

【问题讨论】：

我从未尝试将 R 用于 fasta 文件...但是快速的 Google 搜索会返回几个具有功能的网页。您是否尝试过 Biostrings 包中的 readFASTA（来自 Bioconductor）？还是来自seqinr 包的read.fasta？
seq = BSgenome::getSeq(BSgenome.Ggallus.UCSC.galGal4, olaps); Biostrings::writeXStringSet(seq, ...) 但在 Bioconductor 上询问 support forum。

标签： r range bioconductor

【解决方案1】：

这个，应该能解决你的把戏：

第一：

seq = BSgenome::getSeq(BSgenome.Ggallus.UCSC.galGal4, olaps)

为序列添加名称：

names(seq) = paste0("SEQUENCE_", seq_along(seq))

从您的序列中生成“.fasta”：

Biostrings::writeXStringSet(seq, "my.fasta")

【讨论】：