r - 如何根据序列 ID 或名称对 fasta 文件中的序列进行子集化？

Question

我在 fasta 文件中有一组 EST 序列。这里，如何根据序列ID或名称进行子集化？

>gi|296783888|gb|GW992815.1|GW992815 UAS-Mi10 Complementary DNA of mulberry (Morus indica) Morus indica cDNA 5' similar to Putative phosphoribosyltransferase/phosphoribosylanthranilate-like gene from Morus indica, mRNA sequence
GCAGCCGTCGGATCGTGAGCGTGATCGCGTGGCTAGTCGGGTTGGCGAAATGGTTGGATGATATCCGGAG
GTGGAGGAACCCCATTACCACGGTATTGGTCCACATCTTATATTTAGTGCTTGTTTGGTACCCGGATTTG
ATTGTCCCAACCGGGTTTTTATATGTGTTCCTAATCGGTGTATGGTACTATCGGTTTCGGCCCAAGATAC
CAGCGGGTATGGATACCCGACTCTCACAAGCTGAAGCGGTTGACCCGGATGAGCTTGATGAGGAATTCGA
CACCATACCGAGCTCAAAACCACCCGACATAATCAGGGTCCGGTATGACCGGTTGCGGATATTGGCAGCC
CGGGTTCAAACGGTTTTGGGTGATTTTGCAACACAAGGGGAGCGGGTTCAGGCCTTGGTTAGCTGGAGGG
ACCCAAGGGCCACAAAATTGTTCATAGGCGTGTGCTTGGCCATAACAATAATTCTCTATGTGGTGCCACC
CAAAATGGTTGCCGTGGCACTTGGATTCTACTATTTACGACACCCCATGTTCCGAGACCCCATGCCTCCT
GCAAGCTTGAATTTCTTCAGAAGGCTTCCAAGCCTTTCAGACCGCTTTAATGTAGATTAGAATATTATAT
GATTATTAGTAGGCCCAA

>gi|296783887|gb|GW992814.1|GW992814 UAS-Mi9 Complementary DNA of mulberry (Morus indica) Morus indica cDNA 5' similar to Dehydration-responsive protein RD22, Similar to BURP domain-containing protein like gene from Morus indica, mRNA sequence
AAGCAGTGGTCTAGAACCAGAGTGGCCCCTGCGATGCAGGTATCATCTCTATTATCAAAAGGGATAAGGG
GTGGATCCGTCGGGGATTTGAGTCTCACATGGTCGCTGATAACTTATTGAATGGATATTGGATTGTGTGC
AGTGCGACCTAAACAGGATTGCCGTTGGGGCCTGTGGTCAGAGATACCCCACACTTCTCAACTCCCAAAT
TGGATCTTGTTCCTTGTTTTCCTGTATTAAGCCTGACCCCTGAGGCTTTCGCCACTGCCAACTGGGTGCC
GCCTGCTGACTTCTGATTCCCCGTGCTAACGGTTACTCCCGATTCCTTATCCACATCGAAGATGAACTAT
TGACTTCCGCAAACTCAAAAGGCTGCAAGATATCACTGACCGCTGTCGGGATCCGCGATCGGCATATACG
CGAAATCCGATCCCGGATCCCGGGACTGCAGACGGCTGAA

如使用标题行 >gi|296783888|gb|GW992815.1|GW992815 UAS-Mi10 Complementary DNA of mulberry (Morus indica) Morus indica cDNA 5' similar to Putative phosphoribosyltransferase/phosphoribosylanthranilate-like gene from Morus indica , mRNA 序列 或仅使用 >gi|296783888
如何在 R 中执行此操作？

Answer 1

对于稍微重量级的解决方案，如果这适合 Bioconductor 工作流程，

source("http://bioconductor.org/biocLite.R")
biocLite("Rsamtools")

安装Biostrings和Rsamtools包，然后

library(Rsamtools)
indexFa("foo.fasta")   # create an index of file 'foo.fasta'
fa = FaFile("foo.fasta")  # reference the fasta file and it's index

您可以发现每个序列的坐标(名称和开始/结束)

gr = as(seqinfo(fa), "GRanges")

并通过选择适当的子集来查询序列中的任意序列和范围，例如，第二个序列，然后是示例中的第一个序列

getSeq(fa, gr[2:1])

或者通过部分匹配名称来查找坐标

idx = pmatch("gi|296783888", names(gr))  ## NA's if duplicates or not unique
seq = getSeq(fa, gr[idx])

“seq”是一个 DNAStringSet，可以通过多种方式进行操作；查看包中可用的小插曲

vignette(package="Biostrings")

尤其是 Quick Overview .要将对象保存到相对于当前工作目录的“some”目录中的 fasta 文件“file.fa”，请使用

writeXStringSet(seq, "some/file.fa")