首页 > 解决方案 > 如何读取 DNAstring 的核苷酸序列

问题描述

我是 R 的新手。我有一个组装基因组,我阅读了字符串集

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 中字母的核苷酸序列?假设我有一个更大的基因组,并且我只知道我感兴趣的字符串的长度,但我不知道它的确切位置在没有排序的情况下,我怎么能找到它?

标签: rstring

解决方案


你是这个意思吗?

提取其中所有具有745751 的 s DNAStringxwidth

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"

推荐阅读