I am getting a very strange result for the plot of TSS enrichment (see attached image). I have tried other packages that do QC for ATACseq and they gave a much more expected enrichment plot so I am very confused as to what is going on. Any ideas? I have also attached my code and can send a minimal bam file to reproduce the result if desired.
Hello,
I am getting a very strange result for the plot of TSS enrichment (see attached image). I have tried other packages that do QC for ATACseq and they gave a much more expected enrichment plot so I am very confused as to what is going on. Any ideas? I have also attached my code and can send a minimal bam file to reproduce the result if desired.
Best, Alex.
`library(TxDb.Hsapiens.UCSC.hg38.knownGene) library(ATACseqQC)
bamTop100 <- scanBam(BamFile("rmdup.rmmm.sort2.bam", yieldSize = 100), param = ScanBamParam(tag=unlist(possibleTag)))[[1]]$tag tags <- names(bamTop100)[lengths(bamTop100)>0] tags which <- as(seqinformation[seqlev], "GRanges") gal <- readBamFile("rmdup.rmmm.sort2.bam", tag=tags, which=which, asMates=TRUE, bigFile=TRUE)
files will be output into outPath
outPath <- "splited" dir.create(outPath)
shift the coordinates of 5'ends of alignments in the bam file
seqinformation <- seqinfo(TxDb.Hsapiens.UCSC.hg38.knownGene) which <- as(seqinformation[seqlev], "GRanges") gal <- readBamFile("rmdup.rmmm.sort2.bam", tag=tags, which=which, asMates=TRUE, bigFile=TRUE) shiftedBamfile <- file.path(outPath, "shifted.bam") gal1 <- shiftGAlignmentsList(gal, outbam=shiftedBamfile) txs <- transcripts(TxDb.Hsapiens.UCSC.hg38.knownGene)
tsse <- TSSEscore(gal1, txs) tsse$TSSEscore
plot(100*(-9:10-.5), tsse$values, type="b", xlab="distance to TSS", ylab="aggregate TSS score") `