如何读取DNAstring的核苷酸序列

Question

我是新来的，我有一个组装基因组，我读了字符串集

x <- readDNAStringSet("GCA_000146045.2_R64_genomic.fna", format = "fasta")

然后我得到了这个

    A DNAStringSet instance of length 16
           width seq                                            enter code here   names               
    [1]  230218 CCACACCACACCCACACACCCAC...GGTGTGGGTGTGGTGTGTGTGGG BK006935.2 TPA_in...
    [2]  813184 AAATAGCCCTCATGTACGTCTCC...GTGGGTGTGGTGTGTGGGTGTGT BK006936.2 TPA_in...
    [3]  316620 CCCACACACCACACCCACACCAC...GTGTGGTGGGTGTGGTGTGTGTG BK006937.2 TPA_in...
    [4] 1531933 ACACCACACCCACACCACACCCA...TAAAGGTAGTAAGTAGCTTTTGG BK006938.2 TPA_in...
    [5]  576874 CGTCTCCTCCAAGCCCTGTTGTC...GTTTCATTTTCATTTTTTTTTTT BK006939.2
TPA_in...
...     ... ...
    [12] 1078177 CACACACACACACCACCCACACA...GGAGACGTACATGAGGGCTATTT BK006945.2 TPA_in...
    [13]  924431 CCACACACACACCACACCCACAC...GTGGGTGTGGTGTGTGTGTGGGG BK006946.2 TPA_in...
    [14]  784333 CCGGCTTTCTGACCGAAATTAAA...GTGTGTGTGGGTGTGGTGTGGGT BK006947.3
TPA_in...
     [15] 1091291 ACACCACACCCACACCACACCCA...GAGAGTGTGTGGGTGTGGTGTGT BK006948.2 TPA_in...
     [16]  948066 AAATAGCCCTCATGTACGTCTCC...TTTTTTTTAATTTCGGTCAGAAA BK006949.2 TPA_in...

下一个操作是对宽度进行排序。

 width(x)
[1]  230218  813184  316620 1531933  576874  270161 1090940  562643  439888
[10]  745751  666816 1078177  924431  784333 1091291  948066



 sort(width(x))
[1]  230218  270161  316620  439888  562643  576874  666816  745751  784333
[10]  813184  924431  948066 1078177 1090940 1091291 1531933

现在，我如何读取字符串745751的字母中的核苷酸序列？假设我hAVE一个更大的基因组，我只知道我感兴趣的字符串的长度，但我不知道它的确切位置，白化排序，我怎么能找到它呢？

Answer 1

你是说这个吗？

提取x中所有width为745751的width：

x[width(x) == 745751]

转换为character

as.character(x[width(x) == 745751])

下面是一个基于一些最小样本数据的示例

x <- DNAStringSet(x = c("AGTCTATCTA", "ACTAA", "ATGAA"))

提取所有长度5的序列：

x[width(x) == 5]
#A DNAStringSet instance of length 2
#    width seq
#[1]     5 ACTAA
#[2]     5 ATGAA

转换为character向量

as.character(x[width(x) == 5])
#[1] "ACTAA" "ATGAA"