3

マッピング後に.samファイルから行を読み取ることができるRスクリプトがあり、samファイルの行を文字列に解析して、それらを操作しやすくし、必要なかつらファイルを作成したり、計算したりしたい.私が必要とする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)
4

1 に答える 1

1

サンプルの入力と出力 (ドロップボックス上のデータのサブセットなど) にアクセスせずに詳細な回答を提供することは困難です。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=...)

最終的な目標は、を使用してかつら/ bigWig /ベッドファイルにエクスポートすることです

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

パッケージ ビネット、マニュアル ページ、Bioconductorメーリング リストなど、豊富なヘルプ リソースがあります。

于 2014-05-28T12:46:35.233 に答える