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

问题描述

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

我有R中枯草杆菌的完整DNA序列，用seqinr导入。我使用seqinr；

包的函数read.fast将DNA序列导入到R中

然后，我使用包"genbank"从枯草杆菌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序列)和"intergene"(带有坐标的Granges对象)

我在这个论坛上看到过有人问类似的问题，我试着用各种程序包来回答这些问题，但都没有成功，但我搞不懂。我希望有一个简单的解决方案；任何帮助都是我们的感激之情。

我想要的输出是生成一个包含以下格式的基因间序列的FASTA文件：

>intergenic_SEQ-BEGIN_dnaA
ATATATATATTATTTATTTTTTTTTTTTATTATAT
>intergenic_dnaA_dnaN
ATATATCGCGTCGATCTAGACTCAGGCATG
etc.

基本上是一行，其名称取自My Granges对象的intergenic_id name列中的基因间序列，后跟使用Grange中的坐标从FASTA中提取的序列。

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

推荐答案

我也认为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文件中提取序列片段的文章就介绍到这了，希望我们推荐的答案对大家有所帮助，也希望大家多多支持！