获取限制位点附近300碱基的序列

Question

我有许多包含序列读取的sam文件，我想得到在两个方向上至少300个碱基的限制序列附近存在的读取。
包含限制位置的第一个文件，有两列。
chr01 4957
chr01 6605
chr02 19968
chr02 21055
chr02 208555
chr03 243398号

以SAM文件格式读取的第二个文件。（近2.6米线）

id1995 147 chr03 119509969 42 85M
id1999 83 chr10 131560619 26 81M
id1999 163 chr10 131560429 26 85M
id2099 73 chr10 60627850 42 81M

现在我想比较sam文件的第3列和position文件的第1列，sam文件的第4列和position文件的第2列。 我试着用R语言来做，但是因为数据很大，所以要花很多时间去做。 如果可以通过实现最佳算法来提高R脚本的执行速度。 R代码：

pos = read.csv(file="sites.csv",header=F,sep="\t")
fastq = read.csv(file="reads.sam", header=F,sep="\t")
newFastq = data.frame(fastq)
newFastq = NULL
trim <- function (x) gsub("^\\s+|\\s+$", "", x)

for(i  in 1:nrow(fastq)){
for(j in 1:nrow(pos)){
    if(as.character(fastq[i,3]) == trim(as.character(pos[j,1]))){
      if(fastq[i,4] - pos[j,2] < 300 && fastq[i,4] - pos[j,2] > -300){
       newFastq = rbind(newFastq,fastq[i,])
     }
   }
 }  
}
#Write data into file 
write.table(newFastq, file = "sitesFound.csv",row.names=FALSE, na="",quote=FALSE,col.names=FALSE, sep="\t")

或者可以通过用perl编写代码来改进这段代码。你知道吗

Answer 1

一种总体策略是使用BioconductorRsamtoolsasBam()和indexBam()生成索引的bam文件。将第一个文件读入数据框并构造一个GenomicRangesGRanges()对象。最后，使用GenomicAlignmentsreadGAlignments()读取bam文件，使用GRanges()作为ScanBamParam()的which=参数。如果您决定走这条路线，那么Bioconductor支持站点https://support.bioconductor.org更适合回答Bioconductor问题。你知道吗

看起来您希望读取的值在GRanges对象的+/-300碱基对内。调整GRanges的大小

library(GenomicRanges)
## create gr = GRanges(...)
gr = resize(gr, width = 600, fix="center")

将其用作ScanBamParam()中的which=，然后读取BAM文件

library(GenomicAlignments)
param = ScanBamParam(which = gr)
reads = readGAlignments("your.bam", param = param)

使用what=控制从BAM文件读取的字段，例如

param = ScanBamParam(which = gr, what = "seq")