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

Question

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

我有枯草芽孢杆菌中的R，与seqinr进口的完整DNA序列。 我使用包seqinr的函数read.fasta将dna序列导入到R中;

我然后创建来自枯草芽孢杆菌基因库文件中使用包“genbankr”来提取基因间坐标的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 
    ...  ...    ... ... .      ... 

因此，在RI有： “×”（完整的DNA与seqinr进口序列）和“基因间”（即GRANGES与对象坐标）

我看到有人问类似的问题，在这个论坛上，我曾尝试用多种套餐遵循这些答案没有成功，但我无法弄清楚。我希望有一个简单的解决方案，任何帮助表示赞赏。

我的期望的输出将产生具有以下列格式的基因间序列的FASTA文件：

>intergenic_SEQ-BEGIN_dnaA 
ATATATATATTATTTATTTTTTTTTTTTATTATAT 
>intergenic_dnaA_dnaN 
ATATATCGCGTCGATCTAGACTCAGGCATG 
etc. 

基本上与我的GRANGES对象从intergenic_id名称列采取的基因间序列的名称的行接着是使用GRanges中的坐标从fasta中提取的序列。

注：在所希望的输出I刚刚输入一些随机序列，正如例子。

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')