如何從FASTA文件中爲使用R的BED文件中定義的每個間隔提取序列？

如何從FASTA文件中爲使用R的BED文件中定義的每個間隔提取序列？使用的參照基因組是「原雞」，可以通過以下步驟獲得：如何從FASTA文件中爲使用R的BED文件中定義的每個間隔提取序列？

source("http://bioconductor.org/biocLite.R") 
biocLite("BSgenome.Ggallus.UCSC.galGal4") 
    library(BSgenome.Ggallus.UCSC.galGal4)

我的數據文件是GRANGES包

library("GenomicRanges") 

> olaps 
GRanges object with 2141 ranges and 0 metadata columns: 
     seqnames    ranges strand 
      <Rle>   <IRanges> <Rle> 
    [1] chr14 [ 1665929, 1673673]  * 
    [2] chr14 [ 2587465, 2595209]  * 
    [3] chr14 [ 8143785, 8151529]  * 
    [4] chr14 [ 9779705, 9787449]  * 
    [5] chr14 [10281129, 10288873]  * 
    ...  ...     ... ... 
    [2137] chr24 [3280553, 3288297]  * 
    [2138] chr24 [3330889, 3338633]  * 
    [2139] chr24 [3005641, 3015321]  * 
    [2140] chr24 [3319273, 3327017]  * 
    [2141] chr24 [5549545, 5557289]  * 
    ------- 
    seqinfo: 31 sequences from an unspecified genome; no seqlengths

的結果我可以在data.table

變換要使用

olaps<- as.data.table(olaps)

實施例：

olaps<-"seqnames start  end width strand 
chr1 1665929 1673673 7745  * 
chr1 2587465 2595209 7745  * 
chr1 8143785 8151529 7745  * 
chr2 9779705 9787449 7745  * 
chr2 10281129 10288873 7745  *" 
olaps<-read.table(text=olaps,header=T)

預期成果：像這樣（FASTA格式）：

>SEQUENCE_1 
ACTGACTAGCATCGCAT... 
>SEQUENCE_2 
ACGTAGAGAGGGACATA... 
>SEQUENCE_3...

我試圖用這個包成功到現在爲止：

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

來源

2016-02-01 Fpertille

從來沒有試圖使用R的fasta文件...但快速谷歌搜索返回幾個功能的網頁。您是否嘗試過'Biostrings'軟件包（來自'Bioconductor'）的'readFASTA'？或者'seqinr'包中的'read.fasta'？ – Laterow

'seq = BSgenome :: getSeq（BSgenome.Ggallus.UCSC.galGal4，olaps）; Biostrings :: writeXStringSet（seq，...）'，但詢問Bioconductor [支持論壇]（https://support.bioconductor.org）。 –

這一點，應該可以解決你的絕招：

第一張：

seq = BSgenome::getSeq(BSgenome.Ggallus.UCSC.galGal4, olaps)

至名添加到序列：

names(seq) = paste0("SEQUENCE_", seq_along(seq))

爲了從您的序列的「.fasta」：分別提供

Biostrings::writeXStringSet(seq, "my.fasta")

來源

2016-08-25 16:19:24 Fpertille

如何從FASTA文件中爲使用R的BED文件中定義的每個間隔提取序列？

回答

相關問題