使用 GRanges 对象上的坐标从 FASTA 文件中提取序列片段

Question

我需要提取枯草芽孢杆菌的基因间序列。

我在 R 中有 B.subtilis 的完整 DNA 序列，用 seqinr 导入。我使用包 seqinr 的函数 read.fasta 将 dna 序列导入 R；

然后我使用包 "genbankr" 从枯草芽孢杆菌 genbank 文件创建了一个 GRanges 对象以提取基因间坐标。

这是我的带有基因间坐标的 GRanges 对象的格式：

GRanges object with 3841 ranges and 1 metadata column:
  seqnames             ranges strand |             intergenic_id
     <Rle>          <IRanges>  <Rle> |               <character>
  168      168       [   1,  409]      * | intergenic_SEQ-BEGIN_dnaA
  168      168       [1751, 1938]      * |      intergenic_dnaA_dnaN
  168      168       [3076, 3205]      * |      intergenic_dnaN_yaaA
  168      168       [3422, 3436]      * |      intergenic_yaaA_recF
  168      168       [4550, 4566]      * |      intergenic_recF_yaaB
  ...      ...                ...    ... .                       ...

因此在 R 中我有： "x"（用 seqinr 导入的完整 DNA 序列）和 "intergenic"（带坐标的 GRanges 对象）

我看到有人在这个论坛上问过类似的问题，我试图通过使用各种包来跟踪这些答案，但没有成功，但我无法弄清楚。我希望有一个简单的解决方案；感谢您的帮助。

我想要的输出是生成一个具有以下格式的基因间序列的 fasta 文件：

>intergenic_SEQ-BEGIN_dnaA
ATATATATATTATTTATTTTTTTTTTTTATTATAT
>intergenic_dnaA_dnaN
ATATATCGCGTCGATCTAGACTCAGGCATG
etc.

基本上一行包含取自我的 GRanges 对象的 intergenic_id 名称列的基因间序列名称，后跟使用 GRanges 中的坐标从 fasta 中提取的序列。

注意：在所需的输出中，我只是输入了一些随机序列，就像示例一样。

Answer 1

我也认为 BSGenome 是最好的方法（你可以构建你的 BSgenome），但你也可以使用 seqinr 创建你自己的函数：

require('seqinr')
require('GenomicRanges')

# Create objects
mygenome <- read.fasta('sequence.fasta')[[1]] # I assume it is just one chromosome
mygrs <- GRanges(seqnames=rep('NC_000964',3),
                 ranges=IRanges(c(1,50,100),c(20,55,103)),
                  strand=c('*','*','*'))
mcols(mygrs)$Gene <- c('GenA','GenB','GenC')
mygrs
# GRanges object with 3 ranges and 1 metadata column:
#   seqnames     ranges strand |        Gene
# <Rle>  <IRanges>  <Rle> | <character>
#   [1] NC_000964 [  1,  20]      * |        GenA
# [2] NC_000964 [ 50,  55]      * |        GenB

# Function to subset the seqinr list
 extractSeq <- function(x) {
    if (as.character(strand(x)) == '-') {
      comp(mygenome[end(x):start(x)])
    } else {
      mygenome[start(x):end(x)]
    }
  }    

# Ex
  extractSeq(mygrs[1])
  #  [1] "a" "t" "c" "t" "t" "t" "t" "t" "c" "g" "g" "c" "t" "t" "t" "t" "t" "t" "t" "a"    

# Apply to all
  myseqs <- lapply(mygrs, extractSeq)

# write to a file
  write.fasta(myseqs, mcols(mygrs)$Gene, 'myfile.fa')

使用 GRanges 对象上的坐标从 FASTA 文件中提取序列片段

Extract sequence fragments from FASTA file using coordinates on a GRanges object

r

bioinformatics