我有一个问题，我基于“hg38.knownGene.gtf”对我的配对末端RNA-seq样本进行了修剪和比对，并且我得到了每个样本的.bam文件。当我在GenomicAlignments包中尝试使用此BAM文件进行读数计数时，输出文件为空。以下是我用于计数的代码：有人知道这个问题的原因吗？我使用了相同的参考“hg 38. knownGene. gtf”作为比对步骤。
- 我使用的是最新版本的桌面Rstudio.
input<- ("input.bam")
#total size of BAM file is around 2GB so, I used yieldSize to call .bam file.
data.FL <- BamFileList(input, yieldSize = 3000000)
# reference file
gtf.file<-"hg38.knownGene.gtf" #ucsc
TxDb<-makeTxDbFromGFF(gtf.file,format="gtf")
#GrangesList
exons<-exonsBy(TxDb,by="gene")
transcripts<-transcriptsBy(TxDb,by="gene")
grl<-GRangesList(c(exons,transcripts))
#counting
data<-summarizeOverlaps(grl,data.FL,mode = "Union",singleEnd=F,ignore.strand=TRUE)
#writing the output in .csv:
vector <- as.vector(mcols(data)$data)
write.csv(vector,"vector.csv")提前感谢您的帮助。
1条答案
按热度按时间ne5o7dgx1#
我将作一般性评论，而不涉及技术部分：
使用标准软件对基因、区间或区域的读取进行计数。这可以是类似featureCounts的东西，它既可以作为命令行的独立工具使用，也可以通过一个名为Rsubread的R包从R内部进行计数。不要为这样的标准软件定制方法。将整个BAM文件读取到内存中确实没有意义。featureCounts会做得更好，它占用内存更少，而且几乎不需要任何时间。