如何存储作为for循环输出的DNA字符串实例？

Question

我正在使用库（biomaRt），库（Biostrings）和库（BiocManager）。我已经编写了一个代码来提取人类基因组的22号染色体，并使用IRanges函数在特定位置切割它。然后，使用getSeq（）函数，可以提取所需的序列。我已经写了一个for循环来做到这一点。但是，我的问题是我不知道如何正确保存for循环的输出。理想情况下，我想得到一个DNAStringSet，其中包含我使用for循环提取的所有DNAString。下面是我的代码：

HindIII <- DNAStringSet("AAGCTT")       
pdict0 <- PDict(HindIII) 
params <- new("BSParams", X = Hsapiens, FUN = matchPDict, simplify = TRUE, exclude = c("M", "_", 
"chrU", "alt", "random")) 
count_matches <- bsapply(params, pdict = pdict0)
cut_sites <- as.data.frame(count_matches$chr22)  
cut_sites <- cut_sites[ -c(1:2, 4:5)]

y <- data.frame()
for (i in cut_sites[1:nrow(cut_sites),]+1) {   
    for (y in cut_sites[2:nrow(cut_sites),]) {   
    fragments <- IRanges(start = i + 1, end = y) 
    fragments <- as.data.frame(fragments) 
     store_seqs <- getSeq(Hsapiens, "chr22", start = fragments[,1], 
             end = fragments[,2], strand = "+") #gets the sequence
     y <- rbind(y, store_seqs) #how do I store a DNA-string instance?
  }
}

谢谢你！

Answer 1

据我所知，您似乎有约8000个站点，并且您想要每个唯一组合的序列，即使在输入序列之前，这也为您提供了庞大的数据集，即select（8000,2）给出了约31996000组合，如果您选择更大的染色体，可能会出现问题。

下面，我提出一个解决方案，在其中您可以限制要查看的片段，然后从中获得所有序列。如果要忽略此限制，只需将其设置为染色体的长度即可，但要小心：

首先我们将网站移出：

library(BSgenome.Hsapiens.UCSC.hg38)
library(Biostrings)

frag_lim = 5000
gr = GRanges(seqnames="chr22",count_matches$chr22[[1]])

并且我们将其扩展到极限：

ext_gr = resize(gr,width=5000,fix="end")

基本上与扩展名重叠的任何地方都在另一个站点的5kb之内：

ovlp = as.data.frame(findOverlaps(gr,ext_gr))
#remove self overlap
ovlp = ovlp[ovlp[,1]!=ovlp[,2],]

现在我们有一个站点索引可以连接并获得其顺序：

head(ovlp)
   queryHits subjectHits
2          1           2
4          2           3
5          2           4
7          3           4
9          4           5
11         5           6

现在我们做：

fragment_gr = GRanges(seqnames=seqnames(ext_gr)[1],
IRanges(start=end(gr)[ovlp$queryHits]+1,end=end(gr)[ovlp$subjectHits]+1))

 summary(width(fragment_gr))
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
      7    1208    2483    2480    3749    5000

frag_seq = getSeq(Hsapiens,fragment_gr)

names(frag_seq) = paste0(as.character(seqnames(fragment_gr)),":",start(fragment_gr),"-",end(fragment_gr))

frag_seq
  A DNAStringSet instance of length 12175
        width seq           names               
    [1]  4273 TCCAG...GCTTA chr22:10512065-10...
    [2]   839 AAGGA...GCTTG chr22:10516337-10...
    [3]  2409 AAGGA...GCTTG chr22:10516337-10...
    [4]  1571 GAAGC...GCTTG chr22:10517175-10...
    [5]  4510 GAACT...GCTTG chr22:10518745-10...
    ...   ... ...
[12171]  4555 AGGCA...GCTTA chr22:50787152-50...
[12172]   388 CCTGT...GCTTA chr22:50791319-50...
[12173]  2527 CCTGT...GCTTT chr22:50791319-50...
[12174]  2140 AGATA...GCTTT chr22:50791706-50...
[12175]   364 GGCCA...GCTTT chr22:50798867-50...

关键是收集所需序列的所有范围并一次性执行getSeq，对此已进行了优化。