如何根据序列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 sequence或仅使用>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")

特别是快速概述.要将对象保存到相对于当前工作目录的"some"目录中的fasta文件"file.fa",请使用

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