如何使用更快的 R 脚本从 .sam 文件生成 wig 文件？答案

【问题标题】：How can I produce wig files from .sam files using a faster R script?如何使用更快的 R 脚本从 .sam 文件生成 wig 文件？
【发布时间】：2014-07-17 15:37:36
【问题描述】：

我有一个 R 脚本，在其中我可以在映射后从 .sam 文件中读取行，并且我想将 sam 文件的行解析为字符串，以便更轻松地操作它们并创建我想要的 wig 文件或计算我需要的 cov3 和 cov5。你能帮我让这个脚本更快地工作吗？如何更快地将巨大的 .sam 文件的行解析为数据帧？这是我的脚本：

gc()
rm(list=ls()) 

exptPath <- "/home/dimitris/INDEX3PerfectUnique31cov5.sam"


lines <- readLines(exptPath)
pos = lines
pos
chrom = lines
chrom
pos = ""
chrom = ""
nn = length(lines)
nn

# parse lines of sam file into strings(this part is very very slow)
rr = strsplit(lines,"\t", fixed = TRUE)
rr
trr = do.call(rbind.data.frame, rr)
pos = as.numeric(as.character(trr[8:nn,4]))
# for cov3
#pos = pos+25
#pos
chrom = trr[8:nn,3]
pos = as.numeric(pos)
pos

tab1 = table(chrom,pos, exclude="")
tab1

ftab1 = as.data.frame(tab1)
ftab1 = subset(ftab1, ftab1[3] != 0)
ftab1 = subset(ftab1, ftab1[1] != "<NA>")
oftab1 = ftab1[ order(ftab1[,1]), ]
final.ftab1 = oftab1[,2:3]
write.table(final.ftab1, "ind3_cov5_wig.txt", row.names=FALSE,
            sep="   ", quote=FALSE)

【问题讨论】：

标签： r bioinformatics rscript

【解决方案1】：

如果不访问示例输入和输出（例如，Dropbox 上的数据子集），就很难提供详细的答案。 Bioconductor 解决方案会将 sam 文件转换为 bam

library(Rsamtools)
bam <- "/path/to/new.bam")
asBam("/path/to/old.sam", bam)

然后读取数据，也许直接读取（参见?scanBam 和?ScanBamParam 仅导入感兴趣的字段/区域）

rr <- scanBam(bam)

或者说到底更方便

library(GenomicAlignments)
aln <- readGAlignments(bam)
## maybe cvg <- coverage(bam) ?

进行操作需要几个步骤，以 GRanges 对象（有点像 data.frame，但行具有基因组坐标）或相关对象结束

## ...???
## gr <- GRanges(seqnames, IRanges(start, end), strand=..., score=...)

最终目标是使用

导出到 wig / bigWig / bed 文件

library(rtracklayer)
export(gr, "/path/to.wig")

有大量的帮助资源，包括包小插曲、手册页和 Bioconductor mailing list

【讨论】：