将DNAStringSet转换为R中的元素列表？(seq[[1]][["seq"]]中出错：R中下标出界)

Question

我有一个床文件，其中包含DNA序列信息如下。

**

track name="194" description="194 methylation (sites)" color=0,60,120 useScore=1
chr1    15864   15866   FALSE   894 +
chr1    534241  534243  FALSE   921 -
chr1    710096  710098  FALSE   729 +
chr1    714176  714178  FALSE   12  -
chr1    720864  720866  FALSE   988 -

**

我在R中加载了床文件，并将矩阵命名为 数据集.我使用下面的代码来获取序列。

mydataSet_Test1<-dataSet[,1:3]

library(BSgenome.Hsapiens.UCSC.hg19)
genome <- BSgenome.Hsapiens.UCSC.hg19
chr<-as.matrix(as.character(mydataSet_Test1[,1]))

#50
start<-as.matrix(as.integer(as.character(mydataSet_Test1[,2]))-50)
end<-as.matrix(as.integer(as.character(mydataSet_Test1[,3]))+50)
Seqs50_Test1<-getSeq(genome,chr,start=start,end=end)

现在，Seqs50_Test1是： 大DNAStringSet.我现在要加载的是 BioSeqClass包 在R中，做一个 同源还原 在我的序列中，我想使用hr()函数。

我想使用hr()函数，根据软件包手册，它是这样的。

描述

通过序列相似性.hr(seq, method, identity, cdhit.path)过滤同源序列。

论据

序号 一个列表，每个蛋白酶序列有一个元素。这些元素分为两部分，一部分是描述（"desc"），另一部分是生物序列的字符串（"seq"）。

特性 一个介于0到1之间的数值，它被用作输入序列中的最大标识截止值。

方法 同源重结方法的字符串。这必须是 "cdhit "或 "aligndis "中的一个字符串。

我的问题是，如何将我的DNAStringSet转换为函数hr()想要的元素列表？我试着使用list()函数，但当我运行hr()函数时，它给了我一个错误提示 seq[]中出现错误。1][[["seq"]]：下标出界

完整的代码。

mydataSet = dataSet[,1:3]

library(BSgenome.Hsapiens.UCSC.hg19)
genome = BSgenome.Hsapiens.UCSC.hg19
chr = as.matrix(as.character(mydataSet[,1]))
start = as.matrix(as.integer(as.character(mydataSet[,2]))-200)
end = as.matrix(as.integer(as.character(mydataSet[,3]))+200)
Seqs = getSeq(genome,chr,start=start,end=end)
writeXStringSet(Seqs, "C:\\Users\\JL009\\Desktop\\Seqs.fasta", append=FALSE, format = "fasta")

#if (!requireNamespace("BiocManager", quietly = TRUE))
#  install.packages("BiocManager")

#BiocManager::install("BioSeqClass")

library(BioSeqClass)
library(Biostrings)
seq = as.character(readAAStringSet("C:\\Users\\JL009\\Desktop\\Seqs.fasta"))
reducSeqs = hr(seq, method="aligndis", identity=0.4)

Answer 1

那就试试这样吧。

library(BioSeqClass)
library(Biostrings)
library(BSgenome.Hsapiens.UCSC.hg19)

gr = GRanges(seqnames="chr1",IRanges(start=seq(10e6,11e6,length.out=10),width=50))
S = getSeq(BSgenome.Hsapiens.UCSC.hg19,gr)
names(S) = paste("seq",1:length(S))

input = lapply(seq_along(S),function(i)list(desc=names(S)[i],seq=as.character(S[[i]])))

hr(input,method="aligndis",identity=0.5)